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: j1(x),y1(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 j1 = __j1 40 #pragma weak y1 = __y1 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-20, 53 one = 1.0, 54 invsqrtpi = 5.641895835477562869480794515607725858441e-0001, 55 tpi = 0.636619772367581343075535053490057448; 56 57 static GENERIC pone(GENERIC), qone(GENERIC); 58 static const GENERIC r0[4] = { 59 -6.250000000000002203053200981413218949548e-0002, 60 1.600998455640072901321605101981501263762e-0003, 61 -1.963888815948313758552511884390162864930e-0005, 62 8.263917341093549759781339713418201620998e-0008, 63 }; 64 static const GENERIC s0[7] = { 65 1.0e0, 66 1.605069137643004242395356851797873766927e-0002, 67 1.149454623251299996428500249509098499383e-0004, 68 3.849701673735260970379681807910852327825e-0007, 69 }; 70 static const GENERIC r1[12] = { 71 4.999999999999999995517408894340485471724e-0001, 72 -6.003825028120475684835384519945468075423e-0002, 73 2.301719899263321828388344461995355419832e-0003, 74 -4.208494869238892934859525221654040304068e-0005, 75 4.377745135188837783031540029700282443388e-0007, 76 -2.854106755678624335145364226735677754179e-0009, 77 1.234002865443952024332943901323798413689e-0011, 78 -3.645498437039791058951273508838177134310e-0014, 79 7.404320596071797459925377103787837414422e-0017, 80 -1.009457448277522275262808398517024439084e-0019, 81 8.520158355824819796968771418801019930585e-0023, 82 -3.458159926081163274483854614601091361424e-0026, 83 }; 84 static const GENERIC s1[5] = { 85 1.0e0, 86 4.923499437590484879081138588998986303306e-0003, 87 1.054389489212184156499666953501976688452e-0005, 88 1.180768373106166527048240364872043816050e-0008, 89 5.942665743476099355323245707680648588540e-0012, 90 }; 91 92 GENERIC 93 j1(GENERIC x) { 94 GENERIC z, d, s, c, ss, cc, r; 95 int i, sgn; 96 97 if (!finite(x)) 98 return (one/x); 99 sgn = signbit(x); 100 x = fabs(x); 101 if (x > 8.00) { 102 s = sin(x); 103 c = cos(x); 104 /* 105 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0)) 106 * where x0 = x-3pi/4 107 * Better formula: 108 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) 109 * = 1/sqrt(2) * (sin(x) - cos(x)) 110 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 111 * = -1/sqrt(2) * (cos(x) + sin(x)) 112 * To avoid cancellation, use 113 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 114 * to compute the worse one. 115 */ 116 if (x > 8.9e307) { /* x+x may overflow */ 117 ss = -s-c; 118 cc = s-c; 119 } else if (signbit(s) != signbit(c)) { 120 cc = s - c; 121 ss = cos(x+x)/cc; 122 } else { 123 ss = -s-c; 124 cc = cos(x+x)/ss; 125 } 126 /* 127 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss) 128 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc) 129 */ 130 if (x > 1.0e40) 131 d = (invsqrtpi*cc)/sqrt(x); 132 else 133 d = invsqrtpi*(pone(x)*cc-qone(x)*ss)/sqrt(x); 134 135 if (x > X_TLOSS) { 136 if (sgn != 0) { d = -d; x = -x; } 137 return (_SVID_libm_err(x, d, 36)); 138 } else 139 if (sgn == 0) 140 return (d); 141 else 142 return (-d); 143 } 144 if (x <= small) { 145 if (x <= tiny) 146 d = 0.5*x; 147 else 148 d = x*(0.5-x*x*0.125); 149 if (sgn == 0) 150 return (d); 151 else 152 return (-d); 153 } 154 z = x*x; 155 if (x < 1.28) { 156 r = r0[3]; 157 s = s0[3]; 158 for (i = 2; i >= 0; i--) { 159 r = r*z + r0[i]; 160 s = s*z + s0[i]; 161 } 162 d = x*0.5+x*(z*(r/s)); 163 } else { 164 r = r1[11]; 165 for (i = 10; i >= 0; i--) r = r*z + r1[i]; 166 s = s1[0]+z*(s1[1]+z*(s1[2]+z*(s1[3]+z*s1[4]))); 167 d = x*(r/s); 168 } 169 if (sgn == 0) 170 return (d); 171 else 172 return (-d); 173 } 174 175 static const GENERIC u0[4] = { 176 -1.960570906462389461018983259589655961560e-0001, 177 4.931824118350661953459180060007970291139e-0002, 178 -1.626975871565393656845930125424683008677e-0003, 179 1.359657517926394132692884168082224258360e-0005, 180 }; 181 static const GENERIC v0[5] = { 182 1.0e0, 183 2.565807214838390835108224713630901653793e-0002, 184 3.374175208978404268650522752520906231508e-0004, 185 2.840368571306070719539936935220728843177e-0006, 186 1.396387402048998277638900944415752207592e-0008, 187 }; 188 static const GENERIC u1[12] = { 189 -1.960570906462389473336339614647555351626e-0001, 190 5.336268030335074494231369159933012844735e-0002, 191 -2.684137504382748094149184541866332033280e-0003, 192 5.737671618979185736981543498580051903060e-0005, 193 -6.642696350686335339171171785557663224892e-0007, 194 4.692417922568160354012347591960362101664e-0009, 195 -2.161728635907789319335231338621412258355e-0011, 196 6.727353419738316107197644431844194668702e-0014, 197 -1.427502986803861372125234355906790573422e-0016, 198 2.020392498726806769468143219616642940371e-0019, 199 -1.761371948595104156753045457888272716340e-0022, 200 7.352828391941157905175042420249225115816e-0026, 201 }; 202 static const GENERIC v1[5] = { 203 1.0e0, 204 5.029187436727947764916247076102283399442e-0003, 205 1.102693095808242775074856548927801750627e-0005, 206 1.268035774543174837829534603830227216291e-0008, 207 6.579416271766610825192542295821308730206e-0012, 208 }; 209 210 211 GENERIC 212 y1(GENERIC x) { 213 GENERIC z, d, s, c, ss, cc, u, v; 214 int i; 215 216 if (isnan(x)) 217 return (x*x); /* + -> * for Cheetah */ 218 if (x <= zero) { 219 if (x == zero) 220 /* return -one/zero; */ 221 return (_SVID_libm_err(x, x, 10)); 222 else 223 /* return zero/zero; */ 224 return (_SVID_libm_err(x, x, 11)); 225 } 226 if (x > 8.0) { 227 if (!finite(x)) 228 return (zero); 229 s = sin(x); 230 c = cos(x); 231 /* 232 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0)) 233 * where x0 = x-3pi/4 234 * Better formula: 235 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) 236 * = 1/sqrt(2) * (sin(x) - cos(x)) 237 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 238 * = -1/sqrt(2) * (cos(x) + sin(x)) 239 * To avoid cancellation, use 240 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 241 * to compute the worse one. 242 */ 243 if (x > 8.9e307) { /* x+x may overflow */ 244 ss = -s-c; 245 cc = s-c; 246 } else if (signbit(s) != signbit(c)) { 247 cc = s - c; 248 ss = cos(x+x)/cc; 249 } else { 250 ss = -s-c; 251 cc = cos(x+x)/ss; 252 } 253 /* 254 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss) 255 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc) 256 */ 257 if (x > 1.0e91) 258 d = (invsqrtpi*ss)/sqrt(x); 259 else 260 d = invsqrtpi*(pone(x)*ss+qone(x)*cc)/sqrt(x); 261 262 if (x > X_TLOSS) 263 return (_SVID_libm_err(x, d, 37)); 264 else 265 return (d); 266 } 267 if (x <= tiny) { 268 return (-tpi/x); 269 } 270 z = x*x; 271 if (x < 1.28) { 272 u = u0[3]; v = v0[3]+z*v0[4]; 273 for (i = 2; i >= 0; i--) { 274 u = u*z + u0[i]; 275 v = v*z + v0[i]; 276 } 277 } else { 278 for (u = u1[11], i = 10; i >= 0; i--) u = u*z+u1[i]; 279 v = v1[0]+z*(v1[1]+z*(v1[2]+z*(v1[3]+z*v1[4]))); 280 } 281 return (x*(u/v) + tpi*(j1(x)*log(x)-one/x)); 282 } 283 284 static const GENERIC pr0[6] = { 285 -.4435757816794127857114720794e7, 286 -.9942246505077641195658377899e7, 287 -.6603373248364939109255245434e7, 288 -.1523529351181137383255105722e7, 289 -.1098240554345934672737413139e6, 290 -.1611616644324610116477412898e4, 291 }; 292 static const GENERIC ps0[6] = { 293 -.4435757816794127856828016962e7, 294 -.9934124389934585658967556309e7, 295 -.6585339479723087072826915069e7, 296 -.1511809506634160881644546358e7, 297 -.1072638599110382011903063867e6, 298 -.1455009440190496182453565068e4, 299 }; 300 static const GENERIC huge = 1.0e10; 301 302 static GENERIC 303 pone(GENERIC x) { 304 GENERIC s, r, t, z; 305 int i; 306 /* assume x > 8 */ 307 if (x > huge) 308 return (one); 309 310 t = 8.0/x; z = t*t; 311 r = pr0[5]; s = ps0[5]+z; 312 for (i = 4; i >= 0; i--) { 313 r = z*r + pr0[i]; 314 s = z*s + ps0[i]; 315 } 316 return (r/s); 317 } 318 319 320 static const GENERIC qr0[6] = { 321 0.3322091340985722351859704442e5, 322 0.8514516067533570196555001171e5, 323 0.6617883658127083517939992166e5, 324 0.1849426287322386679652009819e5, 325 0.1706375429020768002061283546e4, 326 0.3526513384663603218592175580e2, 327 }; 328 static const GENERIC qs0[6] = { 329 0.7087128194102874357377502472e6, 330 0.1819458042243997298924553839e7, 331 0.1419460669603720892855755253e7, 332 0.4002944358226697511708610813e6, 333 0.3789022974577220264142952256e5, 334 0.8638367769604990967475517183e3, 335 }; 336 337 static GENERIC 338 qone(GENERIC x) { 339 GENERIC s, r, t, z; 340 int i; 341 if (x > huge) 342 return (0.375/x); 343 344 t = 8.0/x; z = t*t; 345 /* assume x > 8 */ 346 r = qr0[5]; s = qs0[5]+z; 347 for (i = 4; i >= 0; i--) { 348 r = z*r + qr0[i]; 349 s = z*s + qs0[i]; 350 } 351 return (t*(r/s)); 352 } 353