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