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

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

Defined functions

Defined variables