git.gag.com Git - fw/altos/blob - ao-tools/lib/i0.c

   1 /*
   2  * This file comes from the cephes math library, which was
   3  * released under the GPLV2+ license as a part of the Debian labplot
   4  * package (I've included the GPLV2 license reference here to make
   5  * this clear) - Keith Packard <keithp@keithp.com>
   6  *
   7  * Cephes Math Library Release 2.0:  April, 1987
   8  * Copyright 1984, 1987 by Stephen L. Moshier
   9  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  10  *
  11  * This program is free software; you can redistribute it and/or modify
  12  * it under the terms of the GNU General Public License as published by
  13  * the Free Software Foundation; either version 2 of the License, or
  14  * (at your option) any later version.
  15  *
  16  * This program is distributed in the hope that it will be useful, but
  17  * WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU General Public License along
  22  * with this program; if not, write to the Free Software Foundation, Inc.,
  23  * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
  24  */
  25
  26 /*                                                      i0.c
  27  *
  28  *      Modified Bessel function of order zero
  29  *
  30  *
  31  *
  32  * SYNOPSIS:
  33  *
  34  * double x, y, i0();
  35  *
  36  * y = i0( x );
  37  *
  38  *
  39  *
  40  * DESCRIPTION:
  41  *
  42  * Returns modified Bessel function of order zero of the
  43  * argument.
  44  *
  45  * The function is defined as i0(x) = j0( ix ).
  46  *
  47  * The range is partitioned into the two intervals [0,8] and
  48  * (8, infinity).  Chebyshev polynomial expansions are employed
  49  * in each interval.
  50  *
  51  *
  52  *
  53  * ACCURACY:
  54  *
  55  *                      Relative error:
  56  * arithmetic   domain     # trials      peak         rms
  57  *    DEC       0,30         6000       8.2e-17     1.9e-17
  58  *    IEEE      0,30        30000       5.8e-16     1.4e-16
  59  *
  60  */
  61 \f/*                                                     i0e.c
  62  *
  63  *      Modified Bessel function of order zero,
  64  *      exponentially scaled
  65  *
  66  *
  67  *
  68  * SYNOPSIS:
  69  *
  70  * double x, y, i0e();
  71  *
  72  * y = i0e( x );
  73  *
  74  *
  75  *
  76  * DESCRIPTION:
  77  *
  78  * Returns exponentially scaled modified Bessel function
  79  * of order zero of the argument.
  80  *
  81  * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
  82  *
  83  *
  84  *
  85  * ACCURACY:
  86  *
  87  *                      Relative error:
  88  * arithmetic   domain     # trials      peak         rms
  89  *    IEEE      0,30        30000       5.4e-16     1.2e-16
  90  * See i0().
  91  *
  92  */
  93 \f
  94 /*                                                      i0.c            */
  95
  96
  97 /*
  98 Cephes Math Library Release 2.0:  April, 1987
  99 Copyright 1984, 1987 by Stephen L. Moshier
 100 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 101 */
 102
 103 #include <math.h>
 104 #include "mconf.h"
 105 #include "cephes.h"
 106
 107 /* Chebyshev coefficients for exp(-x) I0(x)
 108  * in the interval [0,8].
 109  *
 110  * lim(x->0){ exp(-x) I0(x) } = 1.
 111  */
 112
 113 #ifdef UNK
 114 static double A[] =
 115 {
 116 -4.41534164647933937950E-18,
 117  3.33079451882223809783E-17,
 118 -2.43127984654795469359E-16,
 119  1.71539128555513303061E-15,
 120 -1.16853328779934516808E-14,
 121  7.67618549860493561688E-14,
 122 -4.85644678311192946090E-13,
 123  2.95505266312963983461E-12,
 124 -1.72682629144155570723E-11,
 125  9.67580903537323691224E-11,
 126 -5.18979560163526290666E-10,
 127  2.65982372468238665035E-9,
 128 -1.30002500998624804212E-8,
 129  6.04699502254191894932E-8,
 130 -2.67079385394061173391E-7,
 131  1.11738753912010371815E-6,
 132 -4.41673835845875056359E-6,
 133  1.64484480707288970893E-5,
 134 -5.75419501008210370398E-5,
 135  1.88502885095841655729E-4,
 136 -5.76375574538582365885E-4,
 137  1.63947561694133579842E-3,
 138 -4.32430999505057594430E-3,
 139  1.05464603945949983183E-2,
 140 -2.37374148058994688156E-2,
 141  4.93052842396707084878E-2,
 142 -9.49010970480476444210E-2,
 143  1.71620901522208775349E-1,
 144 -3.04682672343198398683E-1,
 145  6.76795274409476084995E-1
 146 };
 147 #endif
 148
 149 #ifdef DEC
 150 static unsigned short A[] = {
 151 0121642,0162671,0004646,0103567,
 152 0022431,0115424,0135755,0026104,
 153 0123214,0023533,0110365,0156635,
 154 0023767,0033304,0117662,0172716,
 155 0124522,0100426,0012277,0157531,
 156 0025254,0155062,0054461,0030465,
 157 0126010,0131143,0013560,0153604,
 158 0026517,0170577,0006336,0114437,
 159 0127227,0162253,0152243,0052734,
 160 0027724,0142766,0061641,0160200,
 161 0130416,0123760,0116564,0125262,
 162 0031066,0144035,0021246,0054641,
 163 0131537,0053664,0060131,0102530,
 164 0032201,0155664,0165153,0020652,
 165 0132617,0061434,0074423,0176145,
 166 0033225,0174444,0136147,0122542,
 167 0133624,0031576,0056453,0020470,
 168 0034211,0175305,0172321,0041314,
 169 0134561,0054462,0147040,0165315,
 170 0035105,0124333,0120203,0162532,
 171 0135427,0013750,0174257,0055221,
 172 0035726,0161654,0050220,0100162,
 173 0136215,0131361,0000325,0041110,
 174 0036454,0145417,0117357,0017352,
 175 0136702,0072367,0104415,0133574,
 176 0037111,0172126,0072505,0014544,
 177 0137302,0055601,0120550,0033523,
 178 0037457,0136543,0136544,0043002,
 179 0137633,0177536,0001276,0066150,
 180 0040055,0041164,0100655,0010521
 181 };
 182 #endif
 183
 184 #ifdef IBMPC
 185 static unsigned short A[] = {
 186 0xd0ef,0x2134,0x5cb7,0xbc54,
 187 0xa589,0x977d,0x3362,0x3c83,
 188 0xbbb4,0x721e,0x84eb,0xbcb1,
 189 0x5eba,0x93f6,0xe6d8,0x3cde,
 190 0xfbeb,0xc297,0x5022,0xbd0a,
 191 0x2627,0x4b26,0x9b46,0x3d35,
 192 0x1af0,0x62ee,0x164c,0xbd61,
 193 0xd324,0xe19b,0xfe2f,0x3d89,
 194 0x6abc,0x7a94,0xfc95,0xbdb2,
 195 0x3c10,0xcc74,0x98be,0x3dda,
 196 0x9556,0x13ae,0xd4fe,0xbe01,
 197 0xcb34,0xa454,0xd903,0x3e26,
 198 0x30ab,0x8c0b,0xeaf6,0xbe4b,
 199 0x6435,0x9d4d,0x3b76,0x3e70,
 200 0x7f8d,0x8f22,0xec63,0xbe91,
 201 0xf4ac,0x978c,0xbf24,0x3eb2,
 202 0x6427,0xcba5,0x866f,0xbed2,
 203 0x2859,0xbe9a,0x3f58,0x3ef1,
 204 0x1d5a,0x59c4,0x2b26,0xbf0e,
 205 0x7cab,0x7410,0xb51b,0x3f28,
 206 0xeb52,0x1f15,0xe2fd,0xbf42,
 207 0x100e,0x8a12,0xdc75,0x3f5a,
 208 0xa849,0x201a,0xb65e,0xbf71,
 209 0xe3dd,0xf3dd,0x9961,0x3f85,
 210 0xb6f0,0xf121,0x4e9e,0xbf98,
 211 0xa32d,0xcea8,0x3e8a,0x3fa9,
 212 0x06ea,0x342d,0x4b70,0xbfb8,
 213 0x88c0,0x77ac,0xf7ac,0x3fc5,
 214 0xcd8d,0xc057,0x7feb,0xbfd3,
 215 0xa22a,0x9035,0xa84e,0x3fe5,
 216 };
 217 #endif
 218
 219 #ifdef MIEEE
 220 static unsigned short A[] = {
 221 0xbc54,0x5cb7,0x2134,0xd0ef,
 222 0x3c83,0x3362,0x977d,0xa589,
 223 0xbcb1,0x84eb,0x721e,0xbbb4,
 224 0x3cde,0xe6d8,0x93f6,0x5eba,
 225 0xbd0a,0x5022,0xc297,0xfbeb,
 226 0x3d35,0x9b46,0x4b26,0x2627,
 227 0xbd61,0x164c,0x62ee,0x1af0,
 228 0x3d89,0xfe2f,0xe19b,0xd324,
 229 0xbdb2,0xfc95,0x7a94,0x6abc,
 230 0x3dda,0x98be,0xcc74,0x3c10,
 231 0xbe01,0xd4fe,0x13ae,0x9556,
 232 0x3e26,0xd903,0xa454,0xcb34,
 233 0xbe4b,0xeaf6,0x8c0b,0x30ab,
 234 0x3e70,0x3b76,0x9d4d,0x6435,
 235 0xbe91,0xec63,0x8f22,0x7f8d,
 236 0x3eb2,0xbf24,0x978c,0xf4ac,
 237 0xbed2,0x866f,0xcba5,0x6427,
 238 0x3ef1,0x3f58,0xbe9a,0x2859,
 239 0xbf0e,0x2b26,0x59c4,0x1d5a,
 240 0x3f28,0xb51b,0x7410,0x7cab,
 241 0xbf42,0xe2fd,0x1f15,0xeb52,
 242 0x3f5a,0xdc75,0x8a12,0x100e,
 243 0xbf71,0xb65e,0x201a,0xa849,
 244 0x3f85,0x9961,0xf3dd,0xe3dd,
 245 0xbf98,0x4e9e,0xf121,0xb6f0,
 246 0x3fa9,0x3e8a,0xcea8,0xa32d,
 247 0xbfb8,0x4b70,0x342d,0x06ea,
 248 0x3fc5,0xf7ac,0x77ac,0x88c0,
 249 0xbfd3,0x7feb,0xc057,0xcd8d,
 250 0x3fe5,0xa84e,0x9035,0xa22a
 251 };
 252 #endif
 253
 254
 255 /* Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
 256  * in the inverted interval [8,infinity].
 257  *
 258  * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
 259  */
 260
 261 #ifdef UNK
 262 static double B[] =
 263 {
 264 -7.23318048787475395456E-18,
 265 -4.83050448594418207126E-18,
 266  4.46562142029675999901E-17,
 267  3.46122286769746109310E-17,
 268 -2.82762398051658348494E-16,
 269 -3.42548561967721913462E-16,
 270  1.77256013305652638360E-15,
 271  3.81168066935262242075E-15,
 272 -9.55484669882830764870E-15,
 273 -4.15056934728722208663E-14,
 274  1.54008621752140982691E-14,
 275  3.85277838274214270114E-13,
 276  7.18012445138366623367E-13,
 277 -1.79417853150680611778E-12,
 278 -1.32158118404477131188E-11,
 279 -3.14991652796324136454E-11,
 280  1.18891471078464383424E-11,
 281  4.94060238822496958910E-10,
 282  3.39623202570838634515E-9,
 283  2.26666899049817806459E-8,
 284  2.04891858946906374183E-7,
 285  2.89137052083475648297E-6,
 286  6.88975834691682398426E-5,
 287  3.36911647825569408990E-3,
 288  8.04490411014108831608E-1
 289 };
 290 #endif
 291
 292 #ifdef DEC
 293 static unsigned short B[] = {
 294 0122005,0066672,0123124,0054311,
 295 0121662,0033323,0030214,0104602,
 296 0022515,0170300,0113314,0020413,
 297 0022437,0117350,0035402,0007146,
 298 0123243,0000135,0057220,0177435,
 299 0123305,0073476,0144106,0170702,
 300 0023777,0071755,0017527,0154373,
 301 0024211,0052214,0102247,0033270,
 302 0124454,0017763,0171453,0012322,
 303 0125072,0166316,0075505,0154616,
 304 0024612,0133770,0065376,0025045,
 305 0025730,0162143,0056036,0001632,
 306 0026112,0015077,0150464,0063542,
 307 0126374,0101030,0014274,0065457,
 308 0127150,0077271,0125763,0157617,
 309 0127412,0104350,0040713,0120445,
 310 0027121,0023765,0057500,0001165,
 311 0030407,0147146,0003643,0075644,
 312 0031151,0061445,0044422,0156065,
 313 0031702,0132224,0003266,0125551,
 314 0032534,0000076,0147153,0005555,
 315 0033502,0004536,0004016,0026055,
 316 0034620,0076433,0142314,0171215,
 317 0036134,0146145,0013454,0101104,
 318 0040115,0171425,0062500,0047133
 319 };
 320 #endif
 321
 322 #ifdef IBMPC
 323 static unsigned short B[] = {
 324 0x8b19,0x54ca,0xadb7,0xbc60,
 325 0x9130,0x6611,0x46da,0xbc56,
 326 0x8421,0x12d9,0xbe18,0x3c89,
 327 0x41cd,0x0760,0xf3dd,0x3c83,
 328 0x1fe4,0xabd2,0x600b,0xbcb4,
 329 0xde38,0xd908,0xaee7,0xbcb8,
 330 0xfb1f,0xa3ea,0xee7d,0x3cdf,
 331 0xe6d7,0x9094,0x2a91,0x3cf1,
 332 0x629a,0x7e65,0x83fe,0xbd05,
 333 0xbb32,0xcf68,0x5d99,0xbd27,
 334 0xc545,0x0d5f,0x56ff,0x3d11,
 335 0xc073,0x6b83,0x1c8c,0x3d5b,
 336 0x8cec,0xfa26,0x4347,0x3d69,
 337 0x8d66,0x0317,0x9043,0xbd7f,
 338 0x7bf2,0x357e,0x0fd7,0xbdad,
 339 0x7425,0x0839,0x511d,0xbdc1,
 340 0x004f,0xabe8,0x24fe,0x3daa,
 341 0x6f75,0xc0f4,0xf9cc,0x3e00,
 342 0x5b87,0xa922,0x2c64,0x3e2d,
 343 0xd56d,0x80d6,0x5692,0x3e58,
 344 0x616e,0xd9cd,0x8007,0x3e8b,
 345 0xc586,0xc101,0x412b,0x3ec8,
 346 0x9e52,0x7899,0x0fa3,0x3f12,
 347 0x9049,0xa2e5,0x998c,0x3f6b,
 348 0x09cb,0xaca8,0xbe62,0x3fe9
 349 };
 350 #endif
 351
 352 #ifdef MIEEE
 353 static unsigned short B[] = {
 354 0xbc60,0xadb7,0x54ca,0x8b19,
 355 0xbc56,0x46da,0x6611,0x9130,
 356 0x3c89,0xbe18,0x12d9,0x8421,
 357 0x3c83,0xf3dd,0x0760,0x41cd,
 358 0xbcb4,0x600b,0xabd2,0x1fe4,
 359 0xbcb8,0xaee7,0xd908,0xde38,
 360 0x3cdf,0xee7d,0xa3ea,0xfb1f,
 361 0x3cf1,0x2a91,0x9094,0xe6d7,
 362 0xbd05,0x83fe,0x7e65,0x629a,
 363 0xbd27,0x5d99,0xcf68,0xbb32,
 364 0x3d11,0x56ff,0x0d5f,0xc545,
 365 0x3d5b,0x1c8c,0x6b83,0xc073,
 366 0x3d69,0x4347,0xfa26,0x8cec,
 367 0xbd7f,0x9043,0x0317,0x8d66,
 368 0xbdad,0x0fd7,0x357e,0x7bf2,
 369 0xbdc1,0x511d,0x0839,0x7425,
 370 0x3daa,0x24fe,0xabe8,0x004f,
 371 0x3e00,0xf9cc,0xc0f4,0x6f75,
 372 0x3e2d,0x2c64,0xa922,0x5b87,
 373 0x3e58,0x5692,0x80d6,0xd56d,
 374 0x3e8b,0x8007,0xd9cd,0x616e,
 375 0x3ec8,0x412b,0xc101,0xc586,
 376 0x3f12,0x0fa3,0x7899,0x9e52,
 377 0x3f6b,0x998c,0xa2e5,0x9049,
 378 0x3fe9,0xbe62,0xaca8,0x09cb
 379 };
 380 #endif
 381
 382 double i0(double x)
 383 {
 384 double y;
 385
 386 if( x < 0 )
 387         x = -x;
 388 if( x <= 8.0 )
 389         {
 390         y = (x/2.0) - 2.0;
 391         return( exp(x) * chbevl( y, A, 30 ) );
 392         }
 393
 394 return(  exp(x) * chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) );
 395
 396 }
 397
 398
 399
 400
 401 double i0e(double x )
 402 {
 403 double y;
 404
 405 if( x < 0 )
 406         x = -x;
 407 if( x <= 8.0 )
 408         {
 409         y = (x/2.0) - 2.0;
 410         return( chbevl( y, A, 30 ) );
 411         }
 412
 413 return(  chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) );
 414
 415 }