#include #include #include //0.28BE60DB9391054A7F09D5F47D4D377036D8A5664F10E4107F9458EB //0.517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D6 typedef union _f32 { float f; uint32_t i; } _f32; /* 191 bits of 1/pi for Payne-Hanek style argument reduction. */ static const unsigned int one_over_pi_f [] = { 0x00000000, 0x517CC1B7, 0x27220A94, 0xFE13ABE8, 0xFA9A6EE0, 0x6DB14ACC, 0x9E21C821 }; float trig_red_slowpath_f (float a, unsigned int *_q) { unsigned long int p0, ia; unsigned int p, hi, mid, lo, tmp, i; int e, q; double d; _f32 tmp2; float r; a = fabsf(a); if (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 2/pi based on magnitude of argument */ i = (e < 0) ? 0 : (unsigned int)e >> 5; e = (unsigned int)e & 31; hi = one_over_pi_f [i+0]; mid = one_over_pi_f [i+1]; lo = one_over_pi_f [i+2]; tmp = one_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 * 2/pi in 2.62 fixed-point format */ p0 = ia * lo; p0 = ia * mid + (p0 >> 32); p0 = ((ia * hi) << 32) + p0; /* round quotient to nearest */ q = (int)(p0 >> 63); // integral portion = quadrant<1:0> p = (p0 >> 32U); p = p & 0x7fffffffU; // fraction #if 1 if (p & 0x40000000) { // fraction >= 0.5 p = p - 0x80000000U; // fraction - 1.0 q ^= 1; } #endif /* compute remainder of x * (2/pi) */ //d = (double)((int64_t)p); d = (double)((int32_t)p); d = d * 0x1.921fb54442d18p+00; // 1.5707963267948966 r = (float)d; //r = ((3.141571f*((int64_t)p)) + (2.165358979e-05f*((int64_t)p))); //r = (float)p; //r = r * 0x1.921fb54442d18p+00f; // 1.5707963267948966 * 0x1.0p-61 r = r * 0x1.0p-30f; *_q = q; return r; } #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 double pi_1 = 0x1.921fb544p+01, /* 0x3FF921FB, 0x54400000 */ pi_2 = 0x1.0b4611a6p-33, /* 0x3DD0B461, 0x1A600000 */ pi_3 = 0x1.3198a2e0p-68, /* 0x3BA3198A, 0x2E000000 */ pi_2t = 0x1.3198a2e037073p-68, /* 0x3BA3198A, 0x2E037073 */ pi_3t = 0x1.b839a252049c1p-103; /* 0x397B839A, 0x252049C1 */ __attribute__((noinline)) unsigned int __ieee754_rem_pif(float x, float *y) { unsigned int n = 0; float z; x = fabsf(x); n = lrintf(x*inv_pi); if (x < 0x1.921fb54p+00f || x != x) { z = x; } else if (x < 0x1p+22f) { double t = (double)x; double r = (((t-n*pi_1)-n*pi_2)-n*pi_2t); z = (float)r; } else { z = trig_red_slowpath_f(x, &n); } //if (y) *y = z; return n&1; } static const float C1 = 4.1666667908e-02f, /* 0x3d2aaaab */ C2 = -1.3888889225e-03f, /* 0xbab60b61 */ C3 = 2.4801587642e-05f, /* 0x37d00d01 */ C4 = -2.7557314297e-07f, /* 0xb493f27c */ C5 = 2.0875723372e-09f, /* 0x310f74f6 */ C6 = -1.1359647598e-11f, /* 0xad47d74e */ S1 = -0x1.555556p-03f, S2 = 0x1.111111p-07f, S3 = -0x1.a01906p-13f, S4 = 0x1.7196aap-19f, S5 = -0x1.9db448p-26f, S6 = 0x1.15b3a6p-34f; // 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*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); //return 1.0f + (z*r - 0.5f*z); return 0.5f + (z*r - 0.5f*(z-1.0f)); //return 1.0 + ((double)(z*r) - (double)(0.5f*z)); } /* 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+z*S6))))); return x+x*r; } float fast_cosf(float x) { float s, z; unsigned int n = __ieee754_rem_pif(x, &z); s = (n & 1) ? -1.0f : 1.0f; return s*k_cosf(z); } static const float pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */ pio2_1t = 1.0804334124e-05; /* 0x37354443 */ void fast_sincosf(float x, float *cs, float *sn) { float c, s, sgn, z; unsigned int n = __ieee754_rem_pif(x, &z); n <<= 1; if (z > 1.0f) { z = ((z - pio2_1) - pio2_1t); n |= 1; } else if (z < -1.0f) { z = ((z + pio2_1) + pio2_1t); n ^= 3; } sgn = (n & 2) ? -1.0f : 1.0f; c = sgn*k_cosf(z); s = sgn*k_sinf(z); *cs = ((n&1)?-s:c); *sn = ((n&1)? c:s); } #ifdef TEST static __attribute__((always_inline)) inline void disable_denormals(void) { unsigned int a = 0x9f80; __asm__ __volatile__ ("ldmxcsr %0" : "+m" (a)); } int main(void) { _f32 f; unsigned int i, j, nerrs = 0; float c0, c1, s0, s1, diff, maxdiff = 0.0f, maxdiff2 = 0.0f; disable_denormals(); setvbuf(stdout, NULL, _IONBF, 0); nerrs = 0; maxdiff = 0.0f; maxdiff2 = 0.0f; for (j = 120; j < 171; j++) { f.i = (j<<23); printf("start: %.7f (j: %u)\n", f.f, j); for (i = 0; i < (1 << 23U); i++) { f.i = (j<<23)|i; fast_sincosf(f.f, &c0, &s0); c1 = cosf(f.f); s1 = sinf(f.f); if ((diff = fabsf(c0-c1)) >= 0x1.0p-23f) { nerrs++; if (diff > fabsf(maxdiff)) maxdiff = diff; if (!(nerrs & 0x3fff)) printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / %a\n", f.f, c1, f.f, c0, diff, diff); } if ((diff = fabsf(s0-s1)) > 0x1.0p-23f) { nerrs++; if (diff > fabsf(maxdiff2)) maxdiff2 = diff; if (!(nerrs & 0x3fff)) printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / %a\n", f.f, s1, f.f, s0, diff, diff); } } printf("end: %.7f\n", f.f); } printf("Total errors: %u, largest error: %.7e/%.7e\n", nerrs, maxdiff, maxdiff2 ); return 0; } #endif