Welcome to Fractal Forums

1. x = x * 0.15915494309189533576888376337251; x = x - roundint (x);
2. Fold (x, 0.25); // abs(x+Fold) - abs(x-Fold) - x
3. ox =x; x = sqr(x); x2 = x^2 ... etc
4. x = ox* (6.283163944 - 41.3371315 x + 81.34042392 x2 -70.99090501 x3 )

(x,y){
//return sin(x);
#define INV_TWO_PI 0.15915494309189533576888376337251
#define TWO_PI 6.283185307179586476925286766559
   x=x*INV_TWO_PI+(-0.25);
   x=x-floor(x);
   x=abs(x+(-0.5))+(-0.25);
   x=x*TWO_PI;
   x2=x*x;
   //from: http://stackoverflow.com/questions/345085/how-do-trigonometric-functions-work/394512#394512
   x = x * (0.99999660 + x2*(-0.16664824 + x2*(0.00830629 +x2*-0.00018363)));
}

1244240:83 ec 0c                   sub esp, 12
1244243:f2 0f 10 44 24 10          movsd xmm0, [esp+16]
1244249:f2 0f 59 46 30             mulsd xmm0, [esi+48]
124424e:f2 0f 58 46 10             addsd xmm0, [esi+16]
1244253:f2 0f 11 44 24 10          movsd [esp+16], xmm0
1244259:dd 44 24 10                fld qword ptr [esp+16]
124425d:db 5c 24 f8                fistp dword ptr [esp-8]
1244261:f2 0f 2a 4c 24 f8          cvtsi2sd xmm1, [esp-8]
1244267:f2 0f 5c c1                subsd xmm0, xmm1
124426b:f2 0f 11 44 24 10          movsd [esp+16], xmm0
1244271:f2 0f 58 46 08             addsd xmm0, [esi+8]
1244276:0f 54 46 e0                andps xmm0, [esi-32]
124427a:f2 0f 58 46 10             addsd xmm0, [esi+16]
124427f:f2 0f 11 44 24 10          movsd [esp+16], xmm0
1244285:f2 0f 59 46 58             mulsd xmm0, [esi+88]
124428a:f2 0f 11 44 24 10          movsd [esp+16], xmm0
1244290:f2 0f 59 c0                mulsd xmm0, xmm0
1244294:f2 0f 11 04 24             movsd [esp], xmm0
1244299:f2 0f 59 46 20             mulsd xmm0, [esi+32]
124429e:f2 0f 58 46 28             addsd xmm0, [esi+40]
12442a3:f2 0f 59 04 24             mulsd xmm0, [esp]
12442a8:f2 0f 58 46 18             addsd xmm0, [esi+24]
12442ad:f2 0f 59 04 24             mulsd xmm0, [esp]
12442b2:f2 0f 58 46 48             addsd xmm0, [esi+72]
12442b7:f2 0f 59 44 24 10          mulsd xmm0, [esp+16]
12442bd:f2 0f 11 44 24 10          movsd [esp+16], xmm0
12442c3:dd 44 24 10                fld qword ptr [esp+16]
12442c7:83 c4 0c                   add esp, 12
12442ca:c3                         ret

[CONSTANTS]
Double = 0.31830988618379067153776752674503 // [edi]
Double = 3.1415926535897932384626433832795 // [edi+8] ... and so on
Double = -5.167375251323450
Double = 2.543310591226590
Double = -0.557168173477934

// gcc -std=c99 -Wall -Wextra -pedantic -O3 -march=native -ffast-math -funroll-loops -D_BSD_SOURCE -D_POSIX_C_SOURCE=200112L
// note: -ffast-math can break some numerical code, verify results carefully, use it only for translation units where it's known to be ok
#ifndef COSF9_H
#define COSF9_H 1

#include <math.h>

#define   likely(x) __builtin_expect((x),1)
#define unlikely(x) __builtin_expect((x),0)

// postcondition: 0 <= phase <= 1
static inline float cosf9_reduce(float phase) {
  float p = fabsf(phase);
  // p >= 0
  if (likely(p < (1 << 24))) {
    int q = p;
    return p - q;
  } else {
    if (unlikely(isnanf(p) || isinff(p))) {
      // return NaN
      return p - p;
    } else {
      // int could overflow, and it will be integral anyway
      return 0.0f;
    }
  }
}

// precondition: 0 <= phase <= 1
static inline float cosf9_unsafe(float phase) {
  float p = fabsf(4.0f * phase - 2.0f) - 1.0f;
  // p in -1 .. 1
  float s
    = 1.5707963267948965580e+00f * p
    - 6.4596271553942852250e-01f * p * p * p
    + 7.9685048314861006702e-02f * p * p * p * p * p
    - 4.6672571910271187789e-03f * p * p * p * p * p * p * p
    + 1.4859762069630022552e-04f * p * p * p * p * p * p * p * p * p;
  // compiler figures out optimal simd multiplications
  return s;
}

static inline float cosf9(float phase) {
  return cosf9_unsafe(cosf9_reduce(phase));
}

#endif