home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / SIN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  7.9 KB  |  368 lines
  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 P(x**2).
  26.  *
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  *                      Relative error:
  31.  * arithmetic   range       # trials      peak         rms
  32.  *    DEC       0, 10        20000       2.5e-17     7.1e-18
  33.  *    DEC      0, 1.07e9     10000       2.8e-17     7.2e-18
  34.  *    IEEE      0, 10        35000       2.1e-16     5.4e-17
  35.  *    IEEE -1.07e9,+1.07e9   30000       2.1e-16     5.5e-17
  36.  * 
  37.  * ERROR MESSAGES:
  38.  *
  39.  *   message           condition        value returned
  40.  * sin total loss   x > 1.073741824e9      0.0
  41.  *
  42.  * For 43,000 random points between plus and minus 5.36e8, peak
  43.  * relative error was 4.22e-17, rms relative error 9.29e-18.
  44.  *
  45.  * Partial loss of accuracy begins to occur at x = 2**30
  46.  * = 1.074e9.  The loss is not gradual, but jumps suddenly to
  47.  * about 1 part in 10e7.  Results may be meaningless for
  48.  * x > 2**49 = 5.6e14.  The routine as implemented flags a
  49.  * TLOSS error for x > 2**30 and returns 0.0.
  50.  */
  51. /*                            cos.c
  52.  *
  53.  *    Circular cosine
  54.  *
  55.  *
  56.  *
  57.  * SYNOPSIS:
  58.  *
  59.  * double x, y, cos();
  60.  *
  61.  * y = cos( x );
  62.  *
  63.  *
  64.  *
  65.  * DESCRIPTION:
  66.  *
  67.  * Range reduction is into intervals of pi/4.  The reduction
  68.  * error is nearly eliminated by contriving an extended precision
  69.  * modular arithmetic.
  70.  *
  71.  * Two polynomial approximating functions are employed.
  72.  * Between 0 and pi/4 the cosine is approximated by
  73.  *      1  -  x**2 P(x**2).
  74.  * Between pi/4 and pi/2 the sine is represented as
  75.  *      x  +  x**3 P(x**2).
  76.  *
  77.  *
  78.  * ACCURACY:
  79.  *
  80.  *                      Relative error:
  81.  * arithmetic   range       # trials      peak         rms
  82.  *    IEEE -1.07e9,+1.07e9    5000       2.0e-16     5.5e-17
  83.  *    DEC        0,+1.07e9   17000       3.0e-17     7.2e-18
  84.  *  See also sin().
  85.  *
  86.  */
  87.  
  88. /*                            sin.c    */
  89.  
  90. /* Cephes Math Library Release 2.0:  April, 1987
  91.  * Copyright 1985, 1987 by Stephen L. Moshier
  92.  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
  93.  
  94. #include "mconf.h"
  95.  
  96. #ifdef UNK
  97. static double sincof[] = {
  98.  1.58962301572218447952E-10,
  99. -2.50507477628503540135E-8,
  100.  2.75573136213856773549E-6,
  101. -1.98412698295895384658E-4,
  102.  8.33333333332211858862E-3,
  103. -1.66666666666666307295E-1
  104. };
  105. static double coscof[] = {
  106.  1.13678171382044553091E-11,
  107. -2.08758833757683644217E-9,
  108.  2.75573155429816611547E-7,
  109. -2.48015872936186303776E-5,
  110.  1.38888888888806666760E-3,
  111. -4.16666666666666348141E-2,
  112.  4.99999999999999999798E-1
  113. };
  114. static double DP1 =   7.85398125648498535156E-1;
  115. static double DP2 =   3.77489470793079817668E-8;
  116. static double DP3 =   2.69515142907905952645E-15;
  117. static double lossth = 1.073741824e9;
  118. #endif
  119.  
  120. #ifdef DEC
  121. static short sincof[] = {
  122. 0030056,0143750,0177170,0073013,
  123. 0131727,0027455,0044510,0132205,
  124. 0033470,0167432,0131752,0042263,
  125. 0135120,0006400,0146776,0174027,
  126. 0036410,0104210,0104207,0137202,
  127. 0137452,0125252,0125252,0125103
  128. };
  129. static short coscof[] = {
  130. 0027107,0176030,0153315,0110312,
  131. 0131017,0072476,0007450,0123243,
  132. 0032623,0171174,0070066,0146445,
  133. 0134320,0006400,0147355,0163313,
  134. 0035666,0005540,0133012,0165067,
  135. 0137052,0125252,0125252,0125206,
  136. 0040000,0000000,0000000,0000000
  137. };
  138. /*  7.853981629014015197753906250000E-1 */
  139. static short P1[] = {0040111,0007732,0120000,0000000,};
  140. /*  4.960467869796758577649598009884E-10 */
  141. static short P2[] = {0030410,0055060,0100000,0000000,};
  142. /*  2.860594363054915898381331279295E-18 */
  143. static short P3[] = {0021523,0011431,0105056,0001560,};
  144. #define DP1 *(double *)P1
  145. #define DP2 *(double *)P2
  146. #define DP3 *(double *)P3
  147. static double lossth = 1.073741824e9;
  148. #endif
  149.  
  150. #ifdef IBMPC
  151. static short sincof[] = {
  152. 0x0ec1,0x1fcf,0xd8fd,0x3de5,
  153. 0x1691,0xa929,0xe5e5,0xbe5a,
  154. 0x4896,0x567d,0x1de3,0x3ec7,
  155. 0xdf03,0x19bf,0x01a0,0xbf2a,
  156. 0xf7d0,0x1110,0x1111,0x3f81,
  157. 0x5548,0x5555,0x5555,0xbfc5
  158. };
  159. static short coscof[] = {
  160. 0xb219,0x1ad9,0xff83,0x3da8,
  161. 0x14d4,0xc1e5,0xeea7,0xbe21,
  162. 0xd9a5,0x8e06,0x7e4f,0x3e92,
  163. 0xbcd9,0x19dd,0x01a0,0xbefa,
  164. 0x5d47,0x16c1,0xc16c,0x3f56,
  165. 0x5551,0x5555,0x5555,0xbfa5,
  166. 0x0000,0x0000,0x0000,0x3fe0
  167. };
  168.  
  169. /*
  170.   7.85398125648498535156E-1,
  171.   3.77489470793079817668E-8,
  172.   2.69515142907905952645E-15,
  173. */
  174. static short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
  175. static short P2[] = {0x0000,0x0000,0x442d,0x3e64};
  176. static short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
  177. #define DP1 *(double *)P1
  178. #define DP2 *(double *)P2
  179. #define DP3 *(double *)P3
  180. static double lossth = 1.073741824e9;
  181. #endif
  182.  
  183. #ifdef MIEEE
  184. static short sincof[] = {
  185. 0x3de5,0xd8fd,0x1fcf,0x0ec1,
  186. 0xbe5a,0xe5e5,0xa929,0x1691,
  187. 0x3ec7,0x1de3,0x567d,0x4896,
  188. 0xbf2a,0x01a0,0x19bf,0xdf03,
  189. 0x3f81,0x1111,0x1110,0xf7d0,
  190. 0xbfc5,0x5555,0x5555,0x5548
  191. };
  192. static short coscof[] = {
  193. 0x3da8,0xff83,0x1ad9,0xb219,
  194. 0xbe21,0xeea7,0xc1e5,0x14d4,
  195. 0x3e92,0x7e4f,0x8e06,0xd9a5,
  196. 0xbefa,0x01a0,0x19dd,0xbcd9,
  197. 0x3f56,0xc16c,0x16c1,0x5d47,
  198. 0xbfa5,0x5555,0x5555,0x5551,
  199. 0x3fe0,0x0000,0x0000,0x0000
  200. };
  201.  
  202. static short P1[] = {
  203. 0x3fe9,0x21fb,0x4000,0x0000
  204. };
  205. static short P2[] = {
  206. 0x3e64,0x442d,0x0000,0x0000
  207. };
  208. static short P3[] = {
  209. 0x3ce8,0x4698,0x98cc,0x5170
  210. };
  211. #define DP1 *(double *)P1
  212. #define DP2 *(double *)P2
  213. #define DP3 *(double *)P3
  214. static double lossth = 1.073741824e9;
  215. #endif
  216.  
  217. extern double PIO4;
  218.  
  219. double sin(x)
  220. double x;
  221. {
  222. double y, z, zz;
  223. int rflg, j, sign;
  224. double polevl(), floor(), ldexp();
  225.  
  226. /* make argument positive but save the sign */
  227. sign = 1;
  228. if( x < 0 )
  229.     {
  230.     x = -x;
  231.     sign = -1;
  232.     }
  233.  
  234. if( x > lossth )
  235.     {
  236.     mtherr( "sin", TLOSS );
  237.     return(0.0);
  238.     }
  239.  
  240. y = floor( x/PIO4 ); /* integer part of x/PIO4 */
  241.  
  242. /* strip high bits of integer part to prevent integer overflow */
  243. z = ldexp( y, -4 );
  244. z = floor(z);           /* integer part of y/8 */
  245. z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
  246.  
  247. j = z; /* convert to integer for tests on the phase angle */
  248. /* map zeros to origin */
  249. if( j & 1 )
  250.     {
  251.     j += 1;
  252.     y += 1.0;
  253.     }
  254. j = j & 07; /* octant modulo 360 degrees */
  255. /* reflect in x axis */
  256. if( j > 3)
  257.     {
  258.     sign = -sign;
  259.     j -= 4;
  260.     }
  261.  
  262. /* Extended precision modular arithmetic */
  263. z = ((x - y * DP1) - y * DP2) - y * DP3;
  264.  
  265. zz = z * z;
  266.  
  267. if( (j==1) || (j==2) )
  268.     {
  269.     y = 1.0 - zz * polevl( zz, coscof, 6 );
  270.     }
  271. else
  272.     {
  273.     y = z  +  z * (zz * polevl( zz, sincof, 5 ));
  274.     }
  275.  
  276. if(sign < 0)
  277.     y = -y;
  278.  
  279. return(y);
  280. }
  281.  
  282.  
  283.  
  284.  
  285.  
  286. double cos(x)
  287. double x;
  288. {
  289. double y, z, zz;
  290. long i;
  291. int j, sign, refl;
  292. double polevl(), floor(), ldexp();
  293.  
  294.  
  295. /* make argument positive */
  296. sign = 1;
  297. if( x < 0 )
  298.     x = -x;
  299.  
  300. if( x > lossth )
  301.     {
  302.     mtherr( "cos", TLOSS );
  303.     return(0.0);
  304.     }
  305.  
  306. y = floor( x/PIO4 );
  307. z = ldexp( y, -4 );
  308. z = floor(z);        /* integer part of y/8 */
  309. z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
  310.  
  311. /* integer and fractional part modulo one octant */
  312. i = z;
  313. if( i & 1 )    /* map zeros to origin */
  314.     {
  315.     i += 1;
  316.     y += 1.0;
  317.     }
  318. j = i & 07;
  319. if( j > 3)
  320.     {
  321.     j -=4;
  322.     sign = -sign;
  323.     }
  324.  
  325. if( j > 1 )
  326.     sign = -sign;
  327.  
  328. /* Extended precision modular arithmetic */
  329. z = ((x - y * DP1) - y * DP2) - y * DP3;
  330.  
  331. zz = z * z;
  332.  
  333. if( (j==1) || (j==2) )
  334.     {
  335.     y = z  +  z * (zz * polevl( zz, sincof, 5 ));
  336.     }
  337. else
  338.     {
  339.     y = 1.0 - zz * polevl( zz, coscof, 6 );
  340.     }
  341.  
  342. if(sign < 0)
  343.     y = -y;
  344.  
  345. return(y);
  346. }
  347.  
  348.  
  349.  
  350.  
  351.  
  352. /* Degrees, minutes, seconds to radians: */
  353.  
  354. /* 1 arc second, in radians = 4.8481368110953599358991410e-5 */
  355. #ifdef DEC
  356. static short P648[] = {034513,054170,0176773,0116043,};
  357. #define P64800 *(double *)P648
  358. #else
  359. static double P64800 = 4.8481368110953599358991410e-5;
  360. #endif
  361.  
  362. double radian(d,m,s)
  363. double d,m,s;
  364. {
  365.  
  366. return( ((d*60.0 + m)*60.0 + s)*P64800 );
  367. }
  368.