1 /* 2 * Copyright (c) 1985 Regents of the University of California. 3 * All rights reserved. The Berkeley software License Agreement 4 * specifies the terms and conditions for redistribution. 5 */ 6 7 #ifndef lint 8 static char sccsid[] = "@(#)j0.c 5.3 (Berkeley) 09/22/88"; 9 #endif /* not lint */ 10 11 /* 12 floating point Bessel's function 13 of the first and second kinds 14 of order zero 15 16 j0(x) returns the value of J0(x) 17 for all real values of x. 18 19 There are no error returns. 20 Calls sin, cos, sqrt. 21 22 There is a niggling bug in J0 which 23 causes errors up to 2e-16 for x in the 24 interval [-8,8]. 25 The bug is caused by an inappropriate order 26 of summation of the series. rhm will fix it 27 someday. 28 29 Coefficients are from Hart & Cheney. 30 #5849 (19.22D) 31 #6549 (19.25D) 32 #6949 (19.41D) 33 34 y0(x) returns the value of Y0(x) 35 for positive real values of x. 36 For x<=0, if on the VAX, error number EDOM is set and 37 the reserved operand fault is generated; 38 otherwise (an IEEE machine) an invalid operation is performed. 39 40 Calls sin, cos, sqrt, log, j0. 41 42 The values of Y0 have not been checked 43 to more than ten places. 44 45 Coefficients are from Hart & Cheney. 46 #6245 (18.78D) 47 #6549 (19.25D) 48 #6949 (19.41D) 49 */ 50 51 #include "mathimpl.h" 52 53 #if defined(vax)||defined(tahoe) 54 #include <errno.h> 55 #else /* defined(vax)||defined(tahoe) */ 56 static const double zero = 0.e0; 57 #endif /* defined(vax)||defined(tahoe) */ 58 59 static double pzero, qzero; 60 61 static const double tpi = .6366197723675813430755350535e0; 62 static const double pio4 = .7853981633974483096156608458e0; 63 static const double p1[] = { 64 0.4933787251794133561816813446e21, 65 -.1179157629107610536038440800e21, 66 0.6382059341072356562289432465e19, 67 -.1367620353088171386865416609e18, 68 0.1434354939140344111664316553e16, 69 -.8085222034853793871199468171e13, 70 0.2507158285536881945555156435e11, 71 -.4050412371833132706360663322e8, 72 0.2685786856980014981415848441e5, 73 }; 74 static const double q1[] = { 75 0.4933787251794133562113278438e21, 76 0.5428918384092285160200195092e19, 77 0.3024635616709462698627330784e17, 78 0.1127756739679798507056031594e15, 79 0.3123043114941213172572469442e12, 80 0.6699987672982239671814028660e9, 81 0.1114636098462985378182402543e7, 82 0.1363063652328970604442810507e4, 83 1.0 84 }; 85 static const double p2[] = { 86 0.5393485083869438325262122897e7, 87 0.1233238476817638145232406055e8, 88 0.8413041456550439208464315611e7, 89 0.2016135283049983642487182349e7, 90 0.1539826532623911470917825993e6, 91 0.2485271928957404011288128951e4, 92 0.0, 93 }; 94 static const double q2[] = { 95 0.5393485083869438325560444960e7, 96 0.1233831022786324960844856182e8, 97 0.8426449050629797331554404810e7, 98 0.2025066801570134013891035236e7, 99 0.1560017276940030940592769933e6, 100 0.2615700736920839685159081813e4, 101 1.0, 102 }; 103 static const double p3[] = { 104 -.3984617357595222463506790588e4, 105 -.1038141698748464093880530341e5, 106 -.8239066313485606568803548860e4, 107 -.2365956170779108192723612816e4, 108 -.2262630641933704113967255053e3, 109 -.4887199395841261531199129300e1, 110 0.0, 111 }; 112 static const double q3[] = { 113 0.2550155108860942382983170882e6, 114 0.6667454239319826986004038103e6, 115 0.5332913634216897168722255057e6, 116 0.1560213206679291652539287109e6, 117 0.1570489191515395519392882766e5, 118 0.4087714673983499223402830260e3, 119 1.0, 120 }; 121 static const double p4[] = { 122 -.2750286678629109583701933175e20, 123 0.6587473275719554925999402049e20, 124 -.5247065581112764941297350814e19, 125 0.1375624316399344078571335453e18, 126 -.1648605817185729473122082537e16, 127 0.1025520859686394284509167421e14, 128 -.3436371222979040378171030138e11, 129 0.5915213465686889654273830069e8, 130 -.4137035497933148554125235152e5, 131 }; 132 static const double q4[] = { 133 0.3726458838986165881989980e21, 134 0.4192417043410839973904769661e19, 135 0.2392883043499781857439356652e17, 136 0.9162038034075185262489147968e14, 137 0.2613065755041081249568482092e12, 138 0.5795122640700729537480087915e9, 139 0.1001702641288906265666651753e7, 140 0.1282452772478993804176329391e4, 141 1.0, 142 }; 143 144 static void asympt(); 145 146 double 147 j0(arg) double arg;{ 148 double argsq, n, d; 149 int i; 150 151 if(arg < 0.) arg = -arg; 152 if(arg > 8.){ 153 asympt(arg); 154 n = arg - pio4; 155 return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n))); 156 } 157 argsq = arg*arg; 158 for(n=0,d=0,i=8;i>=0;i--){ 159 n = n*argsq + p1[i]; 160 d = d*argsq + q1[i]; 161 } 162 return(n/d); 163 } 164 165 double 166 y0(arg) double arg;{ 167 double argsq, n, d; 168 int i; 169 170 if(arg <= 0.){ 171 #if defined(vax)||defined(tahoe) 172 return(infnan(EDOM)); /* NaN */ 173 #else /* defined(vax)||defined(tahoe) */ 174 return(zero/zero); /* IEEE machines: invalid operation */ 175 #endif /* defined(vax)||defined(tahoe) */ 176 } 177 if(arg > 8.){ 178 asympt(arg); 179 n = arg - pio4; 180 return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n))); 181 } 182 argsq = arg*arg; 183 for(n=0,d=0,i=8;i>=0;i--){ 184 n = n*argsq + p4[i]; 185 d = d*argsq + q4[i]; 186 } 187 return(n/d + tpi*j0(arg)*log(arg)); 188 } 189 190 static void 191 asympt(arg) double arg;{ 192 double zsq, n, d; 193 int i; 194 zsq = 64./(arg*arg); 195 for(n=0,d=0,i=6;i>=0;i--){ 196 n = n*zsq + p2[i]; 197 d = d*zsq + q2[i]; 198 } 199 pzero = n/d; 200 for(n=0,d=0,i=6;i>=0;i--){ 201 n = n*zsq + p3[i]; 202 d = d*zsq + q3[i]; 203 } 204 qzero = (8./arg)*(n/d); 205 } 206