i1.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271
  1. /* i1.c
  2. *
  3. * Modified Bessel function of order one
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, i1();
  10. *
  11. * y = i1( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns modified Bessel function of order one of the
  18. * argument.
  19. *
  20. * The function is defined as i1(x) = -i j1( ix ).
  21. *
  22. * The range is partitioned into the two intervals [0,8] and
  23. * (8, infinity). Chebyshev polynomial expansions are employed
  24. * in each interval.
  25. *
  26. *
  27. *
  28. * ACCURACY:
  29. *
  30. * Relative error:
  31. * arithmetic domain # trials peak rms
  32. * DEC 0, 30 3400 1.2e-16 2.3e-17
  33. * IEEE 0, 30 30000 1.9e-15 2.1e-16
  34. *
  35. *
  36. */
  37. /* i1e.c
  38. *
  39. * Modified Bessel function of order one,
  40. * exponentially scaled
  41. *
  42. *
  43. *
  44. * SYNOPSIS:
  45. *
  46. * double x, y, i1e();
  47. *
  48. * y = i1e( x );
  49. *
  50. *
  51. *
  52. * DESCRIPTION:
  53. *
  54. * Returns exponentially scaled modified Bessel function
  55. * of order one of the argument.
  56. *
  57. * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
  58. *
  59. *
  60. *
  61. * ACCURACY:
  62. *
  63. * Relative error:
  64. * arithmetic domain # trials peak rms
  65. * IEEE 0, 30 30000 2.0e-15 2.0e-16
  66. * See i1().
  67. *
  68. */
  69. /* i1.c 2 */
  70. /*
  71. Cephes Math Library Release 2.8: June, 2000
  72. Copyright 1985, 1987, 2000 by Stephen L. Moshier
  73. */
  74. #include "mconf.h"
  75. /* Chebyshev coefficients for exp(-x) I1(x) / x
  76. * in the interval [0,8].
  77. *
  78. * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
  79. */
  80. #ifdef UNK
  81. static double A[] = {2.77791411276104639959E-18, -2.11142121435816608115E-17,
  82. 1.55363195773620046921E-16, -1.10559694773538630805E-15,
  83. 7.60068429473540693410E-15, -5.04218550472791168711E-14,
  84. 3.22379336594557470981E-13, -1.98397439776494371520E-12,
  85. 1.17361862988909016308E-11, -6.66348972350202774223E-11,
  86. 3.62559028155211703701E-10, -1.88724975172282928790E-9,
  87. 9.38153738649577178388E-9, -4.44505912879632808065E-8,
  88. 2.00329475355213526229E-7, -8.56872026469545474066E-7,
  89. 3.47025130813767847674E-6, -1.32731636560394358279E-5,
  90. 4.78156510755005422638E-5, -1.61760815825896745588E-4,
  91. 5.12285956168575772895E-4, -1.51357245063125314899E-3,
  92. 4.15642294431288815669E-3, -1.05640848946261981558E-2,
  93. 2.47264490306265168283E-2, -5.29459812080949914269E-2,
  94. 1.02643658689847095384E-1, -1.76416518357834055153E-1,
  95. 2.52587186443633654823E-1};
  96. #endif
  97. #ifdef DEC
  98. static unsigned short A[] = {
  99. 0021514, 0174520, 0060742, 0000241, 0122302, 0137206, 0016120, 0025663,
  100. 0023063, 0017437, 0026235, 0176536, 0123637, 0052523, 0170150, 0125632,
  101. 0024410, 0165770, 0030251, 0044134, 0125143, 0012160, 0162170, 0054727,
  102. 0025665, 0075702, 0035716, 0145247, 0126413, 0116032, 0176670, 0015462,
  103. 0027116, 0073425, 0110351, 0105242, 0127622, 0104034, 0137530, 0037364,
  104. 0030307, 0050645, 0120776, 0175535, 0131001, 0130331, 0043523, 0037455,
  105. 0031441, 0026160, 0010712, 0100174, 0132076, 0164761, 0022706, 0017500,
  106. 0032527, 0015045, 0115076, 0104076, 0133146, 0001714, 0015434, 0144520,
  107. 0033550, 0161166, 0124215, 0077050, 0134136, 0127715, 0143365, 0157170,
  108. 0034510, 0106652, 0013070, 0064130, 0135051, 0117126, 0117264, 0123761,
  109. 0035406, 0045355, 0133066, 0175751, 0135706, 0061420, 0054746, 0122440,
  110. 0036210, 0031232, 0047235, 0006640, 0136455, 0012373, 0144235, 0011523,
  111. 0036712, 0107437, 0036731, 0015111, 0137130, 0156742, 0115744, 0172743,
  112. 0037322, 0033326, 0124667, 0124740, 0137464, 0123210, 0021510, 0144556,
  113. 0037601, 0051433, 0111123, 0177721};
  114. #endif
  115. #ifdef IBMPC
  116. static unsigned short A[] = {
  117. 0x4014, 0x0c3c, 0x9f2a, 0x3c49, 0x0576, 0xc38a, 0x57d0, 0xbc78, 0xbfac,
  118. 0xe593, 0x63e3, 0x3ca6, 0x1573, 0x7e0d, 0xeaaa, 0xbcd3, 0x290c, 0x0615,
  119. 0x1d7f, 0x3d01, 0x0b3b, 0x1c8f, 0x628e, 0xbd2c, 0xd955, 0x4779, 0xaf78,
  120. 0x3d56, 0x0366, 0x5fb7, 0x7383, 0xbd81, 0x3154, 0xb21d, 0xcee2, 0x3da9,
  121. 0x07de, 0x97eb, 0x5103, 0xbdd2, 0xdf6c, 0xb43f, 0xea34, 0x3df8, 0x67e6,
  122. 0x28ea, 0x361b, 0xbe20, 0x5010, 0x0239, 0x258e, 0x3e44, 0xc3e8, 0x24b8,
  123. 0xdd3e, 0xbe67, 0xd108, 0xb347, 0xe344, 0x3e8a, 0x992a, 0x8363, 0xc079,
  124. 0xbeac, 0xafc5, 0xd511, 0x1c4e, 0x3ecd, 0xbbcf, 0xb8de, 0xd5f9, 0xbeeb,
  125. 0x0d0b, 0x42c7, 0x11b5, 0x3f09, 0x94fe, 0xd3d6, 0x33ca, 0xbf25, 0xdf7d,
  126. 0xb6c6, 0xc95d, 0x3f40, 0xd4a4, 0x0b3c, 0xcc62, 0xbf58, 0xa1b4, 0x49d3,
  127. 0x0653, 0x3f71, 0xa26a, 0x7913, 0xa29f, 0xbf85, 0x2349, 0xe7bb, 0x51e3,
  128. 0x3f99, 0x9ebc, 0x537c, 0x1bbc, 0xbfab, 0xf53c, 0xd536, 0x46da, 0x3fba,
  129. 0x192e, 0x0469, 0x94d1, 0xbfc6, 0x7ffa, 0x724a, 0x2a63, 0x3fd0};
  130. #endif
  131. #ifdef MIEEE
  132. static unsigned short A[] = {
  133. 0x3c49, 0x9f2a, 0x0c3c, 0x4014, 0xbc78, 0x57d0, 0xc38a, 0x0576, 0x3ca6,
  134. 0x63e3, 0xe593, 0xbfac, 0xbcd3, 0xeaaa, 0x7e0d, 0x1573, 0x3d01, 0x1d7f,
  135. 0x0615, 0x290c, 0xbd2c, 0x628e, 0x1c8f, 0x0b3b, 0x3d56, 0xaf78, 0x4779,
  136. 0xd955, 0xbd81, 0x7383, 0x5fb7, 0x0366, 0x3da9, 0xcee2, 0xb21d, 0x3154,
  137. 0xbdd2, 0x5103, 0x97eb, 0x07de, 0x3df8, 0xea34, 0xb43f, 0xdf6c, 0xbe20,
  138. 0x361b, 0x28ea, 0x67e6, 0x3e44, 0x258e, 0x0239, 0x5010, 0xbe67, 0xdd3e,
  139. 0x24b8, 0xc3e8, 0x3e8a, 0xe344, 0xb347, 0xd108, 0xbeac, 0xc079, 0x8363,
  140. 0x992a, 0x3ecd, 0x1c4e, 0xd511, 0xafc5, 0xbeeb, 0xd5f9, 0xb8de, 0xbbcf,
  141. 0x3f09, 0x11b5, 0x42c7, 0x0d0b, 0xbf25, 0x33ca, 0xd3d6, 0x94fe, 0x3f40,
  142. 0xc95d, 0xb6c6, 0xdf7d, 0xbf58, 0xcc62, 0x0b3c, 0xd4a4, 0x3f71, 0x0653,
  143. 0x49d3, 0xa1b4, 0xbf85, 0xa29f, 0x7913, 0xa26a, 0x3f99, 0x51e3, 0xe7bb,
  144. 0x2349, 0xbfab, 0x1bbc, 0x537c, 0x9ebc, 0x3fba, 0x46da, 0xd536, 0xf53c,
  145. 0xbfc6, 0x94d1, 0x0469, 0x192e, 0x3fd0, 0x2a63, 0x724a, 0x7ffa};
  146. #endif
  147. /* i1.c */
  148. /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
  149. * in the inverted interval [8,infinity].
  150. *
  151. * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
  152. */
  153. #ifdef UNK
  154. static double B[] = {7.51729631084210481353E-18, 4.41434832307170791151E-18,
  155. -4.65030536848935832153E-17, -3.20952592199342395980E-17,
  156. 2.96262899764595013876E-16, 3.30820231092092828324E-16,
  157. -1.88035477551078244854E-15, -3.81440307243700780478E-15,
  158. 1.04202769841288027642E-14, 4.27244001671195135429E-14,
  159. -2.10154184277266431302E-14, -4.08355111109219731823E-13,
  160. -7.19855177624590851209E-13, 2.03562854414708950722E-12,
  161. 1.41258074366137813316E-11, 3.25260358301548823856E-11,
  162. -1.89749581235054123450E-11, -5.58974346219658380687E-10,
  163. -3.83538038596423702205E-9, -2.63146884688951950684E-8,
  164. -2.51223623787020892529E-7, -3.88256480887769039346E-6,
  165. -1.10588938762623716291E-4, -9.76109749136146840777E-3,
  166. 7.78576235018280120474E-1};
  167. #endif
  168. #ifdef DEC
  169. static unsigned short B[] = {
  170. 0022012, 0125555, 0115227, 0043456, 0021642, 0156127, 0052075, 0145203,
  171. 0122526, 0072435, 0111231, 0011664, 0122424, 0001544, 0161671, 0114403,
  172. 0023252, 0144257, 0163532, 0142121, 0023276, 0132162, 0174045, 0013204,
  173. 0124007, 0077154, 0057046, 0110517, 0124211, 0066650, 0116127, 0157073,
  174. 0024473, 0133413, 0130551, 0107504, 0025100, 0064741, 0032631, 0040364,
  175. 0124675, 0045101, 0071551, 0012400, 0125745, 0161054, 0071637, 0011247,
  176. 0126112, 0117410, 0035525, 0122231, 0026417, 0037237, 0131034, 0176427,
  177. 0027170, 0100373, 0024742, 0025725, 0027417, 0006417, 0105303, 0141446,
  178. 0127246, 0163716, 0121202, 0060137, 0130431, 0123122, 0120436, 0166000,
  179. 0131203, 0144134, 0153251, 0124500, 0131742, 0005234, 0122732, 0033006,
  180. 0132606, 0157751, 0072362, 0121031, 0133602, 0043372, 0047120, 0015626,
  181. 0134747, 0165774, 0001125, 0046462, 0136437, 0166402, 0117746, 0155137,
  182. 0040107, 0050305, 0125330, 0124241};
  183. #endif
  184. #ifdef IBMPC
  185. static unsigned short B[] = {
  186. 0xe8e6, 0xb352, 0x556d, 0x3c61, 0xb950, 0xea87, 0x5b8a, 0x3c54, 0x2277,
  187. 0xb253, 0xcea3, 0xbc8a, 0x3320, 0x9c77, 0x806c, 0xbc82, 0x588a, 0xfceb,
  188. 0x5915, 0x3cb5, 0xa2d1, 0x5f04, 0xd68e, 0x3cb7, 0xd22a, 0x8bc4, 0xefcd,
  189. 0xbce0, 0xfbc7, 0x138a, 0x2db5, 0xbcf1, 0x31e8, 0x762d, 0x76e1, 0x3d07,
  190. 0x281e, 0x26b3, 0x0d3c, 0x3d28, 0x22a0, 0x2e6d, 0xa948, 0xbd17, 0xe255,
  191. 0x8e73, 0xbc45, 0xbd5c, 0xb493, 0x076a, 0x53e1, 0xbd69, 0x9fa3, 0xf643,
  192. 0xe7d3, 0x3d81, 0x457b, 0x653c, 0x101f, 0x3daf, 0x7865, 0xf158, 0xe1a1,
  193. 0x3dc1, 0x4c0c, 0xd450, 0xdcf9, 0xbdb4, 0xdd80, 0x5423, 0x34ca, 0xbe03,
  194. 0x3528, 0x9ad5, 0x790b, 0xbe30, 0x46c1, 0x94bb, 0x4153, 0xbe5c, 0x5443,
  195. 0x2e9e, 0xdbfd, 0xbe90, 0x0373, 0x49ca, 0x48df, 0xbed0, 0xa9a6, 0x804a,
  196. 0xfd7f, 0xbf1c, 0xdb4c, 0x53fc, 0xfda0, 0xbf83, 0x1514, 0xb55b, 0xea18,
  197. 0x3fe8};
  198. #endif
  199. #ifdef MIEEE
  200. static unsigned short B[] = {
  201. 0x3c61, 0x556d, 0xb352, 0xe8e6, 0x3c54, 0x5b8a, 0xea87, 0xb950, 0xbc8a,
  202. 0xcea3, 0xb253, 0x2277, 0xbc82, 0x806c, 0x9c77, 0x3320, 0x3cb5, 0x5915,
  203. 0xfceb, 0x588a, 0x3cb7, 0xd68e, 0x5f04, 0xa2d1, 0xbce0, 0xefcd, 0x8bc4,
  204. 0xd22a, 0xbcf1, 0x2db5, 0x138a, 0xfbc7, 0x3d07, 0x76e1, 0x762d, 0x31e8,
  205. 0x3d28, 0x0d3c, 0x26b3, 0x281e, 0xbd17, 0xa948, 0x2e6d, 0x22a0, 0xbd5c,
  206. 0xbc45, 0x8e73, 0xe255, 0xbd69, 0x53e1, 0x076a, 0xb493, 0x3d81, 0xe7d3,
  207. 0xf643, 0x9fa3, 0x3daf, 0x101f, 0x653c, 0x457b, 0x3dc1, 0xe1a1, 0xf158,
  208. 0x7865, 0xbdb4, 0xdcf9, 0xd450, 0x4c0c, 0xbe03, 0x34ca, 0x5423, 0xdd80,
  209. 0xbe30, 0x790b, 0x9ad5, 0x3528, 0xbe5c, 0x4153, 0x94bb, 0x46c1, 0xbe90,
  210. 0xdbfd, 0x2e9e, 0x5443, 0xbed0, 0x48df, 0x49ca, 0x0373, 0xbf1c, 0xfd7f,
  211. 0x804a, 0xa9a6, 0xbf83, 0xfda0, 0x53fc, 0xdb4c, 0x3fe8, 0xea18, 0xb55b,
  212. 0x1514};
  213. #endif
  214. /* i1.c */
  215. #ifdef ANSIPROT
  216. extern double chbevl(double, void *, int);
  217. extern double exp(double);
  218. extern double sqrt(double);
  219. extern double fabs(double);
  220. #else
  221. double chbevl(), exp(), sqrt(), fabs();
  222. #endif
  223. double i1(x) double x;
  224. {
  225. double y, z;
  226. z = fabs(x);
  227. if (z <= 8.0) {
  228. y = (z / 2.0) - 2.0;
  229. z = chbevl(y, A, 29) * z * exp(z);
  230. } else {
  231. z = exp(z) * chbevl(32.0 / z - 2.0, B, 25) / sqrt(z);
  232. }
  233. if (x < 0.0)
  234. z = -z;
  235. return (z);
  236. }
  237. /* i1e() */
  238. double i1e(x) double x;
  239. {
  240. double y, z;
  241. z = fabs(x);
  242. if (z <= 8.0) {
  243. y = (z / 2.0) - 2.0;
  244. z = chbevl(y, A, 29) * z;
  245. } else {
  246. z = chbevl(32.0 / z - 2.0, B, 25) / sqrt(z);
  247. }
  248. if (x < 0.0)
  249. z = -z;
  250. return (z);
  251. }