lib/sdk/crt/math/cos.c

   1 /*
   2  * COPYRIGHT:        See COPYING in the top level directory
   3  * PROJECT:          ReactOS CRT
   4  * FILE:             lib/crt/math/cos.c
   5  * PURPOSE:          Generic C Implementation of cos
   6  * PROGRAMMER:       Timo Kreuzer (timo.kreuzer@reactos.org)
   7  */
   8
   9 #define PRECISION 9
  10 #define M_PI 3.141592653589793238462643
  11
  12 static double cos_off_tbl[] = {0.0, -M_PI/2., 0, -M_PI/2.};
  13 static double cos_sign_tbl[] = {1,-1,-1,1};
  14
  15 double
  16 cos(double x)
  17 {
  18     int quadrant;
  19     double x2, result;
  20
  21     /* Calculate the quadrant */
  22     quadrant = x * (2./M_PI);
  23
  24     /* Get offset inside quadrant */
  25     x = x - quadrant * (M_PI/2.);
  26
  27     /* Normalize quadrant to [0..3] */
  28     quadrant = quadrant & 0x3;
  29
  30     /* Fixup value for the generic function */
  31     x += cos_off_tbl[quadrant];
  32
  33     /* Calculate the negative of the square of x */
  34     x2 = - (x * x);
  35
  36     /* This is an unrolled taylor series using <PRECISION> iterations
  37      * Example with 4 iterations:
  38      * result = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8!
  39      * To save multiplications and to keep the precision high, it's performed
  40      * like this:
  41      * result = 1 - x^2 * (1/2! - x^2 * (1/4! - x^2 * (1/6! - x^2 * (1/8!))))
  42      */
  43
  44     /* Start with 0, compiler will optimize this away */
  45     result = 0;
  46
  47 #if (PRECISION >= 10)
  48     result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20);
  49     result *= x2;
  50 #endif
  51 #if (PRECISION >= 9)
  52     result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18);
  53     result *= x2;
  54 #endif
  55 #if (PRECISION >= 8)
  56     result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16);
  57     result *= x2;
  58 #endif
  59 #if (PRECISION >= 7)
  60     result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14);
  61     result *= x2;
  62 #endif
  63 #if (PRECISION >= 6)
  64     result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12);
  65     result *= x2;
  66 #endif
  67 #if (PRECISION >= 5)
  68     result += 1./(1.*2*3*4*5*6*7*8*9*10);
  69     result *= x2;
  70 #endif
  71     result += 1./(1.*2*3*4*5*6*7*8);
  72     result *= x2;
  73
  74     result += 1./(1.*2*3*4*5*6);
  75     result *= x2;
  76
  77     result += 1./(1.*2*3*4);
  78     result *= x2;
  79
  80     result += 1./(1.*2);
  81     result *= x2;
  82
  83     result += 1;
  84
  85     /* Apply correct sign */
  86     result *= cos_sign_tbl[quadrant];
  87
  88     return result;
  89 }