1 /* 2 * Copyright (c) 2006-2018, RT-Thread Development Team 3 * 4 * SPDX-License-Identifier: Apache-2.0 5 * 6 * Change Logs: 7 * Date Author Notes 8 */ 9 #include <math.h> 10 11 /* 12 * COPYRIGHT: See COPYING in the top level directory 13 * PROJECT: ReactOS CRT 14 * FILE: lib/crt/math/cos.c 15 * PURPOSE: Generic C Implementation of cos 16 * PROGRAMMER: Timo Kreuzer ([email protected]) 17 */ 18 19 #define PRECISION 9 20 21 static double cos_off_tbl[] = {0.0, -M_PI/2., 0, -M_PI/2.}; 22 static double cos_sign_tbl[] = {1,-1,-1,1}; 23 24 static double sin_off_tbl[] = {0.0, -M_PI/2., 0, -M_PI/2.}; 25 static double sin_sign_tbl[] = {1,-1,-1,1}; 26 27 double sin(double x) 28 { 29 int quadrant; 30 double x2, result; 31 32 /* Calculate the quadrant */ 33 quadrant = x * (2./M_PI); 34 35 /* Get offset inside quadrant */ 36 x = x - quadrant * (M_PI/2.); 37 38 /* Normalize quadrant to [0..3] */ 39 quadrant = (quadrant - 1) & 0x3; 40 41 /* Fixup value for the generic function */ 42 x += sin_off_tbl[quadrant]; 43 44 /* Calculate the negative of the square of x */ 45 x2 = - (x * x); 46 47 /* This is an unrolled taylor series using <PRECISION> iterations 48 * Example with 4 iterations: 49 * result = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! 50 * To save multiplications and to keep the precision high, it's performed 51 * like this: 52 * result = 1 - x^2 * (1/2! - x^2 * (1/4! - x^2 * (1/6! - x^2 * (1/8!)))) 53 */ 54 55 /* Start with 0, compiler will optimize this away */ 56 result = 0; 57 58 #if (PRECISION >= 10) 59 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20); 60 result *= x2; 61 #endif 62 #if (PRECISION >= 9) 63 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18); 64 result *= x2; 65 #endif 66 #if (PRECISION >= 8) 67 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16); 68 result *= x2; 69 #endif 70 #if (PRECISION >= 7) 71 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14); 72 result *= x2; 73 #endif 74 #if (PRECISION >= 6) 75 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12); 76 result *= x2; 77 #endif 78 #if (PRECISION >= 5) 79 result += 1./(1.*2*3*4*5*6*7*8*9*10); 80 result *= x2; 81 #endif 82 result += 1./(1.*2*3*4*5*6*7*8); 83 result *= x2; 84 85 result += 1./(1.*2*3*4*5*6); 86 result *= x2; 87 88 result += 1./(1.*2*3*4); 89 result *= x2; 90 91 result += 1./(1.*2); 92 result *= x2; 93 94 result += 1; 95 96 /* Apply correct sign */ 97 result *= sin_sign_tbl[quadrant]; 98 99 return result; 100 } 101 102 double cos(double x) 103 { 104 int quadrant; 105 double x2, result; 106 107 /* Calculate the quadrant */ 108 quadrant = x * (2./M_PI); 109 110 /* Get offset inside quadrant */ 111 x = x - quadrant * (M_PI/2.); 112 113 /* Normalize quadrant to [0..3] */ 114 quadrant = quadrant & 0x3; 115 116 /* Fixup value for the generic function */ 117 x += cos_off_tbl[quadrant]; 118 119 /* Calculate the negative of the square of x */ 120 x2 = - (x * x); 121 122 /* This is an unrolled taylor series using <PRECISION> iterations 123 * Example with 4 iterations: 124 * result = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! 125 * To save multiplications and to keep the precision high, it's performed 126 * like this: 127 * result = 1 - x^2 * (1/2! - x^2 * (1/4! - x^2 * (1/6! - x^2 * (1/8!)))) 128 */ 129 130 /* Start with 0, compiler will optimize this away */ 131 result = 0; 132 133 #if (PRECISION >= 10) 134 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20); 135 result *= x2; 136 #endif 137 #if (PRECISION >= 9) 138 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18); 139 result *= x2; 140 #endif 141 #if (PRECISION >= 8) 142 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16); 143 result *= x2; 144 #endif 145 #if (PRECISION >= 7) 146 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14); 147 result *= x2; 148 #endif 149 #if (PRECISION >= 6) 150 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12); 151 result *= x2; 152 #endif 153 #if (PRECISION >= 5) 154 result += 1./(1.*2*3*4*5*6*7*8*9*10); 155 result *= x2; 156 #endif 157 result += 1./(1.*2*3*4*5*6*7*8); 158 result *= x2; 159 160 result += 1./(1.*2*3*4*5*6); 161 result *= x2; 162 163 result += 1./(1.*2*3*4); 164 result *= x2; 165 166 result += 1./(1.*2); 167 result *= x2; 168 169 result += 1; 170 171 /* Apply correct sign */ 172 result *= cos_sign_tbl[quadrant]; 173 174 return result; 175 } 176 177