libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/08_cyl_bessel_i/check_value.cc

   1 // 2007-02-04  Edward Smith-Rowland <3dw4rd@verizon.net>
   2 //
   3 // Copyright (C) 2007, 2009 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 //  cyl_bessel_i
  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 nu=0.0000000000000000.
  43 testcase_cyl_bessel_i<double> data001[] = {
  44   { 1.0000000000000000, 0.0000000000000000, 0.0000000000000000 },
  45   { 27.239871823604439, 0.0000000000000000, 5.0000000000000000 },
  46   { 2815.7166284662558, 0.0000000000000000, 10.000000000000000 },
  47   { 339649.37329791381, 0.0000000000000000, 15.000000000000000 },
  48   { 43558282.559553474, 0.0000000000000000, 20.000000000000000 },
  49   { 5774560606.4663124, 0.0000000000000000, 25.000000000000000 },
  50   { 781672297823.97925, 0.0000000000000000, 30.000000000000000 },
  51   { 107338818494514.42, 0.0000000000000000, 35.000000000000000 },
  52   { 14894774793419918., 0.0000000000000000, 40.000000000000000 },
  53   { 2.0834140751773158e+18, 0.0000000000000000, 45.000000000000000 },
  54   { 2.9325537838493463e+20, 0.0000000000000000, 50.000000000000000 },
  55   { 4.1487895607332160e+22, 0.0000000000000000, 55.000000000000000 },
  56   { 5.8940770556098216e+24, 0.0000000000000000, 60.000000000000000 },
  57   { 8.4030398456255610e+26, 0.0000000000000000, 65.000000000000000 },
  58   { 1.2015889579125422e+29, 0.0000000000000000, 70.000000000000000 },
  59   { 1.7226390780357971e+31, 0.0000000000000000, 75.000000000000000 },
  60   { 2.4751784043341670e+33, 0.0000000000000000, 80.000000000000000 },
  61   { 3.5634776304081418e+35, 0.0000000000000000, 85.000000000000000 },
  62   { 5.1392383455086475e+37, 0.0000000000000000, 90.000000000000000 },
  63   { 7.4233258618752096e+39, 0.0000000000000000, 95.000000000000000 },
  64   { 1.0737517071310986e+42, 0.0000000000000000, 100.00000000000000 },
  65 };
  66
  67 // Test function for nu=0.0000000000000000.
  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_cyl_bessel_i<double>);
  76   for (unsigned int i = 0; i < num_datum; ++i)
  77     {
  78       const Tp f = std::tr1::cyl_bessel_i(Tp(data001[i].nu), 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.5000000000000015e-12));
  92 }
  93
  94 // Test data for nu=0.33333333333333331.
  95 testcase_cyl_bessel_i<double> data002[] = {
  96   { 0.0000000000000000, 0.33333333333333331, 0.0000000000000000 },
  97   { 26.897553069268362, 0.33333333333333331, 5.0000000000000000 },
  98   { 2799.2396097056790, 0.33333333333333331, 10.000000000000000 },
  99   { 338348.63146593655, 0.33333333333333331, 15.000000000000000 },
 100   { 43434263.927938439, 0.33333333333333331, 20.000000000000000 },
 101   { 5761474759.6213636, 0.33333333333333331, 25.000000000000000 },
 102   { 780201111830.30225, 0.33333333333333331, 30.000000000000000 },
 103   { 107166066959051.92, 0.33333333333333331, 35.000000000000000 },
 104   { 14873836574083760., 0.33333333333333331, 40.000000000000000 },
 105   { 2.0808143020217085e+18, 0.33333333333333331, 45.000000000000000 },
 106   { 2.9292639365644229e+20, 0.33333333333333331, 50.000000000000000 },
 107   { 4.1445621624120489e+22, 0.33333333333333331, 55.000000000000000 },
 108   { 5.8885758374365916e+24, 0.33333333333333331, 60.000000000000000 },
 109   { 8.3958047021083955e+26, 0.33333333333333331, 65.000000000000000 },
 110   { 1.2006287819446431e+29, 0.33333333333333331, 70.000000000000000 },
 111   { 1.7213548977150022e+31, 0.33333333333333331, 75.000000000000000 },
 112   { 2.4734492458444449e+33, 0.33333333333333331, 80.000000000000000 },
 113   { 3.5611354547857122e+35, 0.33333333333333331, 85.000000000000000 },
 114   { 5.1360491295551829e+37, 0.33333333333333331, 90.000000000000000 },
 115   { 7.4189629097600431e+39, 0.33333333333333331, 95.000000000000000 },
 116   { 1.0731523308358370e+42, 0.33333333333333331, 100.00000000000000 },
 117 };
 118
 119 // Test function for nu=0.33333333333333331.
 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_cyl_bessel_i<double>);
 128   for (unsigned int i = 0; i < num_datum; ++i)
 129     {
 130       const Tp f = std::tr1::cyl_bessel_i(Tp(data002[i].nu), 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(1.0000000000000008e-12));
 144 }
 145
 146 // Test data for nu=0.50000000000000000.
 147 testcase_cyl_bessel_i<double> data003[] = {
 148   { 0.0000000000000000, 0.50000000000000000, 0.0000000000000000 },
 149   { 26.477547497559065, 0.50000000000000000, 5.0000000000000000 },
 150   { 2778.7846038745711, 0.50000000000000000, 10.000000000000000 },
 151   { 336729.88718706399, 0.50000000000000000, 15.000000000000000 },
 152   { 43279746.272428922, 0.50000000000000000, 20.000000000000000 },
 153   { 5745159748.3464680, 0.50000000000000000, 25.000000000000000 },
 154   { 778366068840.44580, 0.50000000000000000, 30.000000000000000 },
 155   { 106950522408567.66, 0.50000000000000000, 35.000000000000000 },
 156   { 14847705549021962., 0.50000000000000000, 40.000000000000000 },
 157   { 2.0775691824625661e+18, 0.50000000000000000, 45.000000000000000 },
 158   { 2.9251568529912988e+20, 0.50000000000000000, 50.000000000000000 },
 159   { 4.1392840094781220e+22, 0.50000000000000000, 55.000000000000000 },
 160   { 5.8817065760751945e+24, 0.50000000000000000, 60.000000000000000 },
 161   { 8.3867695787277231e+26, 0.50000000000000000, 65.000000000000000 },
 162   { 1.1994296461653203e+29, 0.50000000000000000, 70.000000000000000 },
 163   { 1.7197510246063332e+31, 0.50000000000000000, 75.000000000000000 },
 164   { 2.4712895036230794e+33, 0.50000000000000000, 80.000000000000000 },
 165   { 3.5582099086757769e+35, 0.50000000000000000, 85.000000000000000 },
 166   { 5.1320654031231090e+37, 0.50000000000000000, 90.000000000000000 },
 167   { 7.4135128383495227e+39, 0.50000000000000000, 95.000000000000000 },
 168   { 1.0724035825423179e+42, 0.50000000000000000, 100.00000000000000 },
 169 };
 170
 171 // Test function for nu=0.50000000000000000.
 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_cyl_bessel_i<double>);
 180   for (unsigned int i = 0; i < num_datum; ++i)
 181     {
 182       const Tp f = std::tr1::cyl_bessel_i(Tp(data003[i].nu), 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(1.0000000000000008e-12));
 196 }
 197
 198 // Test data for nu=0.66666666666666663.
 199 testcase_cyl_bessel_i<double> data004[] = {
 200   { 0.0000000000000000, 0.66666666666666663, 0.0000000000000000 },
 201   { 25.902310583215122, 0.66666666666666663, 5.0000000000000000 },
 202   { 2750.4090423459315, 0.66666666666666663, 10.000000000000000 },
 203   { 334476.98138574383, 0.66666666666666663, 15.000000000000000 },
 204   { 43064361.686912313, 0.66666666666666663, 20.000000000000000 },
 205   { 5722397441.9603882, 0.66666666666666663, 25.000000000000000 },
 206   { 775804343498.02661, 0.66666666666666663, 30.000000000000000 },
 207   { 106649495512800.89, 0.66666666666666663, 35.000000000000000 },
 208   { 14811199896983756., 0.66666666666666663, 40.000000000000000 },
 209   { 2.0730345814356961e+18, 0.66666666666666663, 45.000000000000000 },
 210   { 2.9194166755257467e+20, 0.66666666666666663, 50.000000000000000 },
 211   { 4.1319059569935366e+22, 0.66666666666666663, 55.000000000000000 },
 212   { 5.8721031476386222e+24, 0.66666666666666663, 60.000000000000000 },
 213   { 8.3741368248217830e+26, 0.66666666666666663, 65.000000000000000 },
 214   { 1.1977528777008688e+29, 0.66666666666666663, 70.000000000000000 },
 215   { 1.7175081240014333e+31, 0.66666666666666663, 75.000000000000000 },
 216   { 2.4682690458513916e+33, 0.66666666666666663, 80.000000000000000 },
 217   { 3.5541181975850724e+35, 0.66666666666666663, 85.000000000000000 },
 218   { 5.1264933963228864e+37, 0.66666666666666663, 90.000000000000000 },
 219   { 7.4058894880134064e+39, 0.66666666666666663, 95.000000000000000 },
 220   { 1.0713562154788124e+42, 0.66666666666666663, 100.00000000000000 },
 221 };
 222
 223 // Test function for nu=0.66666666666666663.
 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_cyl_bessel_i<double>);
 232   for (unsigned int i = 0; i < num_datum; ++i)
 233     {
 234       const Tp f = std::tr1::cyl_bessel_i(Tp(data004[i].nu), 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(5.0000000000000039e-13));
 248 }
 249
 250 // Test data for nu=1.0000000000000000.
 251 testcase_cyl_bessel_i<double> data005[] = {
 252   { 0.0000000000000000, 1.0000000000000000, 0.0000000000000000 },
 253   { 24.335642142450524, 1.0000000000000000, 5.0000000000000000 },
 254   { 2670.9883037012560, 1.0000000000000000, 10.000000000000000 },
 255   { 328124.92197020649, 1.0000000000000000, 15.000000000000000 },
 256   { 42454973.385127775, 1.0000000000000000, 20.000000000000000 },
 257   { 5657865129.8787022, 1.0000000000000000, 25.000000000000000 },
 258   { 768532038938.95667, 1.0000000000000000, 30.000000000000000 },
 259   { 105794126051896.17, 1.0000000000000000, 35.000000000000000 },
 260   { 14707396163259354., 1.0000000000000000, 40.000000000000000 },
 261   { 2.0601334620815775e+18, 1.0000000000000000, 45.000000000000000 },
 262   { 2.9030785901035635e+20, 1.0000000000000000, 50.000000000000000 },
 263   { 4.1108986452992812e+22, 1.0000000000000000, 55.000000000000000 },
 264   { 5.8447515883904527e+24, 1.0000000000000000, 60.000000000000000 },
 265   { 8.3381485471501302e+26, 1.0000000000000000, 65.000000000000000 },
 266   { 1.1929750788892366e+29, 1.0000000000000000, 70.000000000000000 },
 267   { 1.7111160152965384e+31, 1.0000000000000000, 75.000000000000000 },
 268   { 2.4596595795675343e+33, 1.0000000000000000, 80.000000000000000 },
 269   { 3.5424536064404024e+35, 1.0000000000000000, 85.000000000000000 },
 270   { 5.1106068152566129e+37, 1.0000000000000000, 90.000000000000000 },
 271   { 7.3841518091360157e+39, 1.0000000000000000, 95.000000000000000 },
 272   { 1.0683693903381569e+42, 1.0000000000000000, 100.00000000000000 },
 273 };
 274
 275 // Test function for nu=1.0000000000000000.
 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_cyl_bessel_i<double>);
 284   for (unsigned int i = 0; i < num_datum; ++i)
 285     {
 286       const Tp f = std::tr1::cyl_bessel_i(Tp(data005[i].nu), 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(5.0000000000000039e-13));
 300 }
 301
 302 // Test data for nu=2.0000000000000000.
 303 testcase_cyl_bessel_i<double> data006[] = {
 304   { 0.0000000000000000, 2.0000000000000000, 0.0000000000000000 },
 305   { 17.505614966624233, 2.0000000000000000, 5.0000000000000000 },
 306   { 2281.5189677260046, 2.0000000000000000, 10.000000000000000 },
 307   { 295899.38370188634, 2.0000000000000000, 15.000000000000000 },
 308   { 39312785.221040756, 2.0000000000000000, 20.000000000000000 },
 309   { 5321931396.0760136, 2.0000000000000000, 25.000000000000000 },
 310   { 730436828561.38013, 2.0000000000000000, 30.000000000000000 },
 311   { 101293439862977.19, 2.0000000000000000, 35.000000000000000 },
 312   { 14159404985256920., 2.0000000000000000, 40.000000000000000 },
 313   { 1.9918525879736883e+18, 2.0000000000000000, 45.000000000000000 },
 314   { 2.8164306402451938e+20, 2.0000000000000000, 50.000000000000000 },
 315   { 3.9993023372677515e+22, 2.0000000000000000, 55.000000000000000 },
 316   { 5.6992520026634433e+24, 2.0000000000000000, 60.000000000000000 },
 317   { 8.1464814287900378e+26, 2.0000000000000000, 65.000000000000000 },
 318   { 1.1675039556585663e+29, 2.0000000000000000, 70.000000000000000 },
 319   { 1.6770093176278926e+31, 2.0000000000000000, 75.000000000000000 },
 320   { 2.4136869148449879e+33, 2.0000000000000000, 80.000000000000000 },
 321   { 3.4801257808448186e+35, 2.0000000000000000, 85.000000000000000 },
 322   { 5.0256693051696307e+37, 2.0000000000000000, 90.000000000000000 },
 323   { 7.2678700343145842e+39, 2.0000000000000000, 95.000000000000000 },
 324   { 1.0523843193243042e+42, 2.0000000000000000, 100.00000000000000 },
 325 };
 326
 327 // Test function for nu=2.0000000000000000.
 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_cyl_bessel_i<double>);
 336   for (unsigned int i = 0; i < num_datum; ++i)
 337     {
 338       const Tp f = std::tr1::cyl_bessel_i(Tp(data006[i].nu), 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(5.0000000000000039e-13));
 352 }
 353
 354 // Test data for nu=5.0000000000000000.
 355 testcase_cyl_bessel_i<double> data007[] = {
 356   { 0.0000000000000000, 5.0000000000000000, 0.0000000000000000 },
 357   { 2.1579745473225476, 5.0000000000000000, 5.0000000000000000 },
 358   { 777.18828640326012, 5.0000000000000000, 10.000000000000000 },
 359   { 144572.01120063406, 5.0000000000000000, 15.000000000000000 },
 360   { 23018392.213413671, 5.0000000000000000, 20.000000000000000 },
 361   { 3472466208.7419176, 5.0000000000000000, 25.000000000000000 },
 362   { 512151465476.93494, 5.0000000000000000, 30.000000000000000 },
 363   { 74756743552251.531, 5.0000000000000000, 35.000000000000000 },
 364   { 10858318337624276., 5.0000000000000000, 40.000000000000000 },
 365   { 1.5736087399245906e+18, 5.0000000000000000, 45.000000000000000 },
 366   { 2.2785483079112829e+20, 5.0000000000000000, 50.000000000000000 },
 367   { 3.2989391052963687e+22, 5.0000000000000000, 55.000000000000000 },
 368   { 4.7777652072561742e+24, 5.0000000000000000, 60.000000000000000 },
 369   { 6.9232165147172671e+26, 5.0000000000000000, 65.000000000000000 },
 370   { 1.0038643002095153e+29, 5.0000000000000000, 70.000000000000000 },
 371   { 1.4566328222327068e+31, 5.0000000000000000, 75.000000000000000 },
 372   { 2.1151488565944838e+33, 5.0000000000000000, 80.000000000000000 },
 373   { 3.0735883450768236e+35, 5.0000000000000000, 85.000000000000000 },
 374   { 4.4694790189230327e+37, 5.0000000000000000, 90.000000000000000 },
 375   { 6.5037505570430971e+39, 5.0000000000000000, 95.000000000000000 },
 376   { 9.4700938730355882e+41, 5.0000000000000000, 100.00000000000000 },
 377 };
 378
 379 // Test function for nu=5.0000000000000000.
 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_cyl_bessel_i<double>);
 388   for (unsigned int i = 0; i < num_datum; ++i)
 389     {
 390       const Tp f = std::tr1::cyl_bessel_i(Tp(data007[i].nu), 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 nu=10.000000000000000.
 407 testcase_cyl_bessel_i<double> data008[] = {
 408   { 0.0000000000000000, 10.000000000000000, 0.0000000000000000 },
 409   { 0.0045800444191760525, 10.000000000000000, 5.0000000000000000 },
 410   { 21.891706163723381, 10.000000000000000, 10.000000000000000 },
 411   { 12267.475049806462, 10.000000000000000, 15.000000000000000 },
 412   { 3540200.2090195213, 10.000000000000000, 20.000000000000000 },
 413   { 771298871.17072666, 10.000000000000000, 25.000000000000000 },
 414   { 145831809975.96713, 10.000000000000000, 30.000000000000000 },
 415   { 25449470018534.785, 10.000000000000000, 35.000000000000000 },
 416   { 4228469210516757.5, 10.000000000000000, 40.000000000000000 },
 417   { 6.8049404557505152e+17, 10.000000000000000, 45.000000000000000 },
 418   { 1.0715971594776370e+20, 10.000000000000000, 50.000000000000000 },
 419   { 1.6618215752886714e+22, 10.000000000000000, 55.000000000000000 },
 420   { 2.5486246072566784e+24, 10.000000000000000, 60.000000000000000 },
 421   { 3.8764628702155481e+26, 10.000000000000000, 65.000000000000000 },
 422   { 5.8592538145409686e+28, 10.000000000000000, 70.000000000000000 },
 423   { 8.8135370711317444e+30, 10.000000000000000, 75.000000000000000 },
 424   { 1.3207418268325279e+33, 10.000000000000000, 80.000000000000000 },
 425   { 1.9732791360862186e+35, 10.000000000000000, 85.000000000000000 },
 426   { 2.9411893748384672e+37, 10.000000000000000, 90.000000000000000 },
 427   { 4.3754494922439990e+39, 10.000000000000000, 95.000000000000000 },
 428   { 6.4989755247201446e+41, 10.000000000000000, 100.00000000000000 },
 429 };
 430
 431 // Test function for nu=10.000000000000000.
 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_cyl_bessel_i<double>);
 440   for (unsigned int i = 0; i < num_datum; ++i)
 441     {
 442       const Tp f = std::tr1::cyl_bessel_i(Tp(data008[i].nu), 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 // Test data for nu=20.000000000000000.
 459 testcase_cyl_bessel_i<double> data009[] = {
 460   { 0.0000000000000000, 20.000000000000000, 0.0000000000000000 },
 461   { 5.0242393579718066e-11, 20.000000000000000, 5.0000000000000000 },
 462   { 0.00012507997356449481, 20.000000000000000, 10.000000000000000 },
 463   { 1.6470152535015834, 20.000000000000000, 15.000000000000000 },
 464   { 3188.7503288536154, 20.000000000000000, 20.000000000000000 },
 465   { 2449840.5422952301, 20.000000000000000, 25.000000000000000 },
 466   { 1126985104.4483771, 20.000000000000000, 30.000000000000000 },
 467   { 379617876611.88586, 20.000000000000000, 35.000000000000000 },
 468   { 104459633129479.89, 20.000000000000000, 40.000000000000000 },
 469   { 25039579987216504., 20.000000000000000, 45.000000000000000 },
 470   { 5.4420084027529964e+18, 20.000000000000000, 50.000000000000000 },
 471   { 1.1007498584335492e+21, 20.000000000000000, 55.000000000000000 },
 472   { 2.1091734863057236e+23, 20.000000000000000, 60.000000000000000 },
 473   { 3.8763618091286899e+25, 20.000000000000000, 65.000000000000000 },
 474   { 6.8946130527930859e+27, 20.000000000000000, 70.000000000000000 },
 475   { 1.1946319948836447e+30, 20.000000000000000, 75.000000000000000 },
 476   { 2.0265314377577580e+32, 20.000000000000000, 80.000000000000000 },
 477   { 3.3784665214179971e+34, 20.000000000000000, 85.000000000000000 },
 478   { 5.5516089411796670e+36, 20.000000000000000, 90.000000000000000 },
 479   { 9.0129310795305151e+38, 20.000000000000000, 95.000000000000000 },
 480   { 1.4483461256427176e+41, 20.000000000000000, 100.00000000000000 },
 481 };
 482
 483 // Test function for nu=20.000000000000000.
 484 template <typename Tp>
 485 void test009()
 486 {
 487   const Tp eps = std::numeric_limits<Tp>::epsilon();
 488   Tp max_abs_diff = -Tp(1);
 489   Tp max_abs_frac = -Tp(1);
 490   unsigned int num_datum = sizeof(data009)
 491                          / sizeof(testcase_cyl_bessel_i<double>);
 492   for (unsigned int i = 0; i < num_datum; ++i)
 493     {
 494       const Tp f = std::tr1::cyl_bessel_i(Tp(data009[i].nu), Tp(data009[i].x));
 495       const Tp f0 = data009[i].f0;
 496       const Tp diff = f - f0;
 497       if (std::abs(diff) > max_abs_diff)
 498         max_abs_diff = std::abs(diff);
 499       if (std::abs(f0) > Tp(10) * eps
 500        && std::abs(f) > Tp(10) * eps)
 501         {
 502           const Tp frac = diff / f0;
 503           if (std::abs(frac) > max_abs_frac)
 504             max_abs_frac = std::abs(frac);
 505         }
 506     }
 507   VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
 508 }
 509
 510 // Test data for nu=50.000000000000000.
 511 testcase_cyl_bessel_i<double> data010[] = {
 512   { 0.0000000000000000, 50.000000000000000, 0.0000000000000000 },
 513   { 2.9314696468108517e-45, 50.000000000000000, 5.0000000000000000 },
 514   { 4.7568945607268442e-30, 50.000000000000000, 10.000000000000000 },
 515   { 5.5468372730667069e-21, 50.000000000000000, 15.000000000000000 },
 516   { 2.2551205757604056e-14, 50.000000000000000, 20.000000000000000 },
 517   { 4.5344251866130282e-09, 50.000000000000000, 25.000000000000000 },
 518   { 0.00014590106916468937, 50.000000000000000, 30.000000000000000 },
 519   { 1.3965549457254878, 50.000000000000000, 35.000000000000000 },
 520   { 5726.8656631289878, 50.000000000000000, 40.000000000000000 },
 521   { 12672593.113027776, 50.000000000000000, 45.000000000000000 },
 522   { 17650802430.016705, 50.000000000000000, 50.000000000000000 },
 523   { 17220231607789.926, 50.000000000000000, 55.000000000000000 },
 524   { 12704607933652172., 50.000000000000000, 60.000000000000000 },
 525   { 7.4989491942193766e+18, 50.000000000000000, 65.000000000000000 },
 526   { 3.6944034898904901e+21, 50.000000000000000, 70.000000000000000 },
 527   { 1.5691634774370194e+24, 50.000000000000000, 75.000000000000000 },
 528   { 5.8927749458163587e+26, 50.000000000000000, 80.000000000000000 },
 529   { 1.9958849054749339e+29, 50.000000000000000, 85.000000000000000 },
 530   { 6.1946050361781518e+31, 50.000000000000000, 90.000000000000000 },
 531   { 1.7845429728697110e+34, 50.000000000000000, 95.000000000000000 },
 532   { 4.8219580855940813e+36, 50.000000000000000, 100.00000000000000 },
 533 };
 534
 535 // Test function for nu=50.000000000000000.
 536 template <typename Tp>
 537 void test010()
 538 {
 539   const Tp eps = std::numeric_limits<Tp>::epsilon();
 540   Tp max_abs_diff = -Tp(1);
 541   Tp max_abs_frac = -Tp(1);
 542   unsigned int num_datum = sizeof(data010)
 543                          / sizeof(testcase_cyl_bessel_i<double>);
 544   for (unsigned int i = 0; i < num_datum; ++i)
 545     {
 546       const Tp f = std::tr1::cyl_bessel_i(Tp(data010[i].nu), Tp(data010[i].x));
 547       const Tp f0 = data010[i].f0;
 548       const Tp diff = f - f0;
 549       if (std::abs(diff) > max_abs_diff)
 550         max_abs_diff = std::abs(diff);
 551       if (std::abs(f0) > Tp(10) * eps
 552        && std::abs(f) > Tp(10) * eps)
 553         {
 554           const Tp frac = diff / f0;
 555           if (std::abs(frac) > max_abs_frac)
 556             max_abs_frac = std::abs(frac);
 557         }
 558     }
 559   VERIFY(max_abs_frac < Tp(5.0000000000000029e-12));
 560 }
 561
 562 // Test data for nu=100.00000000000000.
 563 testcase_cyl_bessel_i<double> data011[] = {
 564   { 0.0000000000000000, 100.00000000000000, 0.0000000000000000 },
 565   { 7.0935514885313123e-119, 100.00000000000000, 5.0000000000000000 },
 566   { 1.0823442017492018e-88, 100.00000000000000, 10.000000000000000 },
 567   { 5.9887888536468904e-71, 100.00000000000000, 15.000000000000000 },
 568   { 2.8703193216428771e-58, 100.00000000000000, 20.000000000000000 },
 569   { 2.4426896913122370e-48, 100.00000000000000, 25.000000000000000 },
 570   { 3.9476420053334271e-40, 100.00000000000000, 30.000000000000000 },
 571   { 4.2836596180818801e-33, 100.00000000000000, 35.000000000000000 },
 572   { 6.6249380222596129e-27, 100.00000000000000, 40.000000000000000 },
 573   { 2.3702587262788881e-21, 100.00000000000000, 45.000000000000000 },
 574   { 2.7278879470966907e-16, 100.00000000000000, 50.000000000000000 },
 575   { 1.2763258878228088e-11, 100.00000000000000, 55.000000000000000 },
 576   { 2.8832770906491951e-07, 100.00000000000000, 60.000000000000000 },
 577   { 0.0035805902717061240, 100.00000000000000, 65.000000000000000 },
 578   { 27.017219102595398, 100.00000000000000, 70.000000000000000 },
 579   { 134001.44891209516, 100.00000000000000, 75.000000000000000 },
 580   { 465194832.85061038, 100.00000000000000, 80.000000000000000 },
 581   { 1189280653119.4819, 100.00000000000000, 85.000000000000000 },
 582   { 2334119331258731.0, 100.00000000000000, 90.000000000000000 },
 583   { 3.6399223078502380e+18, 100.00000000000000, 95.000000000000000 },
 584   { 4.6415349416161989e+21, 100.00000000000000, 100.00000000000000 },
 585 };
 586
 587 // Test function for nu=100.00000000000000.
 588 template <typename Tp>
 589 void test011()
 590 {
 591   const Tp eps = std::numeric_limits<Tp>::epsilon();
 592   Tp max_abs_diff = -Tp(1);
 593   Tp max_abs_frac = -Tp(1);
 594   unsigned int num_datum = sizeof(data011)
 595                          / sizeof(testcase_cyl_bessel_i<double>);
 596   for (unsigned int i = 0; i < num_datum; ++i)
 597     {
 598       const Tp f = std::tr1::cyl_bessel_i(Tp(data011[i].nu), Tp(data011[i].x));
 599       const Tp f0 = data011[i].f0;
 600       const Tp diff = f - f0;
 601       if (std::abs(diff) > max_abs_diff)
 602         max_abs_diff = std::abs(diff);
 603       if (std::abs(f0) > Tp(10) * eps
 604        && std::abs(f) > Tp(10) * eps)
 605         {
 606           const Tp frac = diff / f0;
 607           if (std::abs(frac) > max_abs_frac)
 608             max_abs_frac = std::abs(frac);
 609         }
 610     }
 611   VERIFY(max_abs_frac < Tp(2.5000000000000014e-11));
 612 }
 613
 614 int main(int, char**)
 615 {
 616   test001<double>();
 617   test002<double>();
 618   test003<double>();
 619   test004<double>();
 620   test005<double>();
 621   test006<double>();
 622   test007<double>();
 623   test008<double>();
 624   test009<double>();
 625   test010<double>();
 626   test011<double>();
 627   return 0;
 628 }