4.3BSD - /usr/src/usr.lib/libm/j0.c

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

Defined functions

Defined variables