| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.c */ |
| 2 /* | 2 /* |
| 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
| 4 */ | 4 */ |
| 5 /* | 5 /* |
| 6 * ==================================================== | 6 * ==================================================== |
| 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 8 * | 8 * |
| 9 * Developed at SunPro, a Sun Microsystems, Inc. business. | 9 * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 10 * Permission to use, copy, modify, and distribute this | 10 * Permission to use, copy, modify, and distribute this |
| 11 * software is freely granted, provided that this notice | 11 * software is freely granted, provided that this notice |
| 12 * is preserved. | 12 * is preserved. |
| 13 * ==================================================== | 13 * ==================================================== |
| 14 */ | 14 */ |
| 15 | 15 |
| 16 #define _GNU_SOURCE | 16 #define _GNU_SOURCE |
| 17 #include "libm.h" | 17 #include "libm.h" |
| 18 | 18 |
| 19 static float pzerof(float), qzerof(float); | 19 static float pzerof(float), qzerof(float); |
| 20 | 20 |
| 21 static const float | 21 static const float invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ |
| 22 invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ | 22 tpi = 6.3661974669e-01; /* 0x3f22f983 */ |
| 23 tpi = 6.3661974669e-01; /* 0x3f22f983 */ | |
| 24 | 23 |
| 25 static float common(uint32_t ix, float x, int y0) | 24 static float common(uint32_t ix, float x, int y0) { |
| 26 { | 25 float z, s, c, ss, cc; |
| 27 » float z,s,c,ss,cc; | 26 /* |
| 28 » /* | 27 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) |
| 29 » * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) | 28 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) |
| 30 » * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) | 29 */ |
| 31 » */ | 30 s = sinf(x); |
| 32 » s = sinf(x); | 31 c = cosf(x); |
| 33 » c = cosf(x); | 32 if (y0) |
| 34 » if (y0) | 33 c = -c; |
| 35 » » c = -c; | 34 cc = s + c; |
| 36 » cc = s+c; | 35 if (ix < 0x7f000000) { |
| 37 » if (ix < 0x7f000000) { | 36 ss = s - c; |
| 38 » » ss = s-c; | 37 z = -cosf(2 * x); |
| 39 » » z = -cosf(2*x); | 38 if (s * c < 0) |
| 40 » » if (s*c < 0) | 39 cc = z / ss; |
| 41 » » » cc = z/ss; | 40 else |
| 42 » » else | 41 ss = z / cc; |
| 43 » » » ss = z/cc; | 42 if (ix < 0x58800000) { |
| 44 » » if (ix < 0x58800000) { | 43 if (y0) |
| 45 » » » if (y0) | 44 ss = -ss; |
| 46 » » » » ss = -ss; | 45 cc = pzerof(x) * cc - qzerof(x) * ss; |
| 47 » » » cc = pzerof(x)*cc-qzerof(x)*ss; | 46 } |
| 48 » » } | 47 } |
| 49 » } | 48 return invsqrtpi * cc / sqrtf(x); |
| 50 » return invsqrtpi*cc/sqrtf(x); | |
| 51 } | 49 } |
| 52 | 50 |
| 53 /* R0/S0 on [0, 2.00] */ | 51 /* R0/S0 on [0, 2.00] */ |
| 54 static const float | 52 static const float R02 = 1.5625000000e-02, /* 0x3c800000 */ |
| 55 R02 = 1.5625000000e-02, /* 0x3c800000 */ | 53 R03 = -1.8997929874e-04, /* 0xb947352e */ |
| 56 R03 = -1.8997929874e-04, /* 0xb947352e */ | 54 R04 = 1.8295404516e-06, /* 0x35f58e88 */ |
| 57 R04 = 1.8295404516e-06, /* 0x35f58e88 */ | 55 R05 = -4.6183270541e-09, /* 0xb19eaf3c */ |
| 58 R05 = -4.6183270541e-09, /* 0xb19eaf3c */ | 56 S01 = 1.5619102865e-02, /* 0x3c7fe744 */ |
| 59 S01 = 1.5619102865e-02, /* 0x3c7fe744 */ | 57 S02 = 1.1692678527e-04, /* 0x38f53697 */ |
| 60 S02 = 1.1692678527e-04, /* 0x38f53697 */ | 58 S03 = 5.1354652442e-07, /* 0x3509daa6 */ |
| 61 S03 = 5.1354652442e-07, /* 0x3509daa6 */ | 59 S04 = 1.1661400734e-09; /* 0x30a045e8 */ |
| 62 S04 = 1.1661400734e-09; /* 0x30a045e8 */ | |
| 63 | 60 |
| 64 float j0f(float x) | 61 float j0f(float x) { |
| 65 { | 62 float z, r, s; |
| 66 » float z,r,s; | 63 uint32_t ix; |
| 67 » uint32_t ix; | |
| 68 | 64 |
| 69 » GET_FLOAT_WORD(ix, x); | 65 GET_FLOAT_WORD(ix, x); |
| 70 » ix &= 0x7fffffff; | 66 ix &= 0x7fffffff; |
| 71 » if (ix >= 0x7f800000) | 67 if (ix >= 0x7f800000) |
| 72 » » return 1/(x*x); | 68 return 1 / (x * x); |
| 73 » x = fabsf(x); | 69 x = fabsf(x); |
| 74 | 70 |
| 75 » if (ix >= 0x40000000) { /* |x| >= 2 */ | 71 if (ix >= 0x40000000) { /* |x| >= 2 */ |
| 76 » » /* large ulp error near zeros */ | 72 /* large ulp error near zeros */ |
| 77 » » return common(ix, x, 0); | 73 return common(ix, x, 0); |
| 78 » } | 74 } |
| 79 » if (ix >= 0x3a000000) { /* |x| >= 2**-11 */ | 75 if (ix >= 0x3a000000) { /* |x| >= 2**-11 */ |
| 80 » » /* up to 4ulp error near 2 */ | 76 /* up to 4ulp error near 2 */ |
| 81 » » z = x*x; | 77 z = x * x; |
| 82 » » r = z*(R02+z*(R03+z*(R04+z*R05))); | 78 r = z * (R02 + z * (R03 + z * (R04 + z * R05))); |
| 83 » » s = 1+z*(S01+z*(S02+z*(S03+z*S04))); | 79 s = 1 + z * (S01 + z * (S02 + z * (S03 + z * S04))); |
| 84 » » return (1+x/2)*(1-x/2) + z*(r/s); | 80 return (1 + x / 2) * (1 - x / 2) + z * (r / s); |
| 85 » } | 81 } |
| 86 » if (ix >= 0x21800000) /* |x| >= 2**-60 */ | 82 if (ix >= 0x21800000) /* |x| >= 2**-60 */ |
| 87 » » x = 0.25f*x*x; | 83 x = 0.25f * x * x; |
| 88 » return 1 - x; | 84 return 1 - x; |
| 89 } | 85 } |
| 90 | 86 |
| 91 static const float | 87 static const float u00 = -7.3804296553e-02, /* 0xbd9726b5 */ |
| 92 u00 = -7.3804296553e-02, /* 0xbd9726b5 */ | 88 u01 = 1.7666645348e-01, /* 0x3e34e80d */ |
| 93 u01 = 1.7666645348e-01, /* 0x3e34e80d */ | 89 u02 = -1.3818567619e-02, /* 0xbc626746 */ |
| 94 u02 = -1.3818567619e-02, /* 0xbc626746 */ | 90 u03 = 3.4745343146e-04, /* 0x39b62a69 */ |
| 95 u03 = 3.4745343146e-04, /* 0x39b62a69 */ | 91 u04 = -3.8140706238e-06, /* 0xb67ff53c */ |
| 96 u04 = -3.8140706238e-06, /* 0xb67ff53c */ | 92 u05 = 1.9559013964e-08, /* 0x32a802ba */ |
| 97 u05 = 1.9559013964e-08, /* 0x32a802ba */ | 93 u06 = -3.9820518410e-11, /* 0xae2f21eb */ |
| 98 u06 = -3.9820518410e-11, /* 0xae2f21eb */ | 94 v01 = 1.2730483897e-02, /* 0x3c509385 */ |
| 99 v01 = 1.2730483897e-02, /* 0x3c509385 */ | 95 v02 = 7.6006865129e-05, /* 0x389f65e0 */ |
| 100 v02 = 7.6006865129e-05, /* 0x389f65e0 */ | 96 v03 = 2.5915085189e-07, /* 0x348b216c */ |
| 101 v03 = 2.5915085189e-07, /* 0x348b216c */ | 97 v04 = 4.4111031494e-10; /* 0x2ff280c2 */ |
| 102 v04 = 4.4111031494e-10; /* 0x2ff280c2 */ | |
| 103 | 98 |
| 104 float y0f(float x) | 99 float y0f(float x) { |
| 105 { | 100 float z, u, v; |
| 106 » float z,u,v; | 101 uint32_t ix; |
| 107 » uint32_t ix; | |
| 108 | 102 |
| 109 » GET_FLOAT_WORD(ix, x); | 103 GET_FLOAT_WORD(ix, x); |
| 110 » if ((ix & 0x7fffffff) == 0) | 104 if ((ix & 0x7fffffff) == 0) |
| 111 » » return -1/0.0f; | 105 return -1 / 0.0f; |
| 112 » if (ix>>31) | 106 if (ix >> 31) |
| 113 » » return 0/0.0f; | 107 return 0 / 0.0f; |
| 114 » if (ix >= 0x7f800000) | 108 if (ix >= 0x7f800000) |
| 115 » » return 1/x; | 109 return 1 / x; |
| 116 » if (ix >= 0x40000000) { /* |x| >= 2.0 */ | 110 if (ix >= 0x40000000) { /* |x| >= 2.0 */ |
| 117 » » /* large ulp error near zeros */ | 111 /* large ulp error near zeros */ |
| 118 » » return common(ix,x,1); | 112 return common(ix, x, 1); |
| 119 » } | 113 } |
| 120 » if (ix >= 0x39000000) { /* x >= 2**-13 */ | 114 if (ix >= 0x39000000) { /* x >= 2**-13 */ |
| 121 » » /* large ulp error at x ~= 0.89 */ | 115 /* large ulp error at x ~= 0.89 */ |
| 122 » » z = x*x; | 116 z = x * x; |
| 123 » » u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); | 117 u = u00 + |
| 124 » » v = 1+z*(v01+z*(v02+z*(v03+z*v04))); | 118 z * (u01 + z * (u02 + z * (u03 + z * (u04 + z * (u05 + z * u06))))); |
| 125 » » return u/v + tpi*(j0f(x)*logf(x)); | 119 v = 1 + z * (v01 + z * (v02 + z * (v03 + z * v04))); |
| 126 » } | 120 return u / v + tpi * (j0f(x) * logf(x)); |
| 127 » return u00 + tpi*logf(x); | 121 } |
| 122 return u00 + tpi * logf(x); |
| 128 } | 123 } |
| 129 | 124 |
| 130 /* The asymptotic expansions of pzero is | 125 /* The asymptotic expansions of pzero is |
| 131 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. | 126 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. |
| 132 * For x >= 2, We approximate pzero by | 127 * For x >= 2, We approximate pzero by |
| 133 * pzero(x) = 1 + (R/S) | 128 * pzero(x) = 1 + (R/S) |
| 134 * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 | 129 * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 |
| 135 * S = 1 + pS0*s^2 + ... + pS4*s^10 | 130 * S = 1 + pS0*s^2 + ... + pS4*s^10 |
| 136 * and | 131 * and |
| 137 * | pzero(x)-1-R/S | <= 2 ** ( -60.26) | 132 * | pzero(x)-1-R/S | <= 2 ** ( -60.26) |
| 138 */ | 133 */ |
| 139 static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ | 134 static const float pR8[6] = { |
| 140 0.0000000000e+00, /* 0x00000000 */ | 135 /* for x in [inf, 8]=1/[0,0.125] */ |
| 141 -7.0312500000e-02, /* 0xbd900000 */ | 136 0.0000000000e+00, /* 0x00000000 */ |
| 142 -8.0816707611e+00, /* 0xc1014e86 */ | 137 -7.0312500000e-02, /* 0xbd900000 */ |
| 143 -2.5706311035e+02, /* 0xc3808814 */ | 138 -8.0816707611e+00, /* 0xc1014e86 */ |
| 144 -2.4852163086e+03, /* 0xc51b5376 */ | 139 -2.5706311035e+02, /* 0xc3808814 */ |
| 145 -5.2530439453e+03, /* 0xc5a4285a */ | 140 -2.4852163086e+03, /* 0xc51b5376 */ |
| 141 -5.2530439453e+03, /* 0xc5a4285a */ |
| 146 }; | 142 }; |
| 147 static const float pS8[5] = { | 143 static const float pS8[5] = { |
| 148 1.1653436279e+02, /* 0x42e91198 */ | 144 1.1653436279e+02, /* 0x42e91198 */ |
| 149 3.8337448730e+03, /* 0x456f9beb */ | 145 3.8337448730e+03, /* 0x456f9beb */ |
| 150 4.0597855469e+04, /* 0x471e95db */ | 146 4.0597855469e+04, /* 0x471e95db */ |
| 151 1.1675296875e+05, /* 0x47e4087c */ | 147 1.1675296875e+05, /* 0x47e4087c */ |
| 152 4.7627726562e+04, /* 0x473a0bba */ | 148 4.7627726562e+04, /* 0x473a0bba */ |
| 153 }; | 149 }; |
| 154 static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ | 150 static const float pR5[6] = { |
| 155 -1.1412546255e-11, /* 0xad48c58a */ | 151 /* for x in [8,4.5454]=1/[0.125,0.22001] */ |
| 156 -7.0312492549e-02, /* 0xbd8fffff */ | 152 -1.1412546255e-11, /* 0xad48c58a */ |
| 157 -4.1596107483e+00, /* 0xc0851b88 */ | 153 -7.0312492549e-02, /* 0xbd8fffff */ |
| 158 -6.7674766541e+01, /* 0xc287597b */ | 154 -4.1596107483e+00, /* 0xc0851b88 */ |
| 159 -3.3123129272e+02, /* 0xc3a59d9b */ | 155 -6.7674766541e+01, /* 0xc287597b */ |
| 160 -3.4643338013e+02, /* 0xc3ad3779 */ | 156 -3.3123129272e+02, /* 0xc3a59d9b */ |
| 157 -3.4643338013e+02, /* 0xc3ad3779 */ |
| 161 }; | 158 }; |
| 162 static const float pS5[5] = { | 159 static const float pS5[5] = { |
| 163 6.0753936768e+01, /* 0x42730408 */ | 160 6.0753936768e+01, /* 0x42730408 */ |
| 164 1.0512523193e+03, /* 0x44836813 */ | 161 1.0512523193e+03, /* 0x44836813 */ |
| 165 5.9789707031e+03, /* 0x45bad7c4 */ | 162 5.9789707031e+03, /* 0x45bad7c4 */ |
| 166 9.6254453125e+03, /* 0x461665c8 */ | 163 9.6254453125e+03, /* 0x461665c8 */ |
| 167 2.4060581055e+03, /* 0x451660ee */ | 164 2.4060581055e+03, /* 0x451660ee */ |
| 168 }; | 165 }; |
| 169 | 166 |
| 170 static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ | 167 static const float pR3[6] = { |
| 171 -2.5470459075e-09, /* 0xb12f081b */ | 168 /* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ |
| 172 -7.0311963558e-02, /* 0xbd8fffb8 */ | 169 -2.5470459075e-09, /* 0xb12f081b */ |
| 173 -2.4090321064e+00, /* 0xc01a2d95 */ | 170 -7.0311963558e-02, /* 0xbd8fffb8 */ |
| 174 -2.1965976715e+01, /* 0xc1afba52 */ | 171 -2.4090321064e+00, /* 0xc01a2d95 */ |
| 175 -5.8079170227e+01, /* 0xc2685112 */ | 172 -2.1965976715e+01, /* 0xc1afba52 */ |
| 176 -3.1447946548e+01, /* 0xc1fb9565 */ | 173 -5.8079170227e+01, /* 0xc2685112 */ |
| 174 -3.1447946548e+01, /* 0xc1fb9565 */ |
| 177 }; | 175 }; |
| 178 static const float pS3[5] = { | 176 static const float pS3[5] = { |
| 179 3.5856033325e+01, /* 0x420f6c94 */ | 177 3.5856033325e+01, /* 0x420f6c94 */ |
| 180 3.6151397705e+02, /* 0x43b4c1ca */ | 178 3.6151397705e+02, /* 0x43b4c1ca */ |
| 181 1.1936077881e+03, /* 0x44953373 */ | 179 1.1936077881e+03, /* 0x44953373 */ |
| 182 1.1279968262e+03, /* 0x448cffe6 */ | 180 1.1279968262e+03, /* 0x448cffe6 */ |
| 183 1.7358093262e+02, /* 0x432d94b8 */ | 181 1.7358093262e+02, /* 0x432d94b8 */ |
| 184 }; | 182 }; |
| 185 | 183 |
| 186 static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ | 184 static const float pR2[6] = { |
| 187 -8.8753431271e-08, /* 0xb3be98b7 */ | 185 /* for x in [2.8570,2]=1/[0.3499,0.5] */ |
| 188 -7.0303097367e-02, /* 0xbd8ffb12 */ | 186 -8.8753431271e-08, /* 0xb3be98b7 */ |
| 189 -1.4507384300e+00, /* 0xbfb9b1cc */ | 187 -7.0303097367e-02, /* 0xbd8ffb12 */ |
| 190 -7.6356959343e+00, /* 0xc0f4579f */ | 188 -1.4507384300e+00, /* 0xbfb9b1cc */ |
| 191 -1.1193166733e+01, /* 0xc1331736 */ | 189 -7.6356959343e+00, /* 0xc0f4579f */ |
| 192 -3.2336456776e+00, /* 0xc04ef40d */ | 190 -1.1193166733e+01, /* 0xc1331736 */ |
| 191 -3.2336456776e+00, /* 0xc04ef40d */ |
| 193 }; | 192 }; |
| 194 static const float pS2[5] = { | 193 static const float pS2[5] = { |
| 195 2.2220300674e+01, /* 0x41b1c32d */ | 194 2.2220300674e+01, /* 0x41b1c32d */ |
| 196 1.3620678711e+02, /* 0x430834f0 */ | 195 1.3620678711e+02, /* 0x430834f0 */ |
| 197 2.7047027588e+02, /* 0x43873c32 */ | 196 2.7047027588e+02, /* 0x43873c32 */ |
| 198 1.5387539673e+02, /* 0x4319e01a */ | 197 1.5387539673e+02, /* 0x4319e01a */ |
| 199 1.4657617569e+01, /* 0x416a859a */ | 198 1.4657617569e+01, /* 0x416a859a */ |
| 200 }; | 199 }; |
| 201 | 200 |
| 202 static float pzerof(float x) | 201 static float pzerof(float x) { |
| 203 { | 202 const float *p, *q; |
| 204 » const float *p,*q; | 203 float_t z, r, s; |
| 205 » float_t z,r,s; | 204 uint32_t ix; |
| 206 » uint32_t ix; | |
| 207 | 205 |
| 208 » GET_FLOAT_WORD(ix, x); | 206 GET_FLOAT_WORD(ix, x); |
| 209 » ix &= 0x7fffffff; | 207 ix &= 0x7fffffff; |
| 210 » if (ix >= 0x41000000){p = pR8; q = pS8;} | 208 if (ix >= 0x41000000) { |
| 211 » else if (ix >= 0x40f71c58){p = pR5; q = pS5;} | 209 p = pR8; |
| 212 » else if (ix >= 0x4036db68){p = pR3; q = pS3;} | 210 q = pS8; |
| 213 » else /*ix >= 0x40000000*/ {p = pR2; q = pS2;} | 211 } else if (ix >= 0x40f71c58) { |
| 214 » z = 1.0f/(x*x); | 212 p = pR5; |
| 215 » r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); | 213 q = pS5; |
| 216 » s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); | 214 } else if (ix >= 0x4036db68) { |
| 217 » return 1.0f + r/s; | 215 p = pR3; |
| 216 q = pS3; |
| 217 } else /*ix >= 0x40000000*/ { |
| 218 p = pR2; |
| 219 q = pS2; |
| 220 } |
| 221 z = 1.0f / (x * x); |
| 222 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); |
| 223 s = 1.0f + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4])))); |
| 224 return 1.0f + r / s; |
| 218 } | 225 } |
| 219 | 226 |
| 220 | |
| 221 /* For x >= 8, the asymptotic expansions of qzero is | 227 /* For x >= 8, the asymptotic expansions of qzero is |
| 222 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. | 228 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. |
| 223 * We approximate pzero by | 229 * We approximate pzero by |
| 224 * qzero(x) = s*(-1.25 + (R/S)) | 230 * qzero(x) = s*(-1.25 + (R/S)) |
| 225 * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 | 231 * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 |
| 226 * S = 1 + qS0*s^2 + ... + qS5*s^12 | 232 * S = 1 + qS0*s^2 + ... + qS5*s^12 |
| 227 * and | 233 * and |
| 228 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) | 234 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) |
| 229 */ | 235 */ |
| 230 static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ | 236 static const float qR8[6] = { |
| 231 0.0000000000e+00, /* 0x00000000 */ | 237 /* for x in [inf, 8]=1/[0,0.125] */ |
| 232 7.3242187500e-02, /* 0x3d960000 */ | 238 0.0000000000e+00, /* 0x00000000 */ |
| 233 1.1768206596e+01, /* 0x413c4a93 */ | 239 7.3242187500e-02, /* 0x3d960000 */ |
| 234 5.5767340088e+02, /* 0x440b6b19 */ | 240 1.1768206596e+01, /* 0x413c4a93 */ |
| 235 8.8591972656e+03, /* 0x460a6cca */ | 241 5.5767340088e+02, /* 0x440b6b19 */ |
| 236 3.7014625000e+04, /* 0x471096a0 */ | 242 8.8591972656e+03, /* 0x460a6cca */ |
| 243 3.7014625000e+04, /* 0x471096a0 */ |
| 237 }; | 244 }; |
| 238 static const float qS8[6] = { | 245 static const float qS8[6] = { |
| 239 1.6377603149e+02, /* 0x4323c6aa */ | 246 1.6377603149e+02, /* 0x4323c6aa */ |
| 240 8.0983447266e+03, /* 0x45fd12c2 */ | 247 8.0983447266e+03, /* 0x45fd12c2 */ |
| 241 1.4253829688e+05, /* 0x480b3293 */ | 248 1.4253829688e+05, /* 0x480b3293 */ |
| 242 8.0330925000e+05, /* 0x49441ed4 */ | 249 8.0330925000e+05, /* 0x49441ed4 */ |
| 243 8.4050156250e+05, /* 0x494d3359 */ | 250 8.4050156250e+05, /* 0x494d3359 */ |
| 244 -3.4389928125e+05, /* 0xc8a7eb69 */ | 251 -3.4389928125e+05, /* 0xc8a7eb69 */ |
| 245 }; | 252 }; |
| 246 | 253 |
| 247 static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ | 254 static const float qR5[6] = { |
| 248 1.8408595828e-11, /* 0x2da1ec79 */ | 255 /* for x in [8,4.5454]=1/[0.125,0.22001] */ |
| 249 7.3242180049e-02, /* 0x3d95ffff */ | 256 1.8408595828e-11, /* 0x2da1ec79 */ |
| 250 5.8356351852e+00, /* 0x40babd86 */ | 257 7.3242180049e-02, /* 0x3d95ffff */ |
| 251 1.3511157227e+02, /* 0x43071c90 */ | 258 5.8356351852e+00, /* 0x40babd86 */ |
| 252 1.0272437744e+03, /* 0x448067cd */ | 259 1.3511157227e+02, /* 0x43071c90 */ |
| 253 1.9899779053e+03, /* 0x44f8bf4b */ | 260 1.0272437744e+03, /* 0x448067cd */ |
| 261 1.9899779053e+03, /* 0x44f8bf4b */ |
| 254 }; | 262 }; |
| 255 static const float qS5[6] = { | 263 static const float qS5[6] = { |
| 256 8.2776611328e+01, /* 0x42a58da0 */ | 264 8.2776611328e+01, /* 0x42a58da0 */ |
| 257 2.0778142090e+03, /* 0x4501dd07 */ | 265 2.0778142090e+03, /* 0x4501dd07 */ |
| 258 1.8847289062e+04, /* 0x46933e94 */ | 266 1.8847289062e+04, /* 0x46933e94 */ |
| 259 5.6751113281e+04, /* 0x475daf1d */ | 267 5.6751113281e+04, /* 0x475daf1d */ |
| 260 3.5976753906e+04, /* 0x470c88c1 */ | 268 3.5976753906e+04, /* 0x470c88c1 */ |
| 261 -5.3543427734e+03, /* 0xc5a752be */ | 269 -5.3543427734e+03, /* 0xc5a752be */ |
| 262 }; | 270 }; |
| 263 | 271 |
| 264 static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ | 272 static const float qR3[6] = { |
| 265 4.3774099900e-09, /* 0x3196681b */ | 273 /* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ |
| 266 7.3241114616e-02, /* 0x3d95ff70 */ | 274 4.3774099900e-09, /* 0x3196681b */ |
| 267 3.3442313671e+00, /* 0x405607e3 */ | 275 7.3241114616e-02, /* 0x3d95ff70 */ |
| 268 4.2621845245e+01, /* 0x422a7cc5 */ | 276 3.3442313671e+00, /* 0x405607e3 */ |
| 269 1.7080809021e+02, /* 0x432acedf */ | 277 4.2621845245e+01, /* 0x422a7cc5 */ |
| 270 1.6673394775e+02, /* 0x4326bbe4 */ | 278 1.7080809021e+02, /* 0x432acedf */ |
| 279 1.6673394775e+02, /* 0x4326bbe4 */ |
| 271 }; | 280 }; |
| 272 static const float qS3[6] = { | 281 static const float qS3[6] = { |
| 273 4.8758872986e+01, /* 0x42430916 */ | 282 4.8758872986e+01, /* 0x42430916 */ |
| 274 7.0968920898e+02, /* 0x44316c1c */ | 283 7.0968920898e+02, /* 0x44316c1c */ |
| 275 3.7041481934e+03, /* 0x4567825f */ | 284 3.7041481934e+03, /* 0x4567825f */ |
| 276 6.4604252930e+03, /* 0x45c9e367 */ | 285 6.4604252930e+03, /* 0x45c9e367 */ |
| 277 2.5163337402e+03, /* 0x451d4557 */ | 286 2.5163337402e+03, /* 0x451d4557 */ |
| 278 -1.4924745178e+02, /* 0xc3153f59 */ | 287 -1.4924745178e+02, /* 0xc3153f59 */ |
| 279 }; | 288 }; |
| 280 | 289 |
| 281 static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ | 290 static const float qR2[6] = { |
| 282 1.5044444979e-07, /* 0x342189db */ | 291 /* for x in [2.8570,2]=1/[0.3499,0.5] */ |
| 283 7.3223426938e-02, /* 0x3d95f62a */ | 292 1.5044444979e-07, /* 0x342189db */ |
| 284 1.9981917143e+00, /* 0x3fffc4bf */ | 293 7.3223426938e-02, /* 0x3d95f62a */ |
| 285 1.4495602608e+01, /* 0x4167edfd */ | 294 1.9981917143e+00, /* 0x3fffc4bf */ |
| 286 3.1666231155e+01, /* 0x41fd5471 */ | 295 1.4495602608e+01, /* 0x4167edfd */ |
| 287 1.6252708435e+01, /* 0x4182058c */ | 296 3.1666231155e+01, /* 0x41fd5471 */ |
| 297 1.6252708435e+01, /* 0x4182058c */ |
| 288 }; | 298 }; |
| 289 static const float qS2[6] = { | 299 static const float qS2[6] = { |
| 290 3.0365585327e+01, /* 0x41f2ecb8 */ | 300 3.0365585327e+01, /* 0x41f2ecb8 */ |
| 291 2.6934811401e+02, /* 0x4386ac8f */ | 301 2.6934811401e+02, /* 0x4386ac8f */ |
| 292 8.4478375244e+02, /* 0x44533229 */ | 302 8.4478375244e+02, /* 0x44533229 */ |
| 293 8.8293585205e+02, /* 0x445cbbe5 */ | 303 8.8293585205e+02, /* 0x445cbbe5 */ |
| 294 2.1266638184e+02, /* 0x4354aa98 */ | 304 2.1266638184e+02, /* 0x4354aa98 */ |
| 295 -5.3109550476e+00, /* 0xc0a9f358 */ | 305 -5.3109550476e+00, /* 0xc0a9f358 */ |
| 296 }; | 306 }; |
| 297 | 307 |
| 298 static float qzerof(float x) | 308 static float qzerof(float x) { |
| 299 { | 309 const float *p, *q; |
| 300 » const float *p,*q; | 310 float_t s, r, z; |
| 301 » float_t s,r,z; | 311 uint32_t ix; |
| 302 » uint32_t ix; | |
| 303 | 312 |
| 304 » GET_FLOAT_WORD(ix, x); | 313 GET_FLOAT_WORD(ix, x); |
| 305 » ix &= 0x7fffffff; | 314 ix &= 0x7fffffff; |
| 306 » if (ix >= 0x41000000){p = qR8; q = qS8;} | 315 if (ix >= 0x41000000) { |
| 307 » else if (ix >= 0x40f71c58){p = qR5; q = qS5;} | 316 p = qR8; |
| 308 » else if (ix >= 0x4036db68){p = qR3; q = qS3;} | 317 q = qS8; |
| 309 » else /*ix >= 0x40000000*/ {p = qR2; q = qS2;} | 318 } else if (ix >= 0x40f71c58) { |
| 310 » z = 1.0f/(x*x); | 319 p = qR5; |
| 311 » r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); | 320 q = qS5; |
| 312 » s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); | 321 } else if (ix >= 0x4036db68) { |
| 313 » return (-.125f + r/s)/x; | 322 p = qR3; |
| 323 q = qS3; |
| 324 } else /*ix >= 0x40000000*/ { |
| 325 p = qR2; |
| 326 q = qS2; |
| 327 } |
| 328 z = 1.0f / (x * x); |
| 329 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); |
| 330 s = 1.0f + |
| 331 z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5]))))); |
| 332 return (-.125f + r / s) / x; |
| 314 } | 333 } |
| OLD | NEW |