home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / sincos.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  8.6 KB  |  362 lines
  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.  * which must be greater than -360.5 degrees and less than
  20.  * +360.5 degrees.  You must supply a modulo 360 function
  21.  * externally for arguments outside this range.  A compile
  22.  * time option for a built in modulo function is supplied,
  23.  * however.
  24.  *
  25.  * Several alternate algorithms are supplied to achieve
  26.  * various tradeoffs between accuracy and speed.
  27.  *
  28.  * sin(i) is internally tabulated for 0 <= i <= 90 degrees.
  29.  * Approximation polynomials compute the sine and cosine
  30.  * of the residual, which is between -0.5 and +0.5 degree.
  31.  * The results are combined with the tabulated values
  32.  * using the trigonometry 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 6e-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, they
  63.  * 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 approximations are designed with a relative error
  84.  * criterion.  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.  
  105. /* Define one of the following to be 1:
  106.  */
  107. #define ACC5 1
  108. #define ACC11 0
  109. #define ACC17 0
  110.  
  111. /* Option for linear interpolation when flg = 1
  112.  */
  113. #define LINTERP 0
  114.  
  115. /* Option for absolute error criterion
  116.  */
  117. #define ABSERR 1
  118.  
  119. /* Option to include modulo 360 function:
  120.  */
  121. #define MOD360 0
  122.  
  123. /*
  124. Cephes Math Library Release 2.1
  125. Copyright 1987 by Stephen L. Moshier
  126. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  127. */
  128.  
  129.  
  130. /* Table of sin(i degrees)
  131.  * for 0 <= i <= 90
  132.  */
  133. static double sintbl[92] = {
  134.   0.00000000000000000000E0,
  135.   1.74524064372835128194E-2,
  136.   3.48994967025009716460E-2,
  137.   5.23359562429438327221E-2,
  138.   6.97564737441253007760E-2,
  139.   8.71557427476581735581E-2,
  140.   1.04528463267653471400E-1,
  141.   1.21869343405147481113E-1,
  142.   1.39173100960065444112E-1,
  143.   1.56434465040230869010E-1,
  144.   1.73648177666930348852E-1,
  145.   1.90808995376544812405E-1,
  146.   2.07911690817759337102E-1,
  147.   2.24951054343864998051E-1,
  148.   2.41921895599667722560E-1,
  149.   2.58819045102520762349E-1,
  150.   2.75637355816999185650E-1,
  151.   2.92371704722736728097E-1,
  152.   3.09016994374947424102E-1,
  153.   3.25568154457156668714E-1,
  154.   3.42020143325668733044E-1,
  155.   3.58367949545300273484E-1,
  156.   3.74606593415912035415E-1,
  157.   3.90731128489273755062E-1,
  158.   4.06736643075800207754E-1,
  159.   4.22618261740699436187E-1,
  160.   4.38371146789077417453E-1,
  161.   4.53990499739546791560E-1,
  162.   4.69471562785890775959E-1,
  163.   4.84809620246337029075E-1,
  164.   5.00000000000000000000E-1,
  165.   5.15038074910054210082E-1,
  166.   5.29919264233204954047E-1,
  167.   5.44639035015027082224E-1,
  168.   5.59192903470746830160E-1,
  169.   5.73576436351046096108E-1,
  170.   5.87785252292473129169E-1,
  171.   6.01815023152048279918E-1,
  172.   6.15661475325658279669E-1,
  173.   6.29320391049837452706E-1,
  174.   6.42787609686539326323E-1,
  175.   6.56059028990507284782E-1,
  176.   6.69130606358858213826E-1,
  177.   6.81998360062498500442E-1,
  178.   6.94658370458997286656E-1,
  179.   7.07106781186547524401E-1,
  180.   7.19339800338651139356E-1,
  181.   7.31353701619170483288E-1,
  182.   7.43144825477394235015E-1,
  183.   7.54709580222771997943E-1,
  184.   7.66044443118978035202E-1,
  185.   7.77145961456970879980E-1,
  186.   7.88010753606721956694E-1,
  187.   7.98635510047292846284E-1,
  188.   8.09016994374947424102E-1,
  189.   8.19152044288991789684E-1,
  190.   8.29037572555041692006E-1,
  191.   8.38670567945424029638E-1,
  192.   8.48048096156425970386E-1,
  193.   8.57167300702112287465E-1,
  194.   8.66025403784438646764E-1,
  195.   8.74619707139395800285E-1,
  196.   8.82947592858926942032E-1,
  197.   8.91006524188367862360E-1,
  198.   8.98794046299166992782E-1,
  199.   9.06307787036649963243E-1,
  200.   9.13545457642600895502E-1,
  201.   9.20504853452440327397E-1,
  202.   9.27183854566787400806E-1,
  203.   9.33580426497201748990E-1,
  204.   9.39692620785908384054E-1,
  205.   9.45518575599316810348E-1,
  206.   9.51056516295153572116E-1,
  207.   9.56304755963035481339E-1,
  208.   9.61261695938318861916E-1,
  209.   9.65925826289068286750E-1,
  210.   9.70295726275996472306E-1,
  211.   9.74370064785235228540E-1,
  212.   9.78147600733805637929E-1,
  213.   9.81627183447663953497E-1,
  214.   9.84807753012208059367E-1,
  215.   9.87688340595137726190E-1,
  216.   9.90268068741570315084E-1,
  217.   9.92546151641322034980E-1,
  218.   9.94521895368273336923E-1,
  219.   9.96194698091745532295E-1,
  220.   9.97564050259824247613E-1,
  221.   9.98629534754573873784E-1,
  222.   9.99390827019095730006E-1,
  223.   9.99847695156391239157E-1,
  224.   1.00000000000000000000E0,
  225.   9.99847695156391239157E-1,
  226. };
  227.  
  228.  
  229. sincos(x, s, c, flg)
  230. double x;
  231. double *s, *c;
  232. int flg;
  233. {
  234. int ix, ssign, csign, xsign;
  235. double y, z, sx, sz, cx, cz;
  236. #if MOD360
  237. double floor();
  238. #endif
  239.  
  240. /* Make argument nonnegative.
  241.  */
  242. xsign = 1;
  243. if( x < 0.0 )
  244.     {
  245.     xsign = -1;
  246.     x = -x;
  247.     }
  248.  
  249.  
  250. #if MOD360
  251. x = x  -  360.0 * floor( x/360.0 );
  252. #endif
  253.  
  254. /* Find nearest integer to x.
  255.  * Note there should be a domain error test here,
  256.  * but this is omitted to gain speed.
  257.  */
  258. ix = x + 0.5;
  259. z = x - ix;        /* the residual */
  260.  
  261. /* Look up the sine and cosine of the integer.
  262.  */
  263. if( ix <= 180 )
  264.     {
  265.     ssign = 1;
  266.     csign = 1;
  267.     }
  268. else
  269.     {
  270.     ssign = -1;
  271.     csign = -1;
  272.     ix -= 180;
  273.     }
  274.  
  275. if( ix > 90 )
  276.     {
  277.     csign = -csign;
  278.     ix = 180 - ix;
  279.     }
  280.  
  281. sx = sintbl[ix];
  282. if( ssign < 0 )
  283.     sx = -sx;
  284. cx = sintbl[ 90-ix ];
  285. if( csign < 0 )
  286.     cx = -cx;
  287.  
  288. /* If the flag argument is set, then just return
  289.  * the tabulated values for arg to the nearest whole degree.
  290.  */
  291. if( flg )
  292.     {
  293. #if LINTERP
  294.  
  295.     y = sx + 1.74531263774940077459e-2 * z * cx;
  296.     cx -= 1.74531263774940077459e-2 * z * sx;
  297.     sx = y;
  298. #endif
  299.     if( xsign < 0 )
  300.         sx = -sx;
  301.     *s = sx;
  302.     *c = cx;
  303.     return;
  304.     }
  305.  
  306.  
  307. if( ssign < 0 )
  308.     sx = -sx;
  309. if( csign < 0 )
  310.     cx = -cx;
  311.  
  312. /* Find sine and cosine
  313.  * of the residual angle between -0.5 and +0.5 degree.
  314.  */
  315. #if ACC5
  316. #if ABSERR
  317. /* absolute error = 2.769e-8: */
  318. sz = 1.74531263774940077459e-2 * z;
  319. /* absolute error = 4.146e-11: */
  320. cz = 1.0 - 1.52307909153324666207e-4 * z * z;
  321. #else
  322. /* relative error = 6.346e-6: */
  323. sz = 1.74531817576426662296e-2 * z;
  324. /* relative error = 3.173e-6: */
  325. cz = 1.0 - 1.52308226602566149927e-4 * z * z;
  326. #endif
  327. #else
  328. y = z * z;
  329. #endif
  330.  
  331.  
  332. #if ACC11
  333. sz = ( -8.86092781698004819918e-7 * y
  334.       + 1.74532925198378577601e-2     ) * z;
  335.  
  336. cz = 1.0 - ( -3.86631403698859047896e-9 * y
  337.             + 1.52308709893047593702e-4     ) * y;
  338. #endif
  339.  
  340.  
  341. #if ACC17
  342. sz = ((  1.34959795251974073996e-11 * y
  343.        - 8.86096155697856783296e-7     ) * y
  344.        + 1.74532925199432957214e-2          ) * z;
  345.  
  346. cz = 1.0 - ((  3.92582397764340914444e-14 * y
  347.              - 3.86632385155548605680e-9     ) * y
  348.              + 1.52308709893354299569e-4          ) * y;
  349. #endif
  350.  
  351.  
  352. /* Combine the tabulated part and the calculated part
  353.  * by trigonometry.
  354.  */
  355. y = sx * cz  +  cx * sz;
  356. if( xsign < 0 )
  357.     y = - y;
  358. *s = y;
  359.  
  360. *c = cx * cz  -  sx * sz;
  361. }
  362.