Print this page
11210 libm should be cstyle(1ONBLD) clean


   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_protos.h"
  44 #include <math.h>
  45 #include <values.h>
  46 
  47 #define GENERIC double
  48 static const GENERIC
  49 zero    = 0.0,
  50 small   = 1.0e-5,
  51 tiny    = 1.0e-20,
  52 one     = 1.0,
  53 invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
  54 tpi     = 0.636619772367581343075535053490057448;
  55 
  56 static GENERIC pone(GENERIC), qone(GENERIC);









  57 static const GENERIC r0[4] = {
  58         -6.250000000000002203053200981413218949548e-0002,
  59         1.600998455640072901321605101981501263762e-0003,
  60         -1.963888815948313758552511884390162864930e-0005,
  61         8.263917341093549759781339713418201620998e-0008,
  62 };

  63 static const GENERIC s0[7] = {
  64         1.0e0,
  65         1.605069137643004242395356851797873766927e-0002,
  66         1.149454623251299996428500249509098499383e-0004,
  67         3.849701673735260970379681807910852327825e-0007,
  68 };

  69 static const GENERIC r1[12] = {
  70         4.999999999999999995517408894340485471724e-0001,
  71         -6.003825028120475684835384519945468075423e-0002,
  72         2.301719899263321828388344461995355419832e-0003,
  73         -4.208494869238892934859525221654040304068e-0005,
  74         4.377745135188837783031540029700282443388e-0007,
  75         -2.854106755678624335145364226735677754179e-0009,
  76         1.234002865443952024332943901323798413689e-0011,
  77         -3.645498437039791058951273508838177134310e-0014,
  78         7.404320596071797459925377103787837414422e-0017,
  79         -1.009457448277522275262808398517024439084e-0019,
  80         8.520158355824819796968771418801019930585e-0023,
  81         -3.458159926081163274483854614601091361424e-0026,
  82 };

  83 static const GENERIC s1[5] = {
  84         1.0e0,
  85         4.923499437590484879081138588998986303306e-0003,
  86         1.054389489212184156499666953501976688452e-0005,
  87         1.180768373106166527048240364872043816050e-0008,
  88         5.942665743476099355323245707680648588540e-0012,
  89 };
  90 
  91 GENERIC
  92 j1(GENERIC x) {

  93         GENERIC z, d, s, c, ss, cc, r;
  94         int i, sgn;
  95 
  96         if (!finite(x))
  97                 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         /*
 104          * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
 105          * where x0 = x-3pi/4
 106          *      Better formula:
 107          *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 108          *                      =  1/sqrt(2) * (sin(x) - cos(x))
 109          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 110          *                      = -1/sqrt(2) * (cos(x) + sin(x))
 111          * To avoid cancellation, use
 112          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 113          * to compute the worse one.
 114          */

 115                 if (x > 8.9e307) {   /* x+x may overflow */
 116                         ss = -s-c;
 117                         cc =  s-c;
 118                 } else if (signbit(s) != signbit(c)) {
 119                         cc = s - c;
 120                         ss = cos(x+x)/cc;
 121                 } else {
 122                         ss = -s-c;
 123                         cc = cos(x+x)/ss;
 124                 }

 125         /*
 126          * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
 127          * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
 128          */
 129                 if (x > 1.0e40)
 130                     d = (invsqrtpi*cc)/sqrt(x);
 131                 else
 132                         d =  invsqrtpi*(pone(x)*cc-qone(x)*ss)/sqrt(x);
 133 
 134                 if (x > X_TLOSS) {
 135                     if (sgn != 0) { d = -d; x = -x; }




 136                         return (_SVID_libm_err(x, d, 36));
 137                 } else
 138                     if (sgn == 0)
 139                                 return (d);
 140                         else
 141                                 return (-d);
 142         }


 143         if (x <= small) {
 144                 if (x <= tiny)
 145                         d = 0.5*x;
 146                 else
 147                         d =  x*(0.5-x*x*0.125);

 148                 if (sgn == 0)
 149                         return (d);
 150                 else
 151                         return (-d);
 152         }
 153         z = x*x;


 154         if (x < 1.28) {
 155             r = r0[3];
 156             s = s0[3];

 157             for (i = 2; i >= 0; i--) {
 158                 r = r*z + r0[i];
 159                 s = s*z + s0[i];
 160             }
 161             d = x*0.5+x*(z*(r/s));

 162         } else {
 163             r = r1[11];
 164             for (i = 10; i >= 0; i--) r = r*z + r1[i];
 165             s = s1[0]+z*(s1[1]+z*(s1[2]+z*(s1[3]+z*s1[4])));
 166             d = x*(r/s);



 167         }

 168         if (sgn == 0)
 169                 return (d);
 170         else
 171                 return (-d);
 172 }
 173 
 174 static const GENERIC u0[4] = {
 175         -1.960570906462389461018983259589655961560e-0001,
 176         4.931824118350661953459180060007970291139e-0002,
 177         -1.626975871565393656845930125424683008677e-0003,
 178         1.359657517926394132692884168082224258360e-0005,
 179 };

 180 static const GENERIC v0[5] = {
 181         1.0e0,
 182         2.565807214838390835108224713630901653793e-0002,
 183         3.374175208978404268650522752520906231508e-0004,
 184         2.840368571306070719539936935220728843177e-0006,
 185         1.396387402048998277638900944415752207592e-0008,
 186 };

 187 static const GENERIC u1[12] = {
 188         -1.960570906462389473336339614647555351626e-0001,
 189         5.336268030335074494231369159933012844735e-0002,
 190         -2.684137504382748094149184541866332033280e-0003,
 191         5.737671618979185736981543498580051903060e-0005,
 192         -6.642696350686335339171171785557663224892e-0007,
 193         4.692417922568160354012347591960362101664e-0009,
 194         -2.161728635907789319335231338621412258355e-0011,
 195         6.727353419738316107197644431844194668702e-0014,
 196         -1.427502986803861372125234355906790573422e-0016,
 197         2.020392498726806769468143219616642940371e-0019,
 198         -1.761371948595104156753045457888272716340e-0022,
 199         7.352828391941157905175042420249225115816e-0026,
 200 };

 201 static const GENERIC v1[5] = {
 202         1.0e0,
 203         5.029187436727947764916247076102283399442e-0003,
 204         1.102693095808242775074856548927801750627e-0005,
 205         1.268035774543174837829534603830227216291e-0008,
 206         6.579416271766610825192542295821308730206e-0012,
 207 };
 208 
 209 
 210 GENERIC
 211 y1(GENERIC x) {

 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         if (x <= zero) {
 218                 if (x == zero)
 219                     /* return -one/zero;  */
 220                     return (_SVID_libm_err(x, x, 10));
 221                 else
 222                     /* return zero/zero; */
 223                     return (_SVID_libm_err(x, x, 11));
 224         }

 225         if (x > 8.0) {
 226                 if (!finite(x))
 227                         return (zero);

 228                 s = sin(x);
 229                 c = cos(x);


 230         /*
 231          * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
 232          * where x0 = x-3pi/4
 233          *      Better formula:
 234          *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 235          *                      =  1/sqrt(2) * (sin(x) - cos(x))
 236          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 237          *                      = -1/sqrt(2) * (cos(x) + sin(x))
 238          * To avoid cancellation, use
 239          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 240          * to compute the worse one.
 241          */

 242                 if (x > 8.9e307) {   /* x+x may overflow */
 243                         ss = -s-c;
 244                         cc =  s-c;
 245                 } else if (signbit(s) != signbit(c)) {
 246                         cc = s - c;
 247                         ss = cos(x+x)/cc;
 248                 } else {
 249                         ss = -s-c;
 250                         cc = cos(x+x)/ss;
 251                 }

 252         /*
 253          * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
 254          * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
 255          */
 256                 if (x > 1.0e91)
 257                     d =  (invsqrtpi*ss)/sqrt(x);
 258                 else
 259                         d = invsqrtpi*(pone(x)*ss+qone(x)*cc)/sqrt(x);
 260 
 261                 if (x > X_TLOSS)
 262                         return (_SVID_libm_err(x, d, 37));
 263                 else
 264                         return (d);
 265         }
 266                 if (x <= tiny) {
 267                         return (-tpi/x);
 268                 }
 269         z = x*x;


 270         if (x < 1.28) {
 271             u = u0[3]; v = v0[3]+z*v0[4];


 272             for (i = 2; i >= 0; i--) {
 273                 u = u*z + u0[i];
 274                 v = v*z + v0[i];
 275             }
 276         } else {
 277             for (u = u1[11], i = 10; i >= 0; i--) u = u*z+u1[i];
 278             v = v1[0]+z*(v1[1]+z*(v1[2]+z*(v1[3]+z*v1[4])));


 279         }
 280         return (x*(u/v) + tpi*(j1(x)*log(x)-one/x));

 281 }
 282 
 283 static const GENERIC pr0[6] = {
 284         -.4435757816794127857114720794e7,
 285         -.9942246505077641195658377899e7,
 286         -.6603373248364939109255245434e7,
 287         -.1523529351181137383255105722e7,
 288         -.1098240554345934672737413139e6,
 289         -.1611616644324610116477412898e4,
 290 };

 291 static const GENERIC ps0[6] = {
 292         -.4435757816794127856828016962e7,
 293         -.9934124389934585658967556309e7,
 294         -.6585339479723087072826915069e7,
 295         -.1511809506634160881644546358e7,
 296         -.1072638599110382011903063867e6,
 297         -.1455009440190496182453565068e4,
 298 };
 299 static const GENERIC huge    = 1.0e10;
 300 

 301 static GENERIC
 302 pone(GENERIC x) {

 303         GENERIC s, r, t, z;
 304         int i;

 305                 /* assume x > 8 */
 306         if (x > huge)
 307                 return (one);
 308 
 309         t = 8.0/x; z = t*t;
 310         r = pr0[5]; s = ps0[5]+z;



 311         for (i = 4; i >= 0; i--) {
 312                 r = z*r + pr0[i];
 313                 s = z*s + ps0[i];
 314         }
 315         return (r/s);
 316 }
 317 


 318 
 319 static const GENERIC qr0[6] = {
 320         0.3322091340985722351859704442e5,
 321         0.8514516067533570196555001171e5,
 322         0.6617883658127083517939992166e5,
 323         0.1849426287322386679652009819e5,
 324         0.1706375429020768002061283546e4,
 325         0.3526513384663603218592175580e2,
 326 };

 327 static const GENERIC qs0[6] = {
 328         0.7087128194102874357377502472e6,
 329         0.1819458042243997298924553839e7,
 330         0.1419460669603720892855755253e7,
 331         0.4002944358226697511708610813e6,
 332         0.3789022974577220264142952256e5,
 333         0.8638367769604990967475517183e3,
 334 };
 335 
 336 static GENERIC
 337 qone(GENERIC x) {

 338         GENERIC s, r, t, z;
 339         int i;

 340         if (x > huge)
 341                 return (0.375/x);
 342 
 343         t = 8.0/x; z = t*t;

 344                 /* assume x > 8 */
 345         r = qr0[5]; s = qs0[5]+z;


 346         for (i = 4; i >= 0; i--) {
 347                 r = z*r + qr0[i];
 348                 s = z*s + qs0[i];
 349         }
 350         return (t*(r/s));

 351 }


   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: j1(x),y1(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 __j1 = j1
  41 #pragma weak __y1 = y1
  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-20,
  53         one = 1.0,
  54         invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
  55         tpi = 0.636619772367581343075535053490057448;
  56 
  57 static GENERIC pone(GENERIC);
  58 static GENERIC qone(GENERIC);
  59 
  60 static const GENERIC r0[4] = {
  61         -6.250000000000002203053200981413218949548e-0002,
  62         1.600998455640072901321605101981501263762e-0003,
  63         -1.963888815948313758552511884390162864930e-0005,
  64         8.263917341093549759781339713418201620998e-0008,
  65 };
  66 
  67 static const GENERIC s0[7] = {
  68         1.0e0,
  69         1.605069137643004242395356851797873766927e-0002,
  70         1.149454623251299996428500249509098499383e-0004,
  71         3.849701673735260970379681807910852327825e-0007,
  72 };
  73 
  74 static const GENERIC r1[12] = {
  75         4.999999999999999995517408894340485471724e-0001,
  76         -6.003825028120475684835384519945468075423e-0002,
  77         2.301719899263321828388344461995355419832e-0003,
  78         -4.208494869238892934859525221654040304068e-0005,
  79         4.377745135188837783031540029700282443388e-0007,
  80         -2.854106755678624335145364226735677754179e-0009,
  81         1.234002865443952024332943901323798413689e-0011,
  82         -3.645498437039791058951273508838177134310e-0014,
  83         7.404320596071797459925377103787837414422e-0017,
  84         -1.009457448277522275262808398517024439084e-0019,
  85         8.520158355824819796968771418801019930585e-0023,
  86         -3.458159926081163274483854614601091361424e-0026,
  87 };
  88 
  89 static const GENERIC s1[5] = {
  90         1.0e0,
  91         4.923499437590484879081138588998986303306e-0003,
  92         1.054389489212184156499666953501976688452e-0005,
  93         1.180768373106166527048240364872043816050e-0008,
  94         5.942665743476099355323245707680648588540e-0012,
  95 };
  96 
  97 GENERIC
  98 j1(GENERIC x)
  99 {
 100         GENERIC z, d, s, c, ss, cc, r;
 101         int i, sgn;
 102 
 103         if (!finite(x))
 104                 return (one / x);
 105 
 106         sgn = signbit(x);
 107         x = fabs(x);
 108 
 109         if (x > 8.00) {
 110                 s = sin(x);
 111                 c = cos(x);
 112 
 113                 /* BEGIN CSTYLED */
 114                 /*
 115                  * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
 116                  * where x0 = x-3pi/4
 117                  *      Better formula:
 118                  *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 119                  *                      =  1/sqrt(2) * (sin(x) - cos(x))
 120                  *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 121                  *                      = -1/sqrt(2) * (cos(x) + sin(x))
 122                  * To avoid cancellation, use
 123                  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 124                  * to compute the worse one.
 125                  */
 126                 /* END CSTYLED */
 127                 if (x > 8.9e307) {   /* x+x may overflow */
 128                         ss = -s - c;
 129                         cc = s - c;
 130                 } else if (signbit(s) != signbit(c)) {
 131                         cc = s - c;
 132                         ss = cos(x + x) / cc;
 133                 } else {
 134                         ss = -s - c;
 135                         cc = cos(x + x) / ss;
 136                 }
 137 
 138                 /*
 139                  * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
 140                  * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
 141                  */
 142                 if (x > 1.0e40)
 143                         d = (invsqrtpi * cc) / sqrt(x);
 144                 else
 145                         d = invsqrtpi * (pone(x) * cc - qone(x) * ss) / sqrt(x);
 146 
 147                 if (x > X_TLOSS) {
 148                         if (sgn != 0) {
 149                                 d = -d;
 150                                 x = -x;
 151                         }
 152 
 153                         return (_SVID_libm_err(x, d, 36));
 154                 } else if (sgn == 0) {

 155                         return (d);
 156                 } else {
 157                         return (-d);
 158                 }
 159         }
 160 
 161         if (x <= small) {
 162                 if (x <= tiny)
 163                         d = 0.5 * x;
 164                 else
 165                         d = x * (0.5 - x * x * 0.125);
 166 
 167                 if (sgn == 0)
 168                         return (d);
 169                 else
 170                         return (-d);
 171         }
 172 
 173         z = x * x;
 174 
 175         if (x < 1.28) {
 176                 r = r0[3];
 177                 s = s0[3];
 178 
 179                 for (i = 2; i >= 0; i--) {
 180                         r = r * z + r0[i];
 181                         s = s * z + s0[i];
 182                 }
 183 
 184                 d = x * 0.5 + x * (z * (r / s));
 185         } else {
 186                 r = r1[11];
 187 
 188                 for (i = 10; i >= 0; i--)
 189                         r = r * z + r1[i];
 190 
 191                 s = s1[0] + z * (s1[1] + z * (s1[2] + z * (s1[3] + z * s1[4])));
 192                 d = x * (r / s);
 193         }
 194 
 195         if (sgn == 0)
 196                 return (d);
 197         else
 198                 return (-d);
 199 }
 200 
 201 static const GENERIC u0[4] = {
 202         -1.960570906462389461018983259589655961560e-0001,
 203         4.931824118350661953459180060007970291139e-0002,
 204         -1.626975871565393656845930125424683008677e-0003,
 205         1.359657517926394132692884168082224258360e-0005,
 206 };
 207 
 208 static const GENERIC v0[5] = {
 209         1.0e0,
 210         2.565807214838390835108224713630901653793e-0002,
 211         3.374175208978404268650522752520906231508e-0004,
 212         2.840368571306070719539936935220728843177e-0006,
 213         1.396387402048998277638900944415752207592e-0008,
 214 };
 215 
 216 static const GENERIC u1[12] = {
 217         -1.960570906462389473336339614647555351626e-0001,
 218         5.336268030335074494231369159933012844735e-0002,
 219         -2.684137504382748094149184541866332033280e-0003,
 220         5.737671618979185736981543498580051903060e-0005,
 221         -6.642696350686335339171171785557663224892e-0007,
 222         4.692417922568160354012347591960362101664e-0009,
 223         -2.161728635907789319335231338621412258355e-0011,
 224         6.727353419738316107197644431844194668702e-0014,
 225         -1.427502986803861372125234355906790573422e-0016,
 226         2.020392498726806769468143219616642940371e-0019,
 227         -1.761371948595104156753045457888272716340e-0022,
 228         7.352828391941157905175042420249225115816e-0026,
 229 };
 230 
 231 static const GENERIC v1[5] = {
 232         1.0e0,
 233         5.029187436727947764916247076102283399442e-0003,
 234         1.102693095808242775074856548927801750627e-0005,
 235         1.268035774543174837829534603830227216291e-0008,
 236         6.579416271766610825192542295821308730206e-0012,
 237 };
 238 

 239 GENERIC
 240 y1(GENERIC x)
 241 {
 242         GENERIC z, d, s, c, ss, cc, u, v;
 243         int i;
 244 
 245         if (isnan(x))
 246                 return (x * x);         /* + -> * for Cheetah */
 247 
 248         if (x <= zero) {
 249                 if (x == zero)
 250                         /* return -one/zero;  */
 251                         return (_SVID_libm_err(x, x, 10));
 252                 else
 253                         /* return zero/zero; */
 254                         return (_SVID_libm_err(x, x, 11));
 255         }
 256 
 257         if (x > 8.0) {
 258                 if (!finite(x))
 259                         return (zero);
 260 
 261                 s = sin(x);
 262                 c = cos(x);
 263 
 264                 /* BEGIN CSTYLED */
 265                 /*
 266                  * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
 267                  * where x0 = x-3pi/4
 268                  *      Better formula:
 269                  *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 270                  *                      =  1/sqrt(2) * (sin(x) - cos(x))
 271                  *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 272                  *                      = -1/sqrt(2) * (cos(x) + sin(x))
 273                  * To avoid cancellation, use
 274                  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 275                  * to compute the worse one.
 276                  */
 277                 /* END CSTYLED */
 278                 if (x > 8.9e307) {   /* x+x may overflow */
 279                         ss = -s - c;
 280                         cc = s - c;
 281                 } else if (signbit(s) != signbit(c)) {
 282                         cc = s - c;
 283                         ss = cos(x + x) / cc;
 284                 } else {
 285                         ss = -s - c;
 286                         cc = cos(x + x) / ss;
 287                 }
 288 
 289                 /*
 290                  * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
 291                  * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
 292                  */
 293                 if (x > 1.0e91)
 294                         d = (invsqrtpi * ss) / sqrt(x);
 295                 else
 296                         d = invsqrtpi * (pone(x) * ss + qone(x) * cc) / sqrt(x);
 297 
 298                 if (x > X_TLOSS)
 299                         return (_SVID_libm_err(x, d, 37));
 300                 else
 301                         return (d);
 302         }
 303 
 304         if (x <= tiny)
 305                 return (-tpi / x);
 306 
 307         z = x * x;
 308 
 309         if (x < 1.28) {
 310                 u = u0[3];
 311                 v = v0[3] + z * v0[4];
 312 
 313                 for (i = 2; i >= 0; i--) {
 314                         u = u * z + u0[i];
 315                         v = v * z + v0[i];
 316                 }
 317         } else {
 318                 for (u = u1[11], i = 10; i >= 0; i--)
 319                         u = u * z + u1[i];
 320 
 321                 v = v1[0] + z * (v1[1] + z * (v1[2] + z * (v1[3] + z * v1[4])));
 322         }
 323 
 324         return (x * (u / v) + tpi * (j1(x) * log(x) - one / x));
 325 }
 326 
 327 static const GENERIC pr0[6] = {
 328         -.4435757816794127857114720794e7, -.9942246505077641195658377899e7,
 329         -.6603373248364939109255245434e7, -.1523529351181137383255105722e7,
 330         -.1098240554345934672737413139e6, -.1611616644324610116477412898e4,



 331 };
 332 
 333 static const GENERIC ps0[6] = {
 334         -.4435757816794127856828016962e7, -.9934124389934585658967556309e7,
 335         -.6585339479723087072826915069e7, -.1511809506634160881644546358e7,
 336         -.1072638599110382011903063867e6, -.1455009440190496182453565068e4,



 337 };

 338 
 339 static const GENERIC huge = 1.0e10;
 340 static GENERIC
 341 pone(GENERIC x)
 342 {
 343         GENERIC s, r, t, z;
 344         int i;
 345 
 346         /* assume x > 8 */
 347         if (x > huge)
 348                 return (one);
 349 
 350         t = 8.0 / x;
 351         z = t * t;
 352         r = pr0[5];
 353         s = ps0[5] + z;
 354 
 355         for (i = 4; i >= 0; i--) {
 356                 r = z * r + pr0[i];
 357                 s = z * s + ps0[i];
 358         }


 359 
 360         return (r / s);
 361 }
 362 
 363 static const GENERIC qr0[6] = {
 364         0.3322091340985722351859704442e5, 0.8514516067533570196555001171e5,
 365         0.6617883658127083517939992166e5, 0.1849426287322386679652009819e5,
 366         0.1706375429020768002061283546e4, 0.3526513384663603218592175580e2,



 367 };
 368 
 369 static const GENERIC qs0[6] = {
 370         0.7087128194102874357377502472e6, 0.1819458042243997298924553839e7,
 371         0.1419460669603720892855755253e7, 0.4002944358226697511708610813e6,
 372         0.3789022974577220264142952256e5, 0.8638367769604990967475517183e3,



 373 };
 374 
 375 static GENERIC
 376 qone(GENERIC x)
 377 {
 378         GENERIC s, r, t, z;
 379         int i;
 380 
 381         if (x > huge)
 382                 return (0.375 / x);
 383 
 384         t = 8.0 / x;
 385         z = t * t;
 386         /* assume x > 8 */
 387         r = qr0[5];
 388         s = qs0[5] + z;
 389 
 390         for (i = 4; i >= 0; i--) {
 391                 r = z * r + qr0[i];
 392                 s = z * s + qs0[i];
 393         }
 394 
 395         return (t * (r / s));
 396 }