lib/libm/sin.c

   1 /*                                                      sin.c
   2  *
   3  *      Circular sine
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, sin();
  10  *
  11  * y = sin( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Range reduction is into intervals of pi/4.  The reduction
  18  * error is nearly eliminated by contriving an extended precision
  19  * modular arithmetic.
  20  *
  21  * Two polynomial approximating functions are employed.
  22  * Between 0 and pi/4 the sine is approximated by
  23  *      x  +  x**3 P(x**2).
  24  * Between pi/4 and pi/2 the cosine is represented as
  25  *      1  -  x**2 Q(x**2).
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  *                      Relative error:
  31  * arithmetic   domain      # trials      peak         rms
  32  *    DEC       0, 10       150000       3.0e-17     7.8e-18
  33  *    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
  34  *
  35  * ERROR MESSAGES:
  36  *
  37  *   message           condition        value returned
  38  * sin total loss   x > 1.073741824e9      0.0
  39  *
  40  * Partial loss of accuracy begins to occur at x = 2**30
  41  * = 1.074e9.  The loss is not gradual, but jumps suddenly to
  42  * about 1 part in 10e7.  Results may be meaningless for
  43  * x > 2**49 = 5.6e14.  The routine as implemented flags a
  44  * TLOSS error for x > 2**30 and returns 0.0.
  45  */
  46 \f/*                                                     cos.c
  47  *
  48  *      Circular cosine
  49  *
  50  *
  51  *
  52  * SYNOPSIS:
  53  *
  54  * double x, y, cos();
  55  *
  56  * y = cos( x );
  57  *
  58  *
  59  *
  60  * DESCRIPTION:
  61  *
  62  * Range reduction is into intervals of pi/4.  The reduction
  63  * error is nearly eliminated by contriving an extended precision
  64  * modular arithmetic.
  65  *
  66  * Two polynomial approximating functions are employed.
  67  * Between 0 and pi/4 the cosine is approximated by
  68  *      1  -  x**2 Q(x**2).
  69  * Between pi/4 and pi/2 the sine is represented as
  70  *      x  +  x**3 P(x**2).
  71  *
  72  *
  73  * ACCURACY:
  74  *
  75  *                      Relative error:
  76  * arithmetic   domain      # trials      peak         rms
  77  *    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
  78  *    DEC        0,+1.07e9   17000       3.0e-17     7.2e-18
  79  */
  80 \f
  81 /*                                                      sin.c   */
  82
  83 /*
  84 Cephes Math Library Release 2.8:  June, 2000
  85 Copyright 1985, 1995, 2000 by Stephen L. Moshier
  86 */
  87
  88 #include "mconf.h"
  89
  90 #ifdef UNK
  91 static double sincof[] = {
  92  1.58962301576546568060E-10,
  93 -2.50507477628578072866E-8,
  94  2.75573136213857245213E-6,
  95 -1.98412698295895385996E-4,
  96  8.33333333332211858878E-3,
  97 -1.66666666666666307295E-1,
  98 };
  99 static double coscof[6] = {
 100 -1.13585365213876817300E-11,
 101  2.08757008419747316778E-9,
 102 -2.75573141792967388112E-7,
 103  2.48015872888517045348E-5,
 104 -1.38888888888730564116E-3,
 105  4.16666666666665929218E-2,
 106 };
 107 static double DP1 =   7.85398125648498535156E-1;
 108 static double DP2 =   3.77489470793079817668E-8;
 109 static double DP3 =   2.69515142907905952645E-15;
 110 /* static double lossth = 1.073741824e9; */
 111 #endif
 112
 113 #ifdef DEC
 114 static unsigned short sincof[] = {
 115 0030056,0143750,0177214,0163153,
 116 0131727,0027455,0044510,0175352,
 117 0033470,0167432,0131752,0042414,
 118 0135120,0006400,0146776,0174027,
 119 0036410,0104210,0104207,0137202,
 120 0137452,0125252,0125252,0125103,
 121 };
 122 static unsigned short coscof[24] = {
 123 0127107,0151115,0002060,0152325,
 124 0031017,0072353,0155161,0174053,
 125 0132623,0171173,0172542,0057056,
 126 0034320,0006400,0147102,0023652,
 127 0135666,0005540,0133012,0076213,
 128 0037052,0125252,0125252,0125126,
 129 };
 130 /*  7.853981629014015197753906250000E-1 */
 131 static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
 132 /*  4.960467869796758577649598009884E-10 */
 133 static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
 134 /*  2.860594363054915898381331279295E-18 */
 135 static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
 136 #define DP1 *(double *)P1
 137 #define DP2 *(double *)P2
 138 #define DP3 *(double *)P3
 139 #endif
 140
 141 #ifdef IBMPC
 142 static unsigned short sincof[] = {
 143 0x9ccd,0x1fd1,0xd8fd,0x3de5,
 144 0x1f5d,0xa929,0xe5e5,0xbe5a,
 145 0x48a1,0x567d,0x1de3,0x3ec7,
 146 0xdf03,0x19bf,0x01a0,0xbf2a,
 147 0xf7d0,0x1110,0x1111,0x3f81,
 148 0x5548,0x5555,0x5555,0xbfc5,
 149 };
 150 static unsigned short coscof[24] = {
 151 0x1a9b,0xa086,0xfa49,0xbda8,
 152 0x3f05,0x7b4e,0xee9d,0x3e21,
 153 0x4bc6,0x7eac,0x7e4f,0xbe92,
 154 0x44f5,0x19c8,0x01a0,0x3efa,
 155 0x4f91,0x16c1,0xc16c,0xbf56,
 156 0x554b,0x5555,0x5555,0x3fa5,
 157 };
 158 /*
 159   7.85398125648498535156E-1,
 160   3.77489470793079817668E-8,
 161   2.69515142907905952645E-15,
 162 */
 163 static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
 164 static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
 165 static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
 166 #define DP1 *(double *)P1
 167 #define DP2 *(double *)P2
 168 #define DP3 *(double *)P3
 169 #endif
 170
 171 #ifdef MIEEE
 172 static unsigned short sincof[] = {
 173 0x3de5,0xd8fd,0x1fd1,0x9ccd,
 174 0xbe5a,0xe5e5,0xa929,0x1f5d,
 175 0x3ec7,0x1de3,0x567d,0x48a1,
 176 0xbf2a,0x01a0,0x19bf,0xdf03,
 177 0x3f81,0x1111,0x1110,0xf7d0,
 178 0xbfc5,0x5555,0x5555,0x5548,
 179 };
 180 static unsigned short coscof[24] = {
 181 0xbda8,0xfa49,0xa086,0x1a9b,
 182 0x3e21,0xee9d,0x7b4e,0x3f05,
 183 0xbe92,0x7e4f,0x7eac,0x4bc6,
 184 0x3efa,0x01a0,0x19c8,0x44f5,
 185 0xbf56,0xc16c,0x16c1,0x4f91,
 186 0x3fa5,0x5555,0x5555,0x554b,
 187 };
 188 static unsigned short P1[] = {0x3fe9,0x21fb,0x4000,0x0000};
 189 static unsigned short P2[] = {0x3e64,0x442d,0x0000,0x0000};
 190 static unsigned short P3[] = {0x3ce8,0x4698,0x98cc,0x5170};
 191 #define DP1 *(double *)P1
 192 #define DP2 *(double *)P2
 193 #define DP3 *(double *)P3
 194 #endif
 195
 196 #ifdef ANSIPROT
 197 extern double polevl ( double, void *, int );
 198 extern double p1evl ( double, void *, int );
 199 extern double floor ( double );
 200 extern double ldexp ( double, int );
 201 extern int isnan ( double );
 202 extern int isfinite ( double );
 203 #else
 204 double polevl(), floor(), ldexp();
 205 int isnan(), isfinite();
 206 #endif
 207 extern double PIO4;
 208 static double lossth = 1.073741824e9;
 209 #ifdef NANS
 210 extern double NAN;
 211 #endif
 212 #ifdef INFINITIES
 213 extern double INFINITY;
 214 #endif
 215
 216
 217 double sin(x)
 218 double x;
 219 {
 220 double y, z, zz;
 221 int j, sign;
 222
 223 #ifdef MINUSZERO
 224 if( x == 0.0 )
 225         return(x);
 226 #endif
 227 #ifdef NANS
 228 if( isnan(x) )
 229         return(x);
 230 if( !isfinite(x) )
 231         {
 232         mtherr( "sin", DOMAIN );
 233         return(NAN);
 234         }
 235 #endif
 236 /* make argument positive but save the sign */
 237 sign = 1;
 238 if( x < 0 )
 239         {
 240         x = -x;
 241         sign = -1;
 242         }
 243
 244 if( x > lossth )
 245         {
 246         mtherr( "sin", TLOSS );
 247         return(0.0);
 248         }
 249
 250 y = floor( x/PIO4 ); /* integer part of x/PIO4 */
 251
 252 /* strip high bits of integer part to prevent integer overflow */
 253 z = ldexp( y, -4 );
 254 z = floor(z);           /* integer part of y/8 */
 255 z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
 256
 257 j = z; /* convert to integer for tests on the phase angle */
 258 /* map zeros to origin */
 259 if( j & 1 )
 260         {
 261         j += 1;
 262         y += 1.0;
 263         }
 264 j = j & 07; /* octant modulo 360 degrees */
 265 /* reflect in x axis */
 266 if( j > 3)
 267         {
 268         sign = -sign;
 269         j -= 4;
 270         }
 271
 272 /* Extended precision modular arithmetic */
 273 z = ((x - y * DP1) - y * DP2) - y * DP3;
 274
 275 zz = z * z;
 276
 277 if( (j==1) || (j==2) )
 278         {
 279         y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
 280         }
 281 else
 282         {
 283 /*      y = z  +  z * (zz * polevl( zz, sincof, 5 ));*/
 284         y = z  +  z * z * z * polevl( zz, sincof, 5 );
 285         }
 286
 287 if(sign < 0)
 288         y = -y;
 289
 290 return(y);
 291 }
 292
 293
 294
 295
 296
 297 double cos(x)
 298 double x;
 299 {
 300 double y, z, zz;
 301 long i;
 302 int j, sign;
 303
 304 #ifdef NANS
 305 if( isnan(x) )
 306         return(x);
 307 if( !isfinite(x) )
 308         {
 309         mtherr( "cos", DOMAIN );
 310         return(NAN);
 311         }
 312 #endif
 313
 314 /* make argument positive */
 315 sign = 1;
 316 if( x < 0 )
 317         x = -x;
 318
 319 if( x > lossth )
 320         {
 321         mtherr( "cos", TLOSS );
 322         return(0.0);
 323         }
 324
 325 y = floor( x/PIO4 );
 326 z = ldexp( y, -4 );
 327 z = floor(z);           /* integer part of y/8 */
 328 z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
 329
 330 /* integer and fractional part modulo one octant */
 331 i = z;
 332 if( i & 1 )     /* map zeros to origin */
 333         {
 334         i += 1;
 335         y += 1.0;
 336         }
 337 j = i & 07;
 338 if( j > 3)
 339         {
 340         j -=4;
 341         sign = -sign;
 342         }
 343
 344 if( j > 1 )
 345         sign = -sign;
 346
 347 /* Extended precision modular arithmetic */
 348 z = ((x - y * DP1) - y * DP2) - y * DP3;
 349
 350 zz = z * z;
 351
 352 if( (j==1) || (j==2) )
 353         {
 354 /*      y = z  +  z * (zz * polevl( zz, sincof, 5 ));*/
 355         y = z  +  z * z * z * polevl( zz, sincof, 5 );
 356         }
 357 else
 358         {
 359         y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
 360         }
 361
 362 if(sign < 0)
 363         y = -y;
 364
 365 return(y);
 366 }
 367
 368
 369
 370
 371
 372 /* Degrees, minutes, seconds to radians: */
 373
 374 /* 1 arc second, in radians = 4.8481368110953599358991410e-5 */
 375 #ifdef DEC
 376 static unsigned short P648[] = {034513,054170,0176773,0116043,};
 377 #define P64800 *(double *)P648
 378 #else
 379 static double P64800 = 4.8481368110953599358991410e-5;
 380 #endif
 381
 382 double radian(d,m,s)
 383 double d,m,s;
 384 {
 385
 386 return( ((d*60.0 + m)*60.0 + s)*P64800 );
 387 }