lib/libm/sindg.c

   1 /*                                                      sindg.c
   2  *
   3  *      Circular sine of angle in degrees
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, sindg();
  10  *
  11  * y = sindg( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Range reduction is into intervals of 45 degrees.
  18  *
  19  * Two polynomial approximating functions are employed.
  20  * Between 0 and pi/4 the sine is approximated by
  21  *      x  +  x**3 P(x**2).
  22  * Between pi/4 and pi/2 the cosine is represented as
  23  *      1  -  x**2 P(x**2).
  24  *
  25  *
  26  *
  27  * ACCURACY:
  28  *
  29  *                      Relative error:
  30  * arithmetic   domain      # trials      peak         rms
  31  *    DEC       +-1000        3100      3.3e-17      9.0e-18
  32  *    IEEE      +-1000       30000      2.3e-16      5.6e-17
  33  *
  34  * ERROR MESSAGES:
  35  *
  36  *   message           condition        value returned
  37  * sindg total loss   x > 8.0e14 (DEC)      0.0
  38  *                    x > 1.0e14 (IEEE)
  39  *
  40  */
  41 \f/*                                                     cosdg.c
  42  *
  43  *      Circular cosine of angle in degrees
  44  *
  45  *
  46  *
  47  * SYNOPSIS:
  48  *
  49  * double x, y, cosdg();
  50  *
  51  * y = cosdg( x );
  52  *
  53  *
  54  *
  55  * DESCRIPTION:
  56  *
  57  * Range reduction is into intervals of 45 degrees.
  58  *
  59  * Two polynomial approximating functions are employed.
  60  * Between 0 and pi/4 the cosine is approximated by
  61  *      1  -  x**2 P(x**2).
  62  * Between pi/4 and pi/2 the sine is represented as
  63  *      x  +  x**3 P(x**2).
  64  *
  65  *
  66  * ACCURACY:
  67  *
  68  *                      Relative error:
  69  * arithmetic   domain      # trials      peak         rms
  70  *    DEC      +-1000         3400       3.5e-17     9.1e-18
  71  *    IEEE     +-1000        30000       2.1e-16     5.7e-17
  72  *  See also sin().
  73  *
  74  */
  75 \f
  76 /* Cephes Math Library Release 2.0:  April, 1987
  77  * Copyright 1985, 1987 by Stephen L. Moshier
  78  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
  79
  80 #include "mconf.h"
  81
  82 #ifdef UNK
  83 static double sincof[] = {
  84  1.58962301572218447952E-10,
  85 -2.50507477628503540135E-8,
  86  2.75573136213856773549E-6,
  87 -1.98412698295895384658E-4,
  88  8.33333333332211858862E-3,
  89 -1.66666666666666307295E-1
  90 };
  91 static double coscof[] = {
  92  1.13678171382044553091E-11,
  93 -2.08758833757683644217E-9,
  94  2.75573155429816611547E-7,
  95 -2.48015872936186303776E-5,
  96  1.38888888888806666760E-3,
  97 -4.16666666666666348141E-2,
  98  4.99999999999999999798E-1
  99 };
 100 static double PI180 = 1.74532925199432957692E-2; /* pi/180 */
 101 static double lossth = 1.0e14;
 102 #endif
 103
 104 #ifdef DEC
 105 static unsigned short sincof[] = {
 106 0030056,0143750,0177170,0073013,
 107 0131727,0027455,0044510,0132205,
 108 0033470,0167432,0131752,0042263,
 109 0135120,0006400,0146776,0174027,
 110 0036410,0104210,0104207,0137202,
 111 0137452,0125252,0125252,0125103
 112 };
 113 static unsigned short coscof[] = {
 114 0027107,0176030,0153315,0110312,
 115 0131017,0072476,0007450,0123243,
 116 0032623,0171174,0070066,0146445,
 117 0134320,0006400,0147355,0163313,
 118 0035666,0005540,0133012,0165067,
 119 0137052,0125252,0125252,0125206,
 120 0040000,0000000,0000000,0000000
 121 };
 122 static unsigned short P1[] = {0036616,0175065,0011224,0164711};
 123 #define PI180 *(double *)P1
 124 static double lossth = 8.0e14;
 125 #endif
 126
 127 #ifdef IBMPC
 128 static unsigned short sincof[] = {
 129 0x0ec1,0x1fcf,0xd8fd,0x3de5,
 130 0x1691,0xa929,0xe5e5,0xbe5a,
 131 0x4896,0x567d,0x1de3,0x3ec7,
 132 0xdf03,0x19bf,0x01a0,0xbf2a,
 133 0xf7d0,0x1110,0x1111,0x3f81,
 134 0x5548,0x5555,0x5555,0xbfc5
 135 };
 136 static unsigned short coscof[] = {
 137 0xb219,0x1ad9,0xff83,0x3da8,
 138 0x14d4,0xc1e5,0xeea7,0xbe21,
 139 0xd9a5,0x8e06,0x7e4f,0x3e92,
 140 0xbcd9,0x19dd,0x01a0,0xbefa,
 141 0x5d47,0x16c1,0xc16c,0x3f56,
 142 0x5551,0x5555,0x5555,0xbfa5,
 143 0x0000,0x0000,0x0000,0x3fe0
 144 };
 145
 146 static unsigned short P1[] = {0x9d39,0xa252,0xdf46,0x3f91};
 147 #define PI180 *(double *)P1
 148 static double lossth = 1.0e14;
 149 #endif
 150
 151 #ifdef MIEEE
 152 static unsigned short sincof[] = {
 153 0x3de5,0xd8fd,0x1fcf,0x0ec1,
 154 0xbe5a,0xe5e5,0xa929,0x1691,
 155 0x3ec7,0x1de3,0x567d,0x4896,
 156 0xbf2a,0x01a0,0x19bf,0xdf03,
 157 0x3f81,0x1111,0x1110,0xf7d0,
 158 0xbfc5,0x5555,0x5555,0x5548
 159 };
 160 static unsigned short coscof[] = {
 161 0x3da8,0xff83,0x1ad9,0xb219,
 162 0xbe21,0xeea7,0xc1e5,0x14d4,
 163 0x3e92,0x7e4f,0x8e06,0xd9a5,
 164 0xbefa,0x01a0,0x19dd,0xbcd9,
 165 0x3f56,0xc16c,0x16c1,0x5d47,
 166 0xbfa5,0x5555,0x5555,0x5551,
 167 0x3fe0,0x0000,0x0000,0x0000
 168 };
 169
 170 static unsigned short P1[] = {
 171 0x3f91,0xdf46,0xa252,0x9d39
 172 };
 173 #define PI180 *(double *)P1
 174 static double lossth = 1.0e14;
 175 #endif
 176
 177 #ifdef ANSIPROT
 178 extern double polevl ( double, void *, int );
 179 extern double floor ( double );
 180 extern double ldexp ( double, int );
 181 #else
 182 double polevl(), floor(), ldexp();
 183 #endif
 184 extern double PIO4;
 185
 186 double sindg(x)
 187 double x;
 188 {
 189 double y, z, zz;
 190 int j, sign;
 191
 192 /* make argument positive but save the sign */
 193 sign = 1;
 194 if( x < 0 )
 195         {
 196         x = -x;
 197         sign = -1;
 198         }
 199
 200 if( x > lossth )
 201         {
 202         mtherr( "sindg", TLOSS );
 203         return(0.0);
 204         }
 205
 206 y = floor( x/45.0 ); /* integer part of x/PIO4 */
 207
 208 /* strip high bits of integer part to prevent integer overflow */
 209 z = ldexp( y, -4 );
 210 z = floor(z);           /* integer part of y/8 */
 211 z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
 212
 213 j = z; /* convert to integer for tests on the phase angle */
 214 /* map zeros to origin */
 215 if( j & 1 )
 216         {
 217         j += 1;
 218         y += 1.0;
 219         }
 220 j = j & 07; /* octant modulo 360 degrees */
 221 /* reflect in x axis */
 222 if( j > 3)
 223         {
 224         sign = -sign;
 225         j -= 4;
 226         }
 227
 228 z = x - y * 45.0; /* x mod 45 degrees */
 229 z *= PI180;     /* multiply by pi/180 to convert to radians */
 230 zz = z * z;
 231
 232 if( (j==1) || (j==2) )
 233         {
 234         y = 1.0 - zz * polevl( zz, coscof, 6 );
 235         }
 236 else
 237         {
 238         y = z  +  z * (zz * polevl( zz, sincof, 5 ));
 239         }
 240
 241 if(sign < 0)
 242         y = -y;
 243
 244 return(y);
 245 }
 246
 247
 248
 249
 250
 251 double cosdg(x)
 252 double x;
 253 {
 254 double y, z, zz;
 255 int j, sign;
 256
 257 /* make argument positive */
 258 sign = 1;
 259 if( x < 0 )
 260         x = -x;
 261
 262 if( x > lossth )
 263         {
 264         mtherr( "cosdg", TLOSS );
 265         return(0.0);
 266         }
 267
 268 y = floor( x/45.0 );
 269 z = ldexp( y, -4 );
 270 z = floor(z);           /* integer part of y/8 */
 271 z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
 272
 273 /* integer and fractional part modulo one octant */
 274 j = z;
 275 if( j & 1 )     /* map zeros to origin */
 276         {
 277         j += 1;
 278         y += 1.0;
 279         }
 280 j = j & 07;
 281 if( j > 3)
 282         {
 283         j -=4;
 284         sign = -sign;
 285         }
 286
 287 if( j > 1 )
 288         sign = -sign;
 289
 290 z = x - y * 45.0; /* x mod 45 degrees */
 291 z *= PI180;     /* multiply by pi/180 to convert to radians */
 292
 293 zz = z * z;
 294
 295 if( (j==1) || (j==2) )
 296         {
 297         y = z  +  z * (zz * polevl( zz, sincof, 5 ));
 298         }
 299 else
 300         {
 301         y = 1.0 - zz * polevl( zz, coscof, 6 );
 302         }
 303
 304 if(sign < 0)
 305         y = -y;
 306
 307 return(y);
 308 }