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