1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 /* 31 * Floating point Bessel's function of the first and second kinds 32 * of order zero: j0(x),y0(x); 33 * 34 * Special cases: 35 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; 36 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. 37 */ 38 39 #pragma weak j0 = __j0 40 #pragma weak y0 = __y0 41 42 #include "libm.h" 43 #include "libm_synonyms.h" 44 #include "libm_protos.h" 45 #include <math.h> 46 #include <values.h> 47 48 #define GENERIC double 49 static const GENERIC 50 zero = 0.0, 51 small = 1.0e-5, 52 tiny = 1.0e-18, 53 one = 1.0, 54 eight = 8.0, 55 invsqrtpi = 5.641895835477562869480794515607725858441e-0001, 56 tpi = 0.636619772367581343075535053490057448; 57 58 static GENERIC pzero(GENERIC), qzero(GENERIC); 59 static const GENERIC r0[4] = { /* [1.e-5, 1.28] */ 60 -2.500000000000003622131880894830476755537e-0001, 61 1.095597547334830263234433855932375353303e-0002, 62 -1.819734750463320921799187258987098087697e-0004, 63 9.977001946806131657544212501069893930846e-0007, 64 }; 65 static const GENERIC s0[4] = { /* [1.e-5, 1.28] */ 66 1.0, 67 1.867609810662950169966782360588199673741e-0002, 68 1.590389206181565490878430827706972074208e-0004, 69 6.520867386742583632375520147714499522721e-0007, 70 }; 71 static const GENERIC r1[9] = { /* [1.28,8] */ 72 9.999999999999999942156495584397047660949e-0001, 73 -2.389887722731319130476839836908143731281e-0001, 74 1.293359476138939027791270393439493640570e-0002, 75 -2.770985642343140122168852400228563364082e-0004, 76 2.905241575772067678086738389169625218912e-0006, 77 -1.636846356264052597969042009265043251279e-0008, 78 5.072306160724884775085431059052611737827e-0011, 79 -8.187060730684066824228914775146536139112e-0014, 80 5.422219326959949863954297860723723423842e-0017, 81 }; 82 static const GENERIC s1[9] = { /* [1.28,8] */ 83 1.0, 84 1.101122772686807702762104741932076228349e-0002, 85 6.140169310641649223411427764669143978228e-0005, 86 2.292035877515152097976946119293215705250e-0007, 87 6.356910426504644334558832036362219583789e-0010, 88 1.366626326900219555045096999553948891401e-0012, 89 2.280399586866739522891837985560481180088e-0015, 90 2.801559820648939665270492520004836611187e-0018, 91 2.073101088320349159764410261466350732968e-0021, 92 }; 93 94 GENERIC 95 j0(GENERIC x) { 96 GENERIC z, s, c, ss, cc, r, u, v, ox; 97 int i; 98 99 if (isnan(x)) 100 return (x*x); /* + -> * for Cheetah */ 101 ox = x; 102 x = fabs(x); 103 if (x > 8.0) { 104 if (!finite(x)) 105 return (zero); 106 s = sin(x); 107 c = cos(x); 108 /* 109 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) 110 * where x0 = x-pi/4 111 * Better formula: 112 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) 113 * = 1/sqrt(2) * (cos(x) + sin(x)) 114 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) 115 * = 1/sqrt(2) * (sin(x) - cos(x)) 116 * To avoid cancellation, use 117 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 118 * to compute the worse one. 119 */ 120 if (x > 8.9e307) { /* x+x may overflow */ 121 ss = s-c; 122 cc = s+c; 123 } else if (signbit(s) != signbit(c)) { 124 ss = s - c; 125 cc = -cos(x+x)/ss; 126 } else { 127 cc = s + c; 128 ss = -cos(x+x)/cc; 129 } 130 /* 131 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 132 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 133 */ 134 if (x > 1.0e40) z = (invsqrtpi*cc)/sqrt(x); 135 else { 136 u = pzero(x); v = qzero(x); 137 z = invsqrtpi*(u*cc-v*ss)/sqrt(x); 138 } 139 /* force to pass SVR4 even the result is wrong (sign) */ 140 if (x > X_TLOSS) 141 return (_SVID_libm_err(ox, z, 34)); 142 else 143 return (z); 144 } 145 if (x <= small) { 146 if (x <= tiny) 147 return (one-x); 148 else 149 return (one-x*x*0.25); 150 } 151 z = x*x; 152 if (x <= 1.28) { 153 r = r0[0]+z*(r0[1]+z*(r0[2]+z*r0[3])); 154 s = s0[0]+z*(s0[1]+z*(s0[2]+z*s0[3])); 155 return (one + z*(r/s)); 156 } else { 157 for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) { 158 r = r*z + r1[i]; 159 s = s*z + s1[i]; 160 } 161 return (r/s); 162 } 163 } 164 165 static const GENERIC u0[13] = { 166 -7.380429510868722526754723020704317641941e-0002, 167 1.772607102684869924301459663049874294814e-0001, 168 -1.524370666542713828604078090970799356306e-0002, 169 4.650819100693891757143771557629924591915e-0004, 170 -7.125768872339528975036316108718239946022e-0006, 171 6.411017001656104598327565004771515257146e-0008, 172 -3.694275157433032553021246812379258781665e-0010, 173 1.434364544206266624252820889648445263842e-0012, 174 -3.852064731859936455895036286874139896861e-0015, 175 7.182052899726138381739945881914874579696e-0018, 176 -9.060556574619677567323741194079797987200e-0021, 177 7.124435467408860515265552217131230511455e-0024, 178 -2.709726774636397615328813121715432044771e-0027, 179 }; 180 static const GENERIC v0[5] = { 181 1.0, 182 4.678678931512549002587702477349214886475e-0003, 183 9.486828955529948534822800829497565178985e-0006, 184 1.001495929158861646659010844136682454906e-0008, 185 4.725338116256021660204443235685358593611e-0012, 186 }; 187 188 GENERIC 189 y0(GENERIC x) { 190 GENERIC z, /* d, */ s, c, ss, cc, u, v; 191 int i; 192 193 if (isnan(x)) 194 return (x*x); /* + -> * for Cheetah */ 195 if (x <= zero) { 196 if (x == zero) 197 /* d= -one/(x-x); */ 198 return (_SVID_libm_err(x, x, 8)); 199 else 200 /* d = zero/(x-x); */ 201 return (_SVID_libm_err(x, x, 9)); 202 } 203 if (x > 8.0) { 204 if (!finite(x)) 205 return (zero); 206 s = sin(x); 207 c = cos(x); 208 /* 209 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) 210 * where x0 = x-pi/4 211 * Better formula: 212 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) 213 * = 1/sqrt(2) * (cos(x) + sin(x)) 214 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) 215 * = 1/sqrt(2) * (sin(x) - cos(x)) 216 * To avoid cancellation, use 217 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 218 * to compute the worse one. 219 */ 220 if (x > 8.9e307) { /* x+x may overflow */ 221 ss = s-c; 222 cc = s+c; 223 } else if (signbit(s) != signbit(c)) { 224 ss = s - c; 225 cc = -cos(x+x)/ss; 226 } else { 227 cc = s + c; 228 ss = -cos(x+x)/cc; 229 } 230 /* 231 * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss) 232 * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc) 233 */ 234 if (x > 1.0e40) 235 z = (invsqrtpi*ss)/sqrt(x); 236 else 237 z = invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x); 238 if (x > X_TLOSS) 239 return (_SVID_libm_err(x, z, 35)); 240 else 241 return (z); 242 243 } 244 if (x <= tiny) { 245 return (u0[0] + tpi*log(x)); 246 } 247 z = x*x; 248 for (u = u0[12], i = 11; i >= 0; i--) u = u*z + u0[i]; 249 v = v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4]))); 250 return (u/v + tpi*(j0(x)*log(x))); 251 } 252 253 static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */ 254 .4861344183386052721391238447e5, 255 .1377662549407112278133438945e6, 256 .1222466364088289731869114004e6, 257 .4107070084315176135583353374e5, 258 .5026073801860637125889039915e4, 259 .1783193659125479654541542419e3, 260 .88010344055383421691677564e0, 261 }; 262 static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */ 263 .4861344183386052721414037058e5, 264 .1378196632630384670477582699e6, 265 .1223967185341006542748936787e6, 266 .4120150243795353639995862617e5, 267 .5068271181053546392490184353e4, 268 .1829817905472769960535671664e3, 269 1.0, 270 }; 271 static const GENERIC huge = 1.0e10; 272 273 static GENERIC 274 pzero(GENERIC x) { 275 GENERIC s, r, t, z; 276 int i; 277 if (x > huge) 278 return (one); 279 t = eight/x; z = t*t; 280 r = pr[5]+z*pr[6]; 281 s = ps[5]+z; 282 for (i = 4; i >= 0; i--) { 283 r = r*z + pr[i]; 284 s = s*z + ps[i]; 285 } 286 return (r/s); 287 } 288 289 static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */ 290 -.1731210995701068539185611951e3, 291 -.5522559165936166961235240613e3, 292 -.5604935606637346590614529613e3, 293 -.2200430300226009379477365011e3, 294 -.323869355375648849771296746e2, 295 -.14294979207907956223499258e1, 296 -.834690374102384988158918e-2, 297 }; 298 static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */ 299 .1107975037248683865326709645e5, 300 .3544581680627082674651471873e5, 301 .3619118937918394132179019059e5, 302 .1439895563565398007471485822e5, 303 .2190277023344363955930226234e4, 304 .106695157020407986137501682e3, 305 1.0, 306 }; 307 308 static GENERIC 309 qzero(GENERIC x) { 310 GENERIC s, r, t, z; 311 int i; 312 if (x > huge) 313 return (-0.125/x); 314 t = eight/x; z = t*t; 315 r = qr[5]+z*qr[6]; 316 s = qs[5]+z; 317 for (i = 4; i >= 0; i--) { 318 r = r*z + qr[i]; 319 s = s*z + qs[i]; 320 } 321 return (t*(r/s)); 322 } 323