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