1/* origin: FreeBSD /usr/src/lib/msun/src/e_j1f.c */ 2/* 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 4 */ 5/* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16#define _GNU_SOURCE 17#include "libm.h" 18 19static float ponef(float), qonef(float); 20 21static const float 22invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ 23tpi = 6.3661974669e-01; /* 0x3f22f983 */ 24 25static float common(uint32_t ix, float x, int y1, int sign) 26{ 27 double z,s,c,ss,cc; 28 29 s = sinf(x); 30 if (y1) 31 s = -s; 32 c = cosf(x); 33 cc = s-c; 34 if (ix < 0x7f000000) { 35 ss = -s-c; 36 z = cosf(2*x); 37 if (s*c > 0) 38 cc = z/ss; 39 else 40 ss = z/cc; 41 if (ix < 0x58800000) { 42 if (y1) 43 ss = -ss; 44 cc = ponef(x)*cc-qonef(x)*ss; 45 } 46 } 47 if (sign) 48 cc = -cc; 49 return invsqrtpi*cc/sqrtf(x); 50} 51 52/* R0/S0 on [0,2] */ 53static const float 54r00 = -6.2500000000e-02, /* 0xbd800000 */ 55r01 = 1.4070566976e-03, /* 0x3ab86cfd */ 56r02 = -1.5995563444e-05, /* 0xb7862e36 */ 57r03 = 4.9672799207e-08, /* 0x335557d2 */ 58s01 = 1.9153760746e-02, /* 0x3c9ce859 */ 59s02 = 1.8594678841e-04, /* 0x3942fab6 */ 60s03 = 1.1771846857e-06, /* 0x359dffc2 */ 61s04 = 5.0463624390e-09, /* 0x31ad6446 */ 62s05 = 1.2354227016e-11; /* 0x2d59567e */ 63 64float j1f(float x) 65{ 66 float z,r,s; 67 uint32_t ix; 68 int sign; 69 70 GET_FLOAT_WORD(ix, x); 71 sign = ix>>31; 72 ix &= 0x7fffffff; 73 if (ix >= 0x7f800000) 74 return 1/(x*x); 75 if (ix >= 0x40000000) /* |x| >= 2 */ 76 return common(ix, fabsf(x), 0, sign); 77 if (ix >= 0x32000000) { /* |x| >= 2**-27 */ 78 z = x*x; 79 r = z*(r00+z*(r01+z*(r02+z*r03))); 80 s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); 81 z = 0.5f + r/s; 82 } else 83 /* raise inexact if x!=0 */ 84 z = 0.5f + x; 85 return z*x; 86} 87 88static const float U0[5] = { 89 -1.9605709612e-01, /* 0xbe48c331 */ 90 5.0443872809e-02, /* 0x3d4e9e3c */ 91 -1.9125689287e-03, /* 0xbafaaf2a */ 92 2.3525259166e-05, /* 0x37c5581c */ 93 -9.1909917899e-08, /* 0xb3c56003 */ 94}; 95static const float V0[5] = { 96 1.9916731864e-02, /* 0x3ca3286a */ 97 2.0255257550e-04, /* 0x3954644b */ 98 1.3560879779e-06, /* 0x35b602d4 */ 99 6.2274145840e-09, /* 0x31d5f8eb */ 100 1.6655924903e-11, /* 0x2d9281cf */ 101}; 102 103float y1f(float x) 104{ 105 float z,u,v; 106 uint32_t ix; 107 108 GET_FLOAT_WORD(ix, x); 109 if ((ix & 0x7fffffff) == 0) 110 return -1/0.0f; 111 if (ix>>31) 112 return 0/0.0f; 113 if (ix >= 0x7f800000) 114 return 1/x; 115 if (ix >= 0x40000000) /* |x| >= 2.0 */ 116 return common(ix,x,1,0); 117 if (ix < 0x32000000) /* x < 2**-27 */ 118 return -tpi/x; 119 z = x*x; 120 u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); 121 v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); 122 return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x); 123} 124 125/* For x >= 8, the asymptotic expansions of pone is 126 * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. 127 * We approximate pone by 128 * pone(x) = 1 + (R/S) 129 * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 130 * S = 1 + ps0*s^2 + ... + ps4*s^10 131 * and 132 * | pone(x)-1-R/S | <= 2 ** ( -60.06) 133 */ 134 135static const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 136 0.0000000000e+00, /* 0x00000000 */ 137 1.1718750000e-01, /* 0x3df00000 */ 138 1.3239480972e+01, /* 0x4153d4ea */ 139 4.1205184937e+02, /* 0x43ce06a3 */ 140 3.8747453613e+03, /* 0x45722bed */ 141 7.9144794922e+03, /* 0x45f753d6 */ 142}; 143static const float ps8[5] = { 144 1.1420736694e+02, /* 0x42e46a2c */ 145 3.6509309082e+03, /* 0x45642ee5 */ 146 3.6956207031e+04, /* 0x47105c35 */ 147 9.7602796875e+04, /* 0x47bea166 */ 148 3.0804271484e+04, /* 0x46f0a88b */ 149}; 150 151static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 152 1.3199052094e-11, /* 0x2d68333f */ 153 1.1718749255e-01, /* 0x3defffff */ 154 6.8027510643e+00, /* 0x40d9b023 */ 155 1.0830818176e+02, /* 0x42d89dca */ 156 5.1763616943e+02, /* 0x440168b7 */ 157 5.2871520996e+02, /* 0x44042dc6 */ 158}; 159static const float ps5[5] = { 160 5.9280597687e+01, /* 0x426d1f55 */ 161 9.9140142822e+02, /* 0x4477d9b1 */ 162 5.3532670898e+03, /* 0x45a74a23 */ 163 7.8446904297e+03, /* 0x45f52586 */ 164 1.5040468750e+03, /* 0x44bc0180 */ 165}; 166 167static const float pr3[6] = { 168 3.0250391081e-09, /* 0x314fe10d */ 169 1.1718686670e-01, /* 0x3defffab */ 170 3.9329774380e+00, /* 0x407bb5e7 */ 171 3.5119403839e+01, /* 0x420c7a45 */ 172 9.1055007935e+01, /* 0x42b61c2a */ 173 4.8559066772e+01, /* 0x42423c7c */ 174}; 175static const float ps3[5] = { 176 3.4791309357e+01, /* 0x420b2a4d */ 177 3.3676245117e+02, /* 0x43a86198 */ 178 1.0468714600e+03, /* 0x4482dbe3 */ 179 8.9081134033e+02, /* 0x445eb3ed */ 180 1.0378793335e+02, /* 0x42cf936c */ 181}; 182 183static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 184 1.0771083225e-07, /* 0x33e74ea8 */ 185 1.1717621982e-01, /* 0x3deffa16 */ 186 2.3685150146e+00, /* 0x401795c0 */ 187 1.2242610931e+01, /* 0x4143e1bc */ 188 1.7693971634e+01, /* 0x418d8d41 */ 189 5.0735230446e+00, /* 0x40a25a4d */ 190}; 191static const float ps2[5] = { 192 2.1436485291e+01, /* 0x41ab7dec */ 193 1.2529022980e+02, /* 0x42fa9499 */ 194 2.3227647400e+02, /* 0x436846c7 */ 195 1.1767937469e+02, /* 0x42eb5bd7 */ 196 8.3646392822e+00, /* 0x4105d590 */ 197}; 198 199static float ponef(float x) 200{ 201 const float *p,*q; 202 float_t z,r,s; 203 uint32_t ix; 204 205 GET_FLOAT_WORD(ix, x); 206 ix &= 0x7fffffff; 207 if (ix >= 0x41000000){p = pr8; q = ps8;} 208 else if (ix >= 0x40f71c58){p = pr5; q = ps5;} 209 else if (ix >= 0x4036db68){p = pr3; q = ps3;} 210 else /*ix >= 0x40000000*/ {p = pr2; q = ps2;} 211 z = 1.0f/(x*x); 212 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 213 s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 214 return 1.0f + r/s; 215} 216 217/* For x >= 8, the asymptotic expansions of qone is 218 * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. 219 * We approximate pone by 220 * qone(x) = s*(0.375 + (R/S)) 221 * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 222 * S = 1 + qs1*s^2 + ... + qs6*s^12 223 * and 224 * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) 225 */ 226 227static const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 228 0.0000000000e+00, /* 0x00000000 */ 229 -1.0253906250e-01, /* 0xbdd20000 */ 230 -1.6271753311e+01, /* 0xc1822c8d */ 231 -7.5960174561e+02, /* 0xc43de683 */ 232 -1.1849806641e+04, /* 0xc639273a */ 233 -4.8438511719e+04, /* 0xc73d3683 */ 234}; 235static const float qs8[6] = { 236 1.6139537048e+02, /* 0x43216537 */ 237 7.8253862305e+03, /* 0x45f48b17 */ 238 1.3387534375e+05, /* 0x4802bcd6 */ 239 7.1965775000e+05, /* 0x492fb29c */ 240 6.6660125000e+05, /* 0x4922be94 */ 241 -2.9449025000e+05, /* 0xc88fcb48 */ 242}; 243 244static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 245 -2.0897993405e-11, /* 0xadb7d219 */ 246 -1.0253904760e-01, /* 0xbdd1fffe */ 247 -8.0564479828e+00, /* 0xc100e736 */ 248 -1.8366960144e+02, /* 0xc337ab6b */ 249 -1.3731937256e+03, /* 0xc4aba633 */ 250 -2.6124443359e+03, /* 0xc523471c */ 251}; 252static const float qs5[6] = { 253 8.1276550293e+01, /* 0x42a28d98 */ 254 1.9917987061e+03, /* 0x44f8f98f */ 255 1.7468484375e+04, /* 0x468878f8 */ 256 4.9851425781e+04, /* 0x4742bb6d */ 257 2.7948074219e+04, /* 0x46da5826 */ 258 -4.7191835938e+03, /* 0xc5937978 */ 259}; 260 261static const float qr3[6] = { 262 -5.0783124372e-09, /* 0xb1ae7d4f */ 263 -1.0253783315e-01, /* 0xbdd1ff5b */ 264 -4.6101160049e+00, /* 0xc0938612 */ 265 -5.7847221375e+01, /* 0xc267638e */ 266 -2.2824453735e+02, /* 0xc3643e9a */ 267 -2.1921012878e+02, /* 0xc35b35cb */ 268}; 269static const float qs3[6] = { 270 4.7665153503e+01, /* 0x423ea91e */ 271 6.7386511230e+02, /* 0x4428775e */ 272 3.3801528320e+03, /* 0x45534272 */ 273 5.5477290039e+03, /* 0x45ad5dd5 */ 274 1.9031191406e+03, /* 0x44ede3d0 */ 275 -1.3520118713e+02, /* 0xc3073381 */ 276}; 277 278static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 279 -1.7838172539e-07, /* 0xb43f8932 */ 280 -1.0251704603e-01, /* 0xbdd1f475 */ 281 -2.7522056103e+00, /* 0xc0302423 */ 282 -1.9663616180e+01, /* 0xc19d4f16 */ 283 -4.2325313568e+01, /* 0xc2294d1f */ 284 -2.1371921539e+01, /* 0xc1aaf9b2 */ 285}; 286static const float qs2[6] = { 287 2.9533363342e+01, /* 0x41ec4454 */ 288 2.5298155212e+02, /* 0x437cfb47 */ 289 7.5750280762e+02, /* 0x443d602e */ 290 7.3939318848e+02, /* 0x4438d92a */ 291 1.5594900513e+02, /* 0x431bf2f2 */ 292 -4.9594988823e+00, /* 0xc09eb437 */ 293}; 294 295static float qonef(float x) 296{ 297 const float *p,*q; 298 float_t s,r,z; 299 uint32_t ix; 300 301 GET_FLOAT_WORD(ix, x); 302 ix &= 0x7fffffff; 303 if (ix >= 0x40200000){p = qr8; q = qs8;} 304 else if (ix >= 0x40f71c58){p = qr5; q = qs5;} 305 else if (ix >= 0x4036db68){p = qr3; q = qs3;} 306 else /*ix >= 0x40000000*/ {p = qr2; q = qs2;} 307 z = 1.0f/(x*x); 308 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 309 s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 310 return (.375f + r/s)/x; 311} 312