2.9BSD - /usr/src/lib/m/j1.c

   1: /*
   2: 	floating point Bessel's function
   3: 	of the first and second kinds
   4: 	of order one
   5: 
   6: 	j1(x) returns the value of J1(x)
   7: 	for all real values of x.
   8: 
   9: 	There are no error returns.
  10: 	Calls sin, cos, sqrt.
  11: 
  12: 	There is a niggling bug in J1 which
  13: 	causes errors up to 2e-16 for x in the
  14: 	interval [-8,8].
  15: 	The bug is caused by an inappropriate order
  16: 	of summation of the series.  rhm will fix it
  17: 	someday.
  18: 
  19: 	Coefficients are from Hart & Cheney.
  20: 	#6050 (20.98D)
  21: 	#6750 (19.19D)
  22: 	#7150 (19.35D)
  23: 
  24: 	y1(x) returns the value of Y1(x)
  25: 	for positive real values of x.
  26: 	For x<=0, error number EDOM is set and a
  27: 	large negative value is returned.
  28: 
  29: 	Calls sin, cos, sqrt, log, j1.
  30: 
  31: 	The values of Y1 have not been checked
  32: 	to more than ten places.
  33: 
  34: 	Coefficients are from Hart & Cheney.
  35: 	#6447 (22.18D)
  36: 	#6750 (19.19D)
  37: 	#7150 (19.35D)
  38: */
  39: 
  40: #include <math.h>
  41: #include <errno.h>
  42: 
  43: int errno;
  44: static double pzero, qzero;
  45: static double tpi   = .6366197723675813430755350535e0;
  46: static double pio4  = .7853981633974483096156608458e0;
  47: static double p1[] = {
  48:     0.581199354001606143928050809e21,
  49:     -.6672106568924916298020941484e20,
  50:     0.2316433580634002297931815435e19,
  51:     -.3588817569910106050743641413e17,
  52:     0.2908795263834775409737601689e15,
  53:     -.1322983480332126453125473247e13,
  54:     0.3413234182301700539091292655e10,
  55:     -.4695753530642995859767162166e7,
  56:     0.2701122710892323414856790990e4,
  57: };
  58: static double q1[] = {
  59:     0.1162398708003212287858529400e22,
  60:     0.1185770712190320999837113348e20,
  61:     0.6092061398917521746105196863e17,
  62:     0.2081661221307607351240184229e15,
  63:     0.5243710262167649715406728642e12,
  64:     0.1013863514358673989967045588e10,
  65:     0.1501793594998585505921097578e7,
  66:     0.1606931573481487801970916749e4,
  67:     1.0,
  68: };
  69: static double p2[] = {
  70:     -.4435757816794127857114720794e7,
  71:     -.9942246505077641195658377899e7,
  72:     -.6603373248364939109255245434e7,
  73:     -.1523529351181137383255105722e7,
  74:     -.1098240554345934672737413139e6,
  75:     -.1611616644324610116477412898e4,
  76:     0.0,
  77: };
  78: static double q2[] = {
  79:     -.4435757816794127856828016962e7,
  80:     -.9934124389934585658967556309e7,
  81:     -.6585339479723087072826915069e7,
  82:     -.1511809506634160881644546358e7,
  83:     -.1072638599110382011903063867e6,
  84:     -.1455009440190496182453565068e4,
  85:     1.0,
  86: };
  87: static double p3[] = {
  88:     0.3322091340985722351859704442e5,
  89:     0.8514516067533570196555001171e5,
  90:     0.6617883658127083517939992166e5,
  91:     0.1849426287322386679652009819e5,
  92:     0.1706375429020768002061283546e4,
  93:     0.3526513384663603218592175580e2,
  94:     0.0,
  95: };
  96: static double q3[] = {
  97:     0.7087128194102874357377502472e6,
  98:     0.1819458042243997298924553839e7,
  99:     0.1419460669603720892855755253e7,
 100:     0.4002944358226697511708610813e6,
 101:     0.3789022974577220264142952256e5,
 102:     0.8638367769604990967475517183e3,
 103:     1.0,
 104: };
 105: static double p4[] = {
 106:     -.9963753424306922225996744354e23,
 107:     0.2655473831434854326894248968e23,
 108:     -.1212297555414509577913561535e22,
 109:     0.2193107339917797592111427556e20,
 110:     -.1965887462722140658820322248e18,
 111:     0.9569930239921683481121552788e15,
 112:     -.2580681702194450950541426399e13,
 113:     0.3639488548124002058278999428e10,
 114:     -.2108847540133123652824139923e7,
 115:     0.0,
 116: };
 117: static double q4[] = {
 118:     0.5082067366941243245314424152e24,
 119:     0.5435310377188854170800653097e22,
 120:     0.2954987935897148674290758119e20,
 121:     0.1082258259408819552553850180e18,
 122:     0.2976632125647276729292742282e15,
 123:     0.6465340881265275571961681500e12,
 124:     0.1128686837169442121732366891e10,
 125:     0.1563282754899580604737366452e7,
 126:     0.1612361029677000859332072312e4,
 127:     1.0,
 128: };
 129: 
 130: double
 131: j1(arg) double arg;{
 132:     double xsq, n, d, x;
 133:     double sin(), cos(), sqrt();
 134:     int i;
 135: 
 136:     x = arg;
 137:     if(x < 0.) x = -x;
 138:     if(x > 8.){
 139:         asympt(x);
 140:         n = x - 3.*pio4;
 141:         n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
 142:         if(arg <0.) n = -n;
 143:         return(n);
 144:     }
 145:     xsq = x*x;
 146:     for(n=0,d=0,i=8;i>=0;i--){
 147:         n = n*xsq + p1[i];
 148:         d = d*xsq + q1[i];
 149:     }
 150:     return(arg*n/d);
 151: }
 152: 
 153: double
 154: y1(arg) double arg;{
 155:     double xsq, n, d, x;
 156:     double sin(), cos(), sqrt(), log(), j1();
 157:     int i;
 158: 
 159:     errno = 0;
 160:     x = arg;
 161:     if(x <= 0.){
 162:         errno = EDOM;
 163:         return(-HUGE);
 164:     }
 165:     if(x > 8.){
 166:         asympt(x);
 167:         n = x - 3*pio4;
 168:         return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
 169:     }
 170:     xsq = x*x;
 171:     for(n=0,d=0,i=9;i>=0;i--){
 172:         n = n*xsq + p4[i];
 173:         d = d*xsq + q4[i];
 174:     }
 175:     return(x*n/d + tpi*(j1(x)*log(x)-1./x));
 176: }
 177: 
 178: static
 179: asympt(arg) double arg;{
 180:     double zsq, n, d;
 181:     int i;
 182:     zsq = 64./(arg*arg);
 183:     for(n=0,d=0,i=6;i>=0;i--){
 184:         n = n*zsq + p2[i];
 185:         d = d*zsq + q2[i];
 186:     }
 187:     pzero = n/d;
 188:     for(n=0,d=0,i=6;i>=0;i--){
 189:         n = n*zsq + p3[i];
 190:         d = d*zsq + q3[i];
 191:     }
 192:     qzero = (8./arg)*(n/d);
 193: }

Defined functions

Defined variables