1 // 2007-02-04  Edward Smith-Rowland <3dw4rd@verizon.net>
2 //
3 // Copyright (C) 2007-2016 Free Software Foundation, Inc.
4 //
5 // This file is part of the GNU ISO C++ Library.  This library is free
6 // software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 3, or (at your option)
9 // any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License along
17 // with this library; see the file COPYING3.  If not see
18 // <http://www.gnu.org/licenses/>.
19 
20 //  legendre
21 
22 
23 //  Compare against values generated by the GNU Scientific Library.
24 //  The GSL can be found on the web: http://www.gnu.org/software/gsl/
25 
26 #include <tr1/cmath>
27 #if defined(__TEST_DEBUG)
28 #include <iostream>
29 #define VERIFY(A) \
30 if (!(A)) \
31   { \
32     std::cout << "line " << __LINE__ \
33       << "  max_abs_frac = " << max_abs_frac \
34       << std::endl; \
35   }
36 #else
37 #include <testsuite_hooks.h>
38 #endif
39 #include "../testcase.h"
40 
41 
42 // Test data for l=0.
43 testcase_legendre<double> data001[] = {
44   { 1.0000000000000000, 0, -1.0000000000000000 },
45   { 1.0000000000000000, 0, -0.90000000000000002 },
46   { 1.0000000000000000, 0, -0.80000000000000004 },
47   { 1.0000000000000000, 0, -0.69999999999999996 },
48   { 1.0000000000000000, 0, -0.59999999999999998 },
49   { 1.0000000000000000, 0, -0.50000000000000000 },
50   { 1.0000000000000000, 0, -0.40000000000000002 },
51   { 1.0000000000000000, 0, -0.30000000000000004 },
52   { 1.0000000000000000, 0, -0.19999999999999996 },
53   { 1.0000000000000000, 0, -0.099999999999999978 },
54   { 1.0000000000000000, 0, 0.0000000000000000 },
55   { 1.0000000000000000, 0, 0.10000000000000009 },
56   { 1.0000000000000000, 0, 0.19999999999999996 },
57   { 1.0000000000000000, 0, 0.30000000000000004 },
58   { 1.0000000000000000, 0, 0.39999999999999991 },
59   { 1.0000000000000000, 0, 0.50000000000000000 },
60   { 1.0000000000000000, 0, 0.60000000000000009 },
61   { 1.0000000000000000, 0, 0.69999999999999996 },
62   { 1.0000000000000000, 0, 0.80000000000000004 },
63   { 1.0000000000000000, 0, 0.89999999999999991 },
64   { 1.0000000000000000, 0, 1.0000000000000000 },
65 };
66 
67 // Test function for l=0.
68 template <typename Tp>
test001()69 void test001()
70 {
71   const Tp eps = std::numeric_limits<Tp>::epsilon();
72   Tp max_abs_diff = -Tp(1);
73   Tp max_abs_frac = -Tp(1);
74   unsigned int num_datum = sizeof(data001)
75                          / sizeof(testcase_legendre<double>);
76   for (unsigned int i = 0; i < num_datum; ++i)
77     {
78       const Tp f = std::tr1::legendre(Tp(data001[i].l), Tp(data001[i].x));
79       const Tp f0 = data001[i].f0;
80       const Tp diff = f - f0;
81       if (std::abs(diff) > max_abs_diff)
82         max_abs_diff = std::abs(diff);
83       if (std::abs(f0) > Tp(10) * eps
84        && std::abs(f) > Tp(10) * eps)
85         {
86           const Tp frac = diff / f0;
87           if (std::abs(frac) > max_abs_frac)
88             max_abs_frac = std::abs(frac);
89         }
90     }
91   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
92 }
93 
94 // Test data for l=1.
95 testcase_legendre<double> data002[] = {
96   { -1.0000000000000000, 1, -1.0000000000000000 },
97   { -0.90000000000000002, 1, -0.90000000000000002 },
98   { -0.80000000000000004, 1, -0.80000000000000004 },
99   { -0.69999999999999996, 1, -0.69999999999999996 },
100   { -0.59999999999999998, 1, -0.59999999999999998 },
101   { -0.50000000000000000, 1, -0.50000000000000000 },
102   { -0.40000000000000002, 1, -0.40000000000000002 },
103   { -0.30000000000000004, 1, -0.30000000000000004 },
104   { -0.19999999999999996, 1, -0.19999999999999996 },
105   { -0.099999999999999978, 1, -0.099999999999999978 },
106   { 0.0000000000000000, 1, 0.0000000000000000 },
107   { 0.10000000000000009, 1, 0.10000000000000009 },
108   { 0.19999999999999996, 1, 0.19999999999999996 },
109   { 0.30000000000000004, 1, 0.30000000000000004 },
110   { 0.39999999999999991, 1, 0.39999999999999991 },
111   { 0.50000000000000000, 1, 0.50000000000000000 },
112   { 0.60000000000000009, 1, 0.60000000000000009 },
113   { 0.69999999999999996, 1, 0.69999999999999996 },
114   { 0.80000000000000004, 1, 0.80000000000000004 },
115   { 0.89999999999999991, 1, 0.89999999999999991 },
116   { 1.0000000000000000, 1, 1.0000000000000000 },
117 };
118 
119 // Test function for l=1.
120 template <typename Tp>
test002()121 void test002()
122 {
123   const Tp eps = std::numeric_limits<Tp>::epsilon();
124   Tp max_abs_diff = -Tp(1);
125   Tp max_abs_frac = -Tp(1);
126   unsigned int num_datum = sizeof(data002)
127                          / sizeof(testcase_legendre<double>);
128   for (unsigned int i = 0; i < num_datum; ++i)
129     {
130       const Tp f = std::tr1::legendre(Tp(data002[i].l), Tp(data002[i].x));
131       const Tp f0 = data002[i].f0;
132       const Tp diff = f - f0;
133       if (std::abs(diff) > max_abs_diff)
134         max_abs_diff = std::abs(diff);
135       if (std::abs(f0) > Tp(10) * eps
136        && std::abs(f) > Tp(10) * eps)
137         {
138           const Tp frac = diff / f0;
139           if (std::abs(frac) > max_abs_frac)
140             max_abs_frac = std::abs(frac);
141         }
142     }
143   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
144 }
145 
146 // Test data for l=2.
147 testcase_legendre<double> data003[] = {
148   { 1.0000000000000000, 2, -1.0000000000000000 },
149   { 0.71500000000000008, 2, -0.90000000000000002 },
150   { 0.46000000000000019, 2, -0.80000000000000004 },
151   { 0.23499999999999988, 2, -0.69999999999999996 },
152   { 0.039999999999999925, 2, -0.59999999999999998 },
153   { -0.12500000000000000, 2, -0.50000000000000000 },
154   { -0.25999999999999995, 2, -0.40000000000000002 },
155   { -0.36499999999999999, 2, -0.30000000000000004 },
156   { -0.44000000000000006, 2, -0.19999999999999996 },
157   { -0.48499999999999999, 2, -0.099999999999999978 },
158   { -0.50000000000000000, 2, 0.0000000000000000 },
159   { -0.48499999999999999, 2, 0.10000000000000009 },
160   { -0.44000000000000006, 2, 0.19999999999999996 },
161   { -0.36499999999999999, 2, 0.30000000000000004 },
162   { -0.26000000000000012, 2, 0.39999999999999991 },
163   { -0.12500000000000000, 2, 0.50000000000000000 },
164   { 0.040000000000000147, 2, 0.60000000000000009 },
165   { 0.23499999999999988, 2, 0.69999999999999996 },
166   { 0.46000000000000019, 2, 0.80000000000000004 },
167   { 0.71499999999999986, 2, 0.89999999999999991 },
168   { 1.0000000000000000, 2, 1.0000000000000000 },
169 };
170 
171 // Test function for l=2.
172 template <typename Tp>
test003()173 void test003()
174 {
175   const Tp eps = std::numeric_limits<Tp>::epsilon();
176   Tp max_abs_diff = -Tp(1);
177   Tp max_abs_frac = -Tp(1);
178   unsigned int num_datum = sizeof(data003)
179                          / sizeof(testcase_legendre<double>);
180   for (unsigned int i = 0; i < num_datum; ++i)
181     {
182       const Tp f = std::tr1::legendre(Tp(data003[i].l), Tp(data003[i].x));
183       const Tp f0 = data003[i].f0;
184       const Tp diff = f - f0;
185       if (std::abs(diff) > max_abs_diff)
186         max_abs_diff = std::abs(diff);
187       if (std::abs(f0) > Tp(10) * eps
188        && std::abs(f) > Tp(10) * eps)
189         {
190           const Tp frac = diff / f0;
191           if (std::abs(frac) > max_abs_frac)
192             max_abs_frac = std::abs(frac);
193         }
194     }
195   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
196 }
197 
198 // Test data for l=5.
199 testcase_legendre<double> data004[] = {
200   { -1.0000000000000000, 5, -1.0000000000000000 },
201   { 0.041141249999999997, 5, -0.90000000000000002 },
202   { 0.39951999999999993, 5, -0.80000000000000004 },
203   { 0.36519874999999991, 5, -0.69999999999999996 },
204   { 0.15263999999999994, 5, -0.59999999999999998 },
205   { -0.089843750000000000, 5, -0.50000000000000000 },
206   { -0.27063999999999988, 5, -0.40000000000000002 },
207   { -0.34538624999999995, 5, -0.30000000000000004 },
208   { -0.30751999999999996, 5, -0.19999999999999996 },
209   { -0.17882874999999995, 5, -0.099999999999999978 },
210   { 0.0000000000000000, 5, 0.0000000000000000 },
211   { 0.17882875000000015, 5, 0.10000000000000009 },
212   { 0.30751999999999996, 5, 0.19999999999999996 },
213   { 0.34538624999999995, 5, 0.30000000000000004 },
214   { 0.27064000000000010, 5, 0.39999999999999991 },
215   { 0.089843750000000000, 5, 0.50000000000000000 },
216   { -0.15264000000000022, 5, 0.60000000000000009 },
217   { -0.36519874999999991, 5, 0.69999999999999996 },
218   { -0.39951999999999993, 5, 0.80000000000000004 },
219   { -0.041141250000000407, 5, 0.89999999999999991 },
220   { 1.0000000000000000, 5, 1.0000000000000000 },
221 };
222 
223 // Test function for l=5.
224 template <typename Tp>
test004()225 void test004()
226 {
227   const Tp eps = std::numeric_limits<Tp>::epsilon();
228   Tp max_abs_diff = -Tp(1);
229   Tp max_abs_frac = -Tp(1);
230   unsigned int num_datum = sizeof(data004)
231                          / sizeof(testcase_legendre<double>);
232   for (unsigned int i = 0; i < num_datum; ++i)
233     {
234       const Tp f = std::tr1::legendre(Tp(data004[i].l), Tp(data004[i].x));
235       const Tp f0 = data004[i].f0;
236       const Tp diff = f - f0;
237       if (std::abs(diff) > max_abs_diff)
238         max_abs_diff = std::abs(diff);
239       if (std::abs(f0) > Tp(10) * eps
240        && std::abs(f) > Tp(10) * eps)
241         {
242           const Tp frac = diff / f0;
243           if (std::abs(frac) > max_abs_frac)
244             max_abs_frac = std::abs(frac);
245         }
246     }
247   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
248 }
249 
250 // Test data for l=10.
251 testcase_legendre<double> data005[] = {
252   { 1.0000000000000000, 10, -1.0000000000000000 },
253   { -0.26314561785585977, 10, -0.90000000000000002 },
254   { 0.30052979559999998, 10, -0.80000000000000004 },
255   { 0.085805795531640333, 10, -0.69999999999999996 },
256   { -0.24366274560000006, 10, -0.59999999999999998 },
257   { -0.18822860717773438, 10, -0.50000000000000000 },
258   { 0.096839064399999869, 10, -0.40000000000000002 },
259   { 0.25147634951601561, 10, -0.30000000000000004 },
260   { 0.12907202559999989, 10, -0.19999999999999996 },
261   { -0.12212499738710947, 10, -0.099999999999999978 },
262   { -0.24609375000000000, 10, 0.0000000000000000 },
263   { -0.12212499738710922, 10, 0.10000000000000009 },
264   { 0.12907202559999989, 10, 0.19999999999999996 },
265   { 0.25147634951601561, 10, 0.30000000000000004 },
266   { 0.096839064400000258, 10, 0.39999999999999991 },
267   { -0.18822860717773438, 10, 0.50000000000000000 },
268   { -0.24366274559999987, 10, 0.60000000000000009 },
269   { 0.085805795531640333, 10, 0.69999999999999996 },
270   { 0.30052979559999998, 10, 0.80000000000000004 },
271   { -0.26314561785585888, 10, 0.89999999999999991 },
272   { 1.0000000000000000, 10, 1.0000000000000000 },
273 };
274 
275 // Test function for l=10.
276 template <typename Tp>
test005()277 void test005()
278 {
279   const Tp eps = std::numeric_limits<Tp>::epsilon();
280   Tp max_abs_diff = -Tp(1);
281   Tp max_abs_frac = -Tp(1);
282   unsigned int num_datum = sizeof(data005)
283                          / sizeof(testcase_legendre<double>);
284   for (unsigned int i = 0; i < num_datum; ++i)
285     {
286       const Tp f = std::tr1::legendre(Tp(data005[i].l), Tp(data005[i].x));
287       const Tp f0 = data005[i].f0;
288       const Tp diff = f - f0;
289       if (std::abs(diff) > max_abs_diff)
290         max_abs_diff = std::abs(diff);
291       if (std::abs(f0) > Tp(10) * eps
292        && std::abs(f) > Tp(10) * eps)
293         {
294           const Tp frac = diff / f0;
295           if (std::abs(frac) > max_abs_frac)
296             max_abs_frac = std::abs(frac);
297         }
298     }
299   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
300 }
301 
302 // Test data for l=20.
303 testcase_legendre<double> data006[] = {
304   { 1.0000000000000000, 20, -1.0000000000000000 },
305   { -0.14930823530984821, 20, -0.90000000000000002 },
306   { 0.22420460541741344, 20, -0.80000000000000004 },
307   { -0.20457394463834172, 20, -0.69999999999999996 },
308   { 0.15916752910098114, 20, -0.59999999999999998 },
309   { -0.048358381067373557, 20, -0.50000000000000000 },
310   { -0.10159261558628156, 20, -0.40000000000000002 },
311   { 0.18028715947998047, 20, -0.30000000000000004 },
312   { -0.098042194344594741, 20, -0.19999999999999996 },
313   { -0.082077130944527649, 20, -0.099999999999999978 },
314   { 0.17619705200195312, 20, 0.0000000000000000 },
315   { -0.082077130944528037, 20, 0.10000000000000009 },
316   { -0.098042194344594741, 20, 0.19999999999999996 },
317   { 0.18028715947998047, 20, 0.30000000000000004 },
318   { -0.10159261558628112, 20, 0.39999999999999991 },
319   { -0.048358381067373557, 20, 0.50000000000000000 },
320   { 0.15916752910098084, 20, 0.60000000000000009 },
321   { -0.20457394463834172, 20, 0.69999999999999996 },
322   { 0.22420460541741344, 20, 0.80000000000000004 },
323   { -0.14930823530984949, 20, 0.89999999999999991 },
324   { 1.0000000000000000, 20, 1.0000000000000000 },
325 };
326 
327 // Test function for l=20.
328 template <typename Tp>
test006()329 void test006()
330 {
331   const Tp eps = std::numeric_limits<Tp>::epsilon();
332   Tp max_abs_diff = -Tp(1);
333   Tp max_abs_frac = -Tp(1);
334   unsigned int num_datum = sizeof(data006)
335                          / sizeof(testcase_legendre<double>);
336   for (unsigned int i = 0; i < num_datum; ++i)
337     {
338       const Tp f = std::tr1::legendre(Tp(data006[i].l), Tp(data006[i].x));
339       const Tp f0 = data006[i].f0;
340       const Tp diff = f - f0;
341       if (std::abs(diff) > max_abs_diff)
342         max_abs_diff = std::abs(diff);
343       if (std::abs(f0) > Tp(10) * eps
344        && std::abs(f) > Tp(10) * eps)
345         {
346           const Tp frac = diff / f0;
347           if (std::abs(frac) > max_abs_frac)
348             max_abs_frac = std::abs(frac);
349         }
350     }
351   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
352 }
353 
354 // Test data for l=50.
355 testcase_legendre<double> data007[] = {
356   { 1.0000000000000000, 50, -1.0000000000000000 },
357   { -0.17003765994383685, 50, -0.90000000000000002 },
358   { 0.13879737345093118, 50, -0.80000000000000004 },
359   { -0.014572731645892805, 50, -0.69999999999999996 },
360   { -0.058860798844002173, 50, -0.59999999999999998 },
361   { -0.031059099239609828, 50, -0.50000000000000000 },
362   { 0.041569033381825368, 50, -0.40000000000000002 },
363   { 0.10911051574714808, 50, -0.30000000000000004 },
364   { 0.083432272204197466, 50, -0.19999999999999996 },
365   { -0.038205812661313579, 50, -0.099999999999999978 },
366   { -0.11227517265921705, 50, 0.0000000000000000 },
367   { -0.038205812661314169, 50, 0.10000000000000009 },
368   { 0.083432272204197466, 50, 0.19999999999999996 },
369   { 0.10911051574714808, 50, 0.30000000000000004 },
370   { 0.041569033381824647, 50, 0.39999999999999991 },
371   { -0.031059099239609828, 50, 0.50000000000000000 },
372   { -0.058860798844001430, 50, 0.60000000000000009 },
373   { -0.014572731645892805, 50, 0.69999999999999996 },
374   { 0.13879737345093118, 50, 0.80000000000000004 },
375   { -0.17003765994383663, 50, 0.89999999999999991 },
376   { 1.0000000000000000, 50, 1.0000000000000000 },
377 };
378 
379 // Test function for l=50.
380 template <typename Tp>
test007()381 void test007()
382 {
383   const Tp eps = std::numeric_limits<Tp>::epsilon();
384   Tp max_abs_diff = -Tp(1);
385   Tp max_abs_frac = -Tp(1);
386   unsigned int num_datum = sizeof(data007)
387                          / sizeof(testcase_legendre<double>);
388   for (unsigned int i = 0; i < num_datum; ++i)
389     {
390       const Tp f = std::tr1::legendre(Tp(data007[i].l), Tp(data007[i].x));
391       const Tp f0 = data007[i].f0;
392       const Tp diff = f - f0;
393       if (std::abs(diff) > max_abs_diff)
394         max_abs_diff = std::abs(diff);
395       if (std::abs(f0) > Tp(10) * eps
396        && std::abs(f) > Tp(10) * eps)
397         {
398           const Tp frac = diff / f0;
399           if (std::abs(frac) > max_abs_frac)
400             max_abs_frac = std::abs(frac);
401         }
402     }
403   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
404 }
405 
406 // Test data for l=100.
407 testcase_legendre<double> data008[] = {
408   { 1.0000000000000000, 100, -1.0000000000000000 },
409   { 0.10226582055871908, 100, -0.90000000000000002 },
410   { 0.050861167913584124, 100, -0.80000000000000004 },
411   { -0.077132507199778780, 100, -0.69999999999999996 },
412   { -0.023747023905133110, 100, -0.59999999999999998 },
413   { -0.060518025961861198, 100, -0.50000000000000000 },
414   { -0.072258202125684429, 100, -0.40000000000000002 },
415   { 0.057127392202801719, 100, -0.30000000000000004 },
416   { 0.014681835355659636, 100, -0.19999999999999996 },
417   { -0.063895098434750303, 100, -0.099999999999999978 },
418   { 0.079589237387178727, 100, 0.0000000000000000 },
419   { -0.063895098434749775, 100, 0.10000000000000009 },
420   { 0.014681835355659636, 100, 0.19999999999999996 },
421   { 0.057127392202801719, 100, 0.30000000000000004 },
422   { -0.072258202125685012, 100, 0.39999999999999991 },
423   { -0.060518025961861198, 100, 0.50000000000000000 },
424   { -0.023747023905134217, 100, 0.60000000000000009 },
425   { -0.077132507199778780, 100, 0.69999999999999996 },
426   { 0.050861167913584124, 100, 0.80000000000000004 },
427   { 0.10226582055871723, 100, 0.89999999999999991 },
428   { 1.0000000000000000, 100, 1.0000000000000000 },
429 };
430 
431 // Test function for l=100.
432 template <typename Tp>
test008()433 void test008()
434 {
435   const Tp eps = std::numeric_limits<Tp>::epsilon();
436   Tp max_abs_diff = -Tp(1);
437   Tp max_abs_frac = -Tp(1);
438   unsigned int num_datum = sizeof(data008)
439                          / sizeof(testcase_legendre<double>);
440   for (unsigned int i = 0; i < num_datum; ++i)
441     {
442       const Tp f = std::tr1::legendre(Tp(data008[i].l), Tp(data008[i].x));
443       const Tp f0 = data008[i].f0;
444       const Tp diff = f - f0;
445       if (std::abs(diff) > max_abs_diff)
446         max_abs_diff = std::abs(diff);
447       if (std::abs(f0) > Tp(10) * eps
448        && std::abs(f) > Tp(10) * eps)
449         {
450           const Tp frac = diff / f0;
451           if (std::abs(frac) > max_abs_frac)
452             max_abs_frac = std::abs(frac);
453         }
454     }
455   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
456 }
457 
main(int,char **)458 int main(int, char**)
459 {
460   test001<double>();
461   test002<double>();
462   test003<double>();
463   test004<double>();
464   test005<double>();
465   test006<double>();
466   test007<double>();
467   test008<double>();
468   return 0;
469 }
470