libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/19_legendre/check_value.cc

   1 // 2007-02-04  Edward Smith-Rowland <3dw4rd@verizon.net>
   2 //
   3 // Copyright (C) 2007-2014 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>
  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>
 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>
 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>
 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>
 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>
 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>
 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>
 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
 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 }