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

   1: #ifndef lint
   2: static char sccsid[] = "@(#)j1.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 one
   9: 
  10: 	j1(x) returns the value of J1(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 J1 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: 	#6050 (20.98D)
  25: 	#6750 (19.19D)
  26: 	#7150 (19.35D)
  27: 
  28: 	y1(x) returns the value of Y1(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, j1.
  35: 
  36: 	The values of Y1 have not been checked
  37: 	to more than ten places.
  38: 
  39: 	Coefficients are from Hart & Cheney.
  40: 	#6447 (22.18D)
  41: 	#6750 (19.19D)
  42: 	#7150 (19.35D)
  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.581199354001606143928050809e21,
  56:     -.6672106568924916298020941484e20,
  57:     0.2316433580634002297931815435e19,
  58:     -.3588817569910106050743641413e17,
  59:     0.2908795263834775409737601689e15,
  60:     -.1322983480332126453125473247e13,
  61:     0.3413234182301700539091292655e10,
  62:     -.4695753530642995859767162166e7,
  63:     0.2701122710892323414856790990e4,
  64: };
  65: static double q1[] = {
  66:     0.1162398708003212287858529400e22,
  67:     0.1185770712190320999837113348e20,
  68:     0.6092061398917521746105196863e17,
  69:     0.2081661221307607351240184229e15,
  70:     0.5243710262167649715406728642e12,
  71:     0.1013863514358673989967045588e10,
  72:     0.1501793594998585505921097578e7,
  73:     0.1606931573481487801970916749e4,
  74:     1.0,
  75: };
  76: static double p2[] = {
  77:     -.4435757816794127857114720794e7,
  78:     -.9942246505077641195658377899e7,
  79:     -.6603373248364939109255245434e7,
  80:     -.1523529351181137383255105722e7,
  81:     -.1098240554345934672737413139e6,
  82:     -.1611616644324610116477412898e4,
  83:     0.0,
  84: };
  85: static double q2[] = {
  86:     -.4435757816794127856828016962e7,
  87:     -.9934124389934585658967556309e7,
  88:     -.6585339479723087072826915069e7,
  89:     -.1511809506634160881644546358e7,
  90:     -.1072638599110382011903063867e6,
  91:     -.1455009440190496182453565068e4,
  92:     1.0,
  93: };
  94: static double p3[] = {
  95:     0.3322091340985722351859704442e5,
  96:     0.8514516067533570196555001171e5,
  97:     0.6617883658127083517939992166e5,
  98:     0.1849426287322386679652009819e5,
  99:     0.1706375429020768002061283546e4,
 100:     0.3526513384663603218592175580e2,
 101:     0.0,
 102: };
 103: static double q3[] = {
 104:     0.7087128194102874357377502472e6,
 105:     0.1819458042243997298924553839e7,
 106:     0.1419460669603720892855755253e7,
 107:     0.4002944358226697511708610813e6,
 108:     0.3789022974577220264142952256e5,
 109:     0.8638367769604990967475517183e3,
 110:     1.0,
 111: };
 112: static double p4[] = {
 113:     -.9963753424306922225996744354e23,
 114:     0.2655473831434854326894248968e23,
 115:     -.1212297555414509577913561535e22,
 116:     0.2193107339917797592111427556e20,
 117:     -.1965887462722140658820322248e18,
 118:     0.9569930239921683481121552788e15,
 119:     -.2580681702194450950541426399e13,
 120:     0.3639488548124002058278999428e10,
 121:     -.2108847540133123652824139923e7,
 122:     0.0,
 123: };
 124: static double q4[] = {
 125:     0.5082067366941243245314424152e24,
 126:     0.5435310377188854170800653097e22,
 127:     0.2954987935897148674290758119e20,
 128:     0.1082258259408819552553850180e18,
 129:     0.2976632125647276729292742282e15,
 130:     0.6465340881265275571961681500e12,
 131:     0.1128686837169442121732366891e10,
 132:     0.1563282754899580604737366452e7,
 133:     0.1612361029677000859332072312e4,
 134:     1.0,
 135: };
 136: 
 137: double
 138: j1(arg) double arg;{
 139:     double xsq, n, d, x;
 140:     double sin(), cos(), sqrt();
 141:     int i;
 142: 
 143:     x = arg;
 144:     if(x < 0.) x = -x;
 145:     if(x > 8.){
 146:         asympt(x);
 147:         n = x - 3.*pio4;
 148:         n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
 149:         if(arg <0.) n = -n;
 150:         return(n);
 151:     }
 152:     xsq = x*x;
 153:     for(n=0,d=0,i=8;i>=0;i--){
 154:         n = n*xsq + p1[i];
 155:         d = d*xsq + q1[i];
 156:     }
 157:     return(arg*n/d);
 158: }
 159: 
 160: double
 161: y1(arg) double arg;{
 162:     double xsq, n, d, x;
 163:     double sin(), cos(), sqrt(), log(), j1();
 164:     int i;
 165: 
 166:     x = arg;
 167:     if(x <= 0.){
 168: #ifdef VAX
 169:         extern double infnan();
 170:         return(infnan(EDOM));       /* NaN */
 171: #else   /* IEEE double */
 172:         return(zero/zero);  /* IEEE machines: invalid operation */
 173: #endif
 174:     }
 175:     if(x > 8.){
 176:         asympt(x);
 177:         n = x - 3*pio4;
 178:         return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
 179:     }
 180:     xsq = x*x;
 181:     for(n=0,d=0,i=9;i>=0;i--){
 182:         n = n*xsq + p4[i];
 183:         d = d*xsq + q4[i];
 184:     }
 185:     return(x*n/d + tpi*(j1(x)*log(x)-1./x));
 186: }
 187: 
 188: static
 189: asympt(arg) double arg;{
 190:     double zsq, n, d;
 191:     int i;
 192:     zsq = 64./(arg*arg);
 193:     for(n=0,d=0,i=6;i>=0;i--){
 194:         n = n*zsq + p2[i];
 195:         d = d*zsq + q2[i];
 196:     }
 197:     pzero = n/d;
 198:     for(n=0,d=0,i=6;i>=0;i--){
 199:         n = n*zsq + p3[i];
 200:         d = d*zsq + q3[i];
 201:     }
 202:     qzero = (8./arg)*(n/d);
 203: }

Defined functions

Defined variables