#include #include #include //0.28BE60DB9391054A7F09D5F47D4D377036D8A5664F10E4107F9458EB //0.517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D6 typedef union _f32 { float f; uint32_t i; } _f32; /* 190 bits of 2/pi for Payne-Hanek style argument reduction. */ static const unsigned int two_over_pi_f [] = { 0x00000000, 0x28be60db, 0x9391054a, 0x7f09d5f4, 0x7d4d3770, 0x36d8a566, 0x4f10e410 }; float trig_red_slowpath_f (float a, int *_q) { unsigned long int p, ia; unsigned int hi, mid, lo, tmp, i; int e, q; _f32 tmp2; float r; a = fabsf(a); if (a <= 0.785398f || a != a) { *_q = 0; return a; } tmp2.f = a; e = (int)(tmp2.i>>23)-127; tmp2.i &= 0x007fffff; tmp2.i |= 0x3f000000; ia = (unsigned int)(tmp2.f * 0x1p+32f); /* extract 96 relevant bits of 1/pi based on magnitude of argument */ i = (e < 0) ? 0 : (unsigned int)e >> 5; e = (unsigned int)e & 31; hi = two_over_pi_f [i+0]; mid = two_over_pi_f [i+1]; lo = two_over_pi_f [i+2]; tmp = two_over_pi_f [i+3]; if (e) { hi = (hi << e) | (mid >> (32 - e)); mid = (mid << e) | (lo >> (32 - e)); lo = (lo << e) | (tmp >> (32 - e)); } /* compute product x * 1/pi in 2.63 fixed-point format */ p = ia * lo; p = ia * mid + (p >> 32); p = ((ia * hi) << 32) + p; /* round quotient to nearest */ q = (int)(p >> 63); // integral portion = quadrant<1:0> p = p & 0x7fffffffffffffffULL; // fraction if (p & 0x4000000000000000ULL) { // fraction >= 0.5 p = p - 0x4000000000000000ULL; // fraction - 1.0 q = q + 1; } /* compute remainder of x * (1/pi) */ r = (float)p; r = r * 0x1.921fb54442d18p-61f; // 1.5707963267948966 * 0x1.0p-61 *_q = q; return r; } #ifdef TEST #include static const float inv_pi = 0.3183098733; /* 0x3ea2f984 */ /* * pio2_1: first 17 bit of pi/2 * pio2_1t: pi/2 - pio2_1 * pio2_2: second 17 bit of pi/2 * pio2_2t: pi/2 - (pio2_1+pio2_2) * pio2_3: third 17 bit of pi/2 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ static const float pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */ pio2_2 = 1.0804273188e-05, /* 0x37354400 */ pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */ pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */ pio2_3t = 6.1232342629e-17; /* 0x248d3132 */ float __ieee754_rem_pio2f(float x, int *q) { float fn,w,r,t = fabsf(x); unsigned int n = (unsigned int)(t*inv_pi); fn = 2*(float)n; r = t-fn*pio2_1; w = fn*pio2_2; r = r-w; t = r; /* will cover all possible cases */ w = fn*pio2_3; r = t-w; w = fn*pio2_3t-((t-r)-w); if (q) *q = n; return (r-w); } unsigned int __ieee754_npif(float x) { float f0, z = fabsf(x); unsigned int n0, n1; n0 = z*(float)M_1_PI; f0 = (z - (float)n0*M_PI); n1 = f0*(float)M_1_PI; return n0+n1; } /* * dpio2_1: first 33 bit of pi * dpio2_1t: pi - pio2_1 * dpio2_2: second 33 bit of pi * dpio2_2t: pi - (pio2_1+pio2_2) * dpio2_3: third 33 bit of pi * dpio2_3t: pi - (pio2_1+pio2_2+pio2_3) */ static const double dpio2_1 = 3.14159265346825122833e+00, /* 0x400921FB, 0x54400000 */ dpio2_1t = 1.21542010130123844986e-10, /* 0x3DE0B461, 0x1A626331 */ dpio2_2 = 1.21542010126079319532e-10, /* 0x3DE0B461, 0x1A600000 */ dpio2_2t = 4.04453249759190126308e-21, /* 0x3BB3198A, 0x2E037073 */ dpio2_3 = 4.04453249742233291159e-21, /* 0x3BB3198A, 0x2E000000 */ dpio2_3t = 1.69568553207377991399e-31; /* 0x398B839A, 0x252049C1 */ double __ieee754_rem_pio2(double x, int *q) { /* 3rd iteration need, 151 bits acc */ double fn,w,r,t = fabs(x); unsigned int n = (uint32_t) (t*M_1_PI); fn = (double)n; r = (t-(fn*dpio2_1))-(fn*dpio2_2); /* 1st round good to 85 bit */ t = r; /* will cover all possible cases */ w = fn*dpio2_3; r = t-w; w = fn*dpio2_3t-((t-r)-w); //y[1] = (r-y[0])-w; if (q) *q = n; return (r-w); } #include static const float C1 = -5.0000000000e-01, C2 = 4.1666667908e-02, /* 0x3d2aaaab */ C3 = -1.3888889225e-03, /* 0xbab60b61 */ C4 = 2.4801587642e-05, /* 0x37d00d01 */ C5 = -2.7557314297e-07, /* 0xb493f27c */ C6 = 2.0875723372e-09, /* 0x310f74f6 */ C7 = -1.1359647598e-11; /* 0xad47d74e */ // Differs from libc cosf on [0, pi/2] at only the point 0x3fc460f8 (~1.5342093) by 0.0000001005828f static inline float k_cosf(float x) { float z,r; z = x*x; r = z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))); return 1.0f + (z*r + C1*z); } static const float S1 = -0x1.555556p-03f, S2 = 0x1.111111p-07f, S3 = -0x1.a01906p-13f, S4 = 0x1.7196aap-19f, S5 = -0x1.9db448p-26f; /* Approximate sine on [-PI/4,+PI/4]. Maximum ulp error = 0.64721 * Returns -0.0f for an argument of -0.0f * Polynomial approximation based on unpublished work by T. Myklebust * Note: S3 & S4 replaced with versions from a Remez Exchange minimax polynomial computation, * and S5 added. S1 & S2 remain the same. Now this has error of at most 1ulp on [-PI/2,PI/2]. * * - bofh */ static inline float k_sinf(float x) { float z,r; z = x*x; r = z*(S1+z*(S2+z*(S3+z*(S4+z*S5)))); return x+x*r; } static const float _pio2_1 = 0x1.921f00p+00, /* 0x3fc90f80 */ _pio2_1t = 0x1.6a8886p-17, /* 0x37354443 */ _pio2_2 = 0x1.6a8800p-17, /* 0x37354400 */ _pio2_2t = 0x1.0b4610p-34; /* 0x2e85a308 */ float fast_cosf(float x) { float y = x, z; if (fabsf(y) >= 0x1.921f00p+00f) { y = (y - _pio2_1) - _pio2_1t; /* near pi/2, use 24+24+24 bit pi */ z = -k_sinf(y); } else { z = k_cosf(y); } return z; } static unsigned int c = 0xdeadbeef; void xorshift32_seed(unsigned int s) { c = !!s ? s : time(0); } unsigned int xorshift32_rand(void) { unsigned int x = c; x ^= x >> 13; // a x ^= x << 5; // b x ^= x >> 19; // c x = x * 747796405U; c = x; return x; } int main(void) { unsigned int i, q0, q1; xorshift32_seed(0); for (i = 4096; i < 1048576; i+=16) { unsigned int r = (xorshift32_rand() >> 8); float r0 = (float)r * 0x1p-24f * (float)M_PI; double f0 = (double)(i<<4)*M_PI + (double)r0; float f = (float)(i<<4)*(float)M_PI + r0; float g1 = trig_red_slowpath_f(f, &q0); float g2 = __ieee754_rem_pio2f(f, &q1); double f1 = (double)(i<<4)*M_PI + (double)g1; float g3 = (float)__ieee754_rem_pio2((double)f, 0); float s = (q0 & 1) ? -1.0f : 1.0f; q0 = __ieee754_npif(f); //float f = (float)i * 0.125f * 4096.0f * (float)M_SQRT2; //printf("%#.7f %#.7e %#.7f %#.7f\n", cosf((float)M_PI * f), k_cosf(g1)-fast_cosf(g1), cosf((float)M_PI*f)-s*fast_cosf(g1), cosf((float)M_PI*f)-s*fast_cosf(g2)); //printf("%#.7f %#.7f %#.1f %#.7f %#.7f / %#.7f %#.7f %#.7f\n", f, r0, s, g1, g2, cosf(f)-s*fast_cosf(g1), cosf(f)-s*cosf(g1), cosf(f)-s*cosf(g2)); printf("%#.7f(%#.7f) %#.7f %#.7f %#.1f %#.7f %#.7f %#.7f / %#.7f %#.7f %#.7f %#.7f\n", f, f0, f1, r0, s, g1, g2, g3, r0-g1, r0-g2, r0-g3, g1-g3); } return 0; } #endif