libm/double/sincos.c

   1 /*                                                      sincos.c
   2  *
   3  *      Circular sine and cosine of argument in degrees
   4  *      Table lookup and interpolation algorithm
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * double x, sine, cosine, flg, sincos();
  11  *
  12  * sincos( x, &sine, &cosine, flg );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns both the sine and the cosine of the argument x.
  19  * Several different compile time options and minimax
  20  * approximations are supplied to permit tailoring the
  21  * tradeoff between computation speed and accuracy.
  22  *
  23  * Since range reduction is time consuming, the reduction
  24  * of x modulo 360 degrees is also made optional.
  25  *
  26  * sin(i) is internally tabulated for 0 <= i <= 90 degrees.
  27  * Approximation polynomials, ranging from linear interpolation
  28  * to cubics in (x-i)**2, compute the sine and cosine
  29  * of the residual x-i which is between -0.5 and +0.5 degree.
  30  * In the case of the high accuracy options, the residual
  31  * and the tabulated values are combined using the trigonometry
  32  * formulas for sin(A+B) and cos(A+B).
  33  *
  34  * Compile time options are supplied for 5, 11, or 17 decimal
  35  * relative accuracy (ACC5, ACC11, ACC17 respectively).
  36  * A subroutine flag argument "flg" chooses betwen this
  37  * accuracy and table lookup only (peak absolute error
  38  * = 0.0087).
  39  *
  40  * If the argument flg = 1, then the tabulated value is
  41  * returned for the nearest whole number of degrees. The
  42  * approximation polynomials are not computed.  At
  43  * x = 0.5 deg, the absolute error is then sin(0.5) = 0.0087.
  44  *
  45  * An intermediate speed and precision can be obtained using
  46  * the compile time option LINTERP and flg = 1.  This yields
  47  * a linear interpolation using a slope estimated from the sine
  48  * or cosine at the nearest integer argument.  The peak absolute
  49  * error with this option is 3.8e-5.  Relative error at small
  50  * angles is about 1e-5.
  51  *
  52  * If flg = 0, then the approximation polynomials are computed
  53  * and applied.
  54  *
  55  *
  56  *
  57  * SPEED:
  58  *
  59  * Relative speed comparisons follow for 6MHz IBM AT clone
  60  * and Microsoft C version 4.0.  These figures include
  61  * software overhead of do loop and function calls.
  62  * Since system hardware and software vary widely, the
  63  * numbers should be taken as representative only.
  64  *
  65  *                      flg=0   flg=0   flg=1   flg=1
  66  *                      ACC11   ACC5    LINTERP Lookup only
  67  * In-line 8087 (/FPi)
  68  * sin(), cos()         1.0     1.0     1.0     1.0
  69  *
  70  * In-line 8087 (/FPi)
  71  * sincos()             1.1     1.4     1.9     3.0
  72  *
  73  * Software (/FPa)
  74  * sin(), cos()         0.19    0.19    0.19    0.19
  75  *
  76  * Software (/FPa)
  77  * sincos()             0.39    0.50    0.73    1.7
  78  *
  79  *
  80  *
  81  * ACCURACY:
  82  *
  83  * The accurate approximations are designed with a relative error
  84  * criterion.  The absolute error is greatest at x = 0.5 degree.
  85  * It decreases from a local maximum at i+0.5 degrees to full
  86  * machine precision at each integer i degrees.  With the
  87  * ACC5 option, the relative error of 6.3e-6 is equivalent to
  88  * an absolute angular error of 0.01 arc second in the argument
  89  * at x = i+0.5 degrees.  For small angles < 0.5 deg, the ACC5
  90  * accuracy is 6.3e-6 (.00063%) of reading; i.e., the absolute
  91  * error decreases in proportion to the argument.  This is true
  92  * for both the sine and cosine approximations, since the latter
  93  * is for the function 1 - cos(x).
  94  *
  95  * If absolute error is of most concern, use the compile time
  96  * option ABSERR to obtain an absolute error of 2.7e-8 for ACC5
  97  * precision.  This is about half the absolute error of the
  98  * relative precision option.  In this case the relative error
  99  * for small angles will increase to 9.5e-6 -- a reasonable
 100  * tradeoff.
 101  */
 102
 103
 104 #include <math.h>
 105
 106 /* Define one of the following to be 1:
 107  */
 108 #define ACC5 1
 109 #define ACC11 0
 110 #define ACC17 0
 111
 112 /* Option for linear interpolation when flg = 1
 113  */
 114 #define LINTERP 1
 115
 116 /* Option for absolute error criterion
 117  */
 118 #define ABSERR 1
 119
 120 /* Option to include modulo 360 function:
 121  */
 122 #define MOD360 0
 123
 124 /*
 125 Cephes Math Library Release 2.1
 126 Copyright 1987 by Stephen L. Moshier
 127 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 128 */
 129
 130
 131 /* Table of sin(i degrees)
 132  * for 0 <= i <= 90
 133  */
 134 static double sintbl[92] = {
 135   0.00000000000000000000E0,
 136   1.74524064372835128194E-2,
 137   3.48994967025009716460E-2,
 138   5.23359562429438327221E-2,
 139   6.97564737441253007760E-2,
 140   8.71557427476581735581E-2,
 141   1.04528463267653471400E-1,
 142   1.21869343405147481113E-1,
 143   1.39173100960065444112E-1,
 144   1.56434465040230869010E-1,
 145   1.73648177666930348852E-1,
 146   1.90808995376544812405E-1,
 147   2.07911690817759337102E-1,
 148   2.24951054343864998051E-1,
 149   2.41921895599667722560E-1,
 150   2.58819045102520762349E-1,
 151   2.75637355816999185650E-1,
 152   2.92371704722736728097E-1,
 153   3.09016994374947424102E-1,
 154   3.25568154457156668714E-1,
 155   3.42020143325668733044E-1,
 156   3.58367949545300273484E-1,
 157   3.74606593415912035415E-1,
 158   3.90731128489273755062E-1,
 159   4.06736643075800207754E-1,
 160   4.22618261740699436187E-1,
 161   4.38371146789077417453E-1,
 162   4.53990499739546791560E-1,
 163   4.69471562785890775959E-1,
 164   4.84809620246337029075E-1,
 165   5.00000000000000000000E-1,
 166   5.15038074910054210082E-1,
 167   5.29919264233204954047E-1,
 168   5.44639035015027082224E-1,
 169   5.59192903470746830160E-1,
 170   5.73576436351046096108E-1,
 171   5.87785252292473129169E-1,
 172   6.01815023152048279918E-1,
 173   6.15661475325658279669E-1,
 174   6.29320391049837452706E-1,
 175   6.42787609686539326323E-1,
 176   6.56059028990507284782E-1,
 177   6.69130606358858213826E-1,
 178   6.81998360062498500442E-1,
 179   6.94658370458997286656E-1,
 180   7.07106781186547524401E-1,
 181   7.19339800338651139356E-1,
 182   7.31353701619170483288E-1,
 183   7.43144825477394235015E-1,
 184   7.54709580222771997943E-1,
 185   7.66044443118978035202E-1,
 186   7.77145961456970879980E-1,
 187   7.88010753606721956694E-1,
 188   7.98635510047292846284E-1,
 189   8.09016994374947424102E-1,
 190   8.19152044288991789684E-1,
 191   8.29037572555041692006E-1,
 192   8.38670567945424029638E-1,
 193   8.48048096156425970386E-1,
 194   8.57167300702112287465E-1,
 195   8.66025403784438646764E-1,
 196   8.74619707139395800285E-1,
 197   8.82947592858926942032E-1,
 198   8.91006524188367862360E-1,
 199   8.98794046299166992782E-1,
 200   9.06307787036649963243E-1,
 201   9.13545457642600895502E-1,
 202   9.20504853452440327397E-1,
 203   9.27183854566787400806E-1,
 204   9.33580426497201748990E-1,
 205   9.39692620785908384054E-1,
 206   9.45518575599316810348E-1,
 207   9.51056516295153572116E-1,
 208   9.56304755963035481339E-1,
 209   9.61261695938318861916E-1,
 210   9.65925826289068286750E-1,
 211   9.70295726275996472306E-1,
 212   9.74370064785235228540E-1,
 213   9.78147600733805637929E-1,
 214   9.81627183447663953497E-1,
 215   9.84807753012208059367E-1,
 216   9.87688340595137726190E-1,
 217   9.90268068741570315084E-1,
 218   9.92546151641322034980E-1,
 219   9.94521895368273336923E-1,
 220   9.96194698091745532295E-1,
 221   9.97564050259824247613E-1,
 222   9.98629534754573873784E-1,
 223   9.99390827019095730006E-1,
 224   9.99847695156391239157E-1,
 225   1.00000000000000000000E0,
 226   9.99847695156391239157E-1,
 227 };
 228
 229 #ifdef ANSIPROT
 230 double floor ( double );
 231 #else
 232 double floor();
 233 #endif
 234
 235 int sincos(x, s, c, flg)
 236 double x;
 237 double *s, *c;
 238 int flg;
 239 {
 240 int ix, ssign, csign, xsign;
 241 double y, z, sx, sz, cx, cz;
 242
 243 /* Make argument nonnegative.
 244  */
 245 xsign = 1;
 246 if( x < 0.0 )
 247         {
 248         xsign = -1;
 249         x = -x;
 250         }
 251
 252
 253 #if MOD360
 254 x = x  -  360.0 * floor( x/360.0 );
 255 #endif
 256
 257 /* Find nearest integer to x.
 258  * Note there should be a domain error test here,
 259  * but this is omitted to gain speed.
 260  */
 261 ix = x + 0.5;
 262 z = x - ix;             /* the residual */
 263
 264 /* Look up the sine and cosine of the integer.
 265  */
 266 if( ix <= 180 )
 267         {
 268         ssign = 1;
 269         csign = 1;
 270         }
 271 else
 272         {
 273         ssign = -1;
 274         csign = -1;
 275         ix -= 180;
 276         }
 277
 278 if( ix > 90 )
 279         {
 280         csign = -csign;
 281         ix = 180 - ix;
 282         }
 283
 284 sx = sintbl[ix];
 285 if( ssign < 0 )
 286         sx = -sx;
 287 cx = sintbl[ 90-ix ];
 288 if( csign < 0 )
 289         cx = -cx;
 290
 291 /* If the flag argument is set, then just return
 292  * the tabulated values for arg to the nearest whole degree.
 293  */
 294 if( flg )
 295         {
 296 #if LINTERP
 297         y = sx + 1.74531263774940077459e-2 * z * cx;
 298         cx -= 1.74531263774940077459e-2 * z * sx;
 299         sx = y;
 300 #endif
 301         if( xsign < 0 )
 302                 sx = -sx;
 303         *s = sx;        /* sine */
 304         *c = cx;        /* cosine */
 305         return 0;
 306         }
 307
 308
 309 if( ssign < 0 )
 310         sx = -sx;
 311 if( csign < 0 )
 312         cx = -cx;
 313
 314 /* Find sine and cosine
 315  * of the residual angle between -0.5 and +0.5 degree.
 316  */
 317 #if ACC5
 318 #if ABSERR
 319 /* absolute error = 2.769e-8: */
 320 sz = 1.74531263774940077459e-2 * z;
 321 /* absolute error = 4.146e-11: */
 322 cz = 1.0 - 1.52307909153324666207e-4 * z * z;
 323 #else
 324 /* relative error = 6.346e-6: */
 325 sz = 1.74531817576426662296e-2 * z;
 326 /* relative error = 3.173e-6: */
 327 cz = 1.0 - 1.52308226602566149927e-4 * z * z;
 328 #endif
 329 #else
 330 y = z * z;
 331 #endif
 332
 333
 334 #if ACC11
 335 sz = ( -8.86092781698004819918e-7 * y
 336       + 1.74532925198378577601e-2     ) * z;
 337
 338 cz = 1.0 - ( -3.86631403698859047896e-9 * y
 339             + 1.52308709893047593702e-4     ) * y;
 340 #endif
 341
 342
 343 #if ACC17
 344 sz = ((  1.34959795251974073996e-11 * y
 345        - 8.86096155697856783296e-7     ) * y
 346        + 1.74532925199432957214e-2          ) * z;
 347
 348 cz = 1.0 - ((  3.92582397764340914444e-14 * y
 349              - 3.86632385155548605680e-9     ) * y
 350              + 1.52308709893354299569e-4          ) * y;
 351 #endif
 352
 353
 354 /* Combine the tabulated part and the calculated part
 355  * by trigonometry.
 356  */
 357 y = sx * cz  +  cx * sz;
 358 if( xsign < 0 )
 359         y = - y;
 360 *s = y; /* sine */
 361
 362 *c = cx * cz  -  sx * sz; /* cosine */
 363 return 0;
 364 }