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