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