home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / SINL.C < prev    next >
C/C++ Source or Header  |  1993-01-30  |  7KB  |  303 lines
  1. /*                            sinl.c
  2.  *
  3.  *    Circular sine, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, sinl();
  10.  *
  11.  * y = sinl( 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 the Cody
  23.  * and Waite polynomial form
  24.  *      x + x**3 P(x**2) .
  25.  * Between pi/4 and pi/2 the cosine is represented as
  26.  *      1 - .5 x**2 + x**4 Q(x**2) .
  27.  *
  28.  *
  29.  * ACCURACY:
  30.  *
  31.  *                      Relative error:
  32.  * arithmetic   domain      # trials      peak         rms
  33.  *    IEEE     +-5.5e11      200,000    1.2e-19     2.9e-20
  34.  * 
  35.  * ERROR MESSAGES:
  36.  *
  37.  *   message           condition        value returned
  38.  * sin total loss   x > 2**39               0.0
  39.  *
  40.  * Loss of precision occurs for x > 2**39 = 5.49755813888e11.
  41.  * The routine as implemented flags a TLOSS error for
  42.  * x > 2**39 and returns 0.0.
  43.  */
  44. /*                            cosl.c
  45.  *
  46.  *    Circular cosine, long double precision
  47.  *
  48.  *
  49.  *
  50.  * SYNOPSIS:
  51.  *
  52.  * long double x, y, cosl();
  53.  *
  54.  * y = cosl( x );
  55.  *
  56.  *
  57.  *
  58.  * DESCRIPTION:
  59.  *
  60.  * Range reduction is into intervals of pi/4.  The reduction
  61.  * error is nearly eliminated by contriving an extended precision
  62.  * modular arithmetic.
  63.  *
  64.  * Two polynomial approximating functions are employed.
  65.  * Between 0 and pi/4 the cosine is approximated by
  66.  *      1 - .5 x**2 + x**4 Q(x**2) .
  67.  * Between pi/4 and pi/2 the sine is represented by the Cody
  68.  * and Waite polynomial form
  69.  *      x  +  x**3 P(x**2) .
  70.  *
  71.  *
  72.  * ACCURACY:
  73.  *
  74.  *                      Relative error:
  75.  * arithmetic   domain      # trials      peak         rms
  76.  *    IEEE     +-5.5e11       50000      1.2e-19     2.9e-20
  77.  */
  78.  
  79. /*                            sin.c    */
  80.  
  81. /*
  82. Cephes Math Library Release 2.2:  December, 1990
  83. Copyright 1985, 1990 by Stephen L. Moshier
  84. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  85. */
  86.  
  87. #ifndef NOINCS
  88. #include "mconf.h"
  89. #endif
  90.  
  91. #ifdef UNK
  92. long double sincof[7] = {
  93. -7.5785404094842805756289E-13L,
  94.  1.6058363167320443249231E-10L,
  95. -2.5052104881870868784055E-8L,
  96.  2.7557319214064922217861E-6L,
  97. -1.9841269841254799668344E-4L,
  98.  8.3333333333333225058715E-3L,
  99. -1.6666666666666666640255E-1L,
  100. };
  101. long double coscof[7] = {
  102.  4.7377507964246204691685E-14L,
  103. -1.1470284843425359765671E-11L,
  104.  2.0876754287081521758361E-9L,
  105. -2.7557319214999787979814E-7L,
  106.  2.4801587301570552304991E-5L,
  107. -1.3888888888888872993737E-3L,
  108.  4.1666666666666666609054E-2L,
  109. };
  110. long double sinDP1 = 7.853981554508209228515625E-1L;
  111. long double sinDP2 = 7.946627356147928367136046290398E-9L;
  112. long double sinDP3 = 3.061616997868382943065164830688E-17L;
  113. long double sinlossth = 5.49755813888e11L; /* 2^39 */
  114. #endif
  115.  
  116. #ifdef IBMPC
  117. short sincof[] = {
  118. 0x4e27,0xe1d6,0x2389,0xd551,0xbfd6, XPD
  119. 0x64d7,0xe706,0x4623,0xb090,0x3fde, XPD
  120. 0x01b1,0xbf34,0x2946,0xd732,0xbfe5, XPD
  121. 0xc8f7,0x9845,0x1d29,0xb8ef,0x3fec, XPD
  122. 0x6514,0x0c53,0x00d0,0xd00d,0xbff2, XPD
  123. 0x569a,0x8888,0x8888,0x8888,0x3ff8, XPD
  124. 0xaa97,0xaaaa,0xaaaa,0xaaaa,0xbffc, XPD
  125. };
  126. short coscof[] = {
  127. 0x7436,0x6f99,0x8c3a,0xd55e,0x3fd2, XPD
  128. 0x2f37,0x58f4,0x920f,0xc9c9,0xbfda, XPD
  129. 0x5350,0x659e,0xc648,0x8f76,0x3fe2, XPD
  130. 0x4d2b,0xf5c6,0x7dba,0x93f2,0xbfe9, XPD
  131. 0x53ed,0x0c66,0x00d0,0xd00d,0x3fef, XPD
  132. 0x7b67,0x0b60,0x60b6,0xb60b,0xbff5, XPD
  133. 0xaa9a,0xaaaa,0xaaaa,0xaaaa,0x3ffa, XPD
  134. };
  135. short sinP1[] = {0x0000,0x0000,0xda80,0xc90f,0x3ffe, XPD};
  136. short sinP2[] = {0x0000,0x0000,0xa300,0x8885,0x3fe4, XPD};
  137. short sinP3[] = {0x3707,0xa2e0,0x3198,0x8d31,0x3fc8, XPD};
  138. #define sinDP1 *(long double *)sinP1
  139. #define sinDP2 *(long double *)sinP2
  140. #define sinDP3 *(long double *)sinP3
  141. short sinlosi[] = {0x0000,0x0000,0x0000,0x8000,0x4026, XPD};
  142. #define sinlossth *(long double *)sinlosi
  143. #endif
  144.  
  145. #ifdef MIEEE
  146. long sincof[] = {
  147. 0xbfd60000,0xd5512389,0xe1d64e27,
  148. 0x3fde0000,0xb0904623,0xe70664d7,
  149. 0xbfe50000,0xd7322946,0xbf3401b1,
  150. 0x3fec0000,0xb8ef1d29,0x9845c8f7,
  151. 0xbff20000,0xd00d00d0,0x0c536514,
  152. 0x3ff80000,0x88888888,0x8888569a,
  153. 0xbffc0000,0xaaaaaaaa,0xaaaaaa97,
  154. };
  155. long coscof[] = {
  156. 0x3fd20000,0xd55e8c3a,0x6f997436,
  157. 0xbfda0000,0xc9c9920f,0x58f42f37,
  158. 0x3fe20000,0x8f76c648,0x659e5350,
  159. 0xbfe90000,0x93f27dba,0xf5c64d2b,
  160. 0x3fef0000,0xd00d00d0,0x0c6653ed,
  161. 0xbff50000,0xb60b60b6,0x0b607b67,
  162. 0x3ffa0000,0xaaaaaaaa,0xaaaaaa9a,
  163. };
  164. long sinP1[] = {0x3ffe0000,0xc90fda80,0x00000000};
  165. long sinP2[] = {0x3fe40000,0x8885a300,0x00000000};
  166. long sinP3[] = {0x3fc80000,0x8d313198,0xa2e03707};
  167. #define sinDP1 *(long double *)sinP1
  168. #define sinDP2 *(long double *)sinP2
  169. #define sinDP3 *(long double *)sinP3
  170. long sinlosi[] = {0x40260000,0x80000000,0x00000000};
  171. #define sinlossth *(long double *)sinlosi
  172. #endif
  173.  
  174. extern long double PIO4L, Zero, One;
  175.  
  176.  
  177. long double sinl(x)
  178. long double x;
  179. {
  180. long double y, z, zz;
  181. int j, sign;
  182. long double polevll(), floorl(), ldexpl();
  183.  
  184. /* make argument positive but save the sign */
  185. sign = 1;
  186. if( x < 0 )
  187.     {
  188.     x = -x;
  189.     sign = -1;
  190.     }
  191.  
  192. if( x > sinlossth )
  193.     {
  194.     mtherr( "sinl", TLOSS );
  195.     return(Zero);
  196.     }
  197.  
  198. y = floorl( x/PIO4L ); /* integer part of x/PIO4 */
  199.  
  200. /* strip high bits of integer part to prevent integer overflow */
  201. z = ldexpl( y, -4 );
  202. z = floorl(z);           /* integer part of y/8 */
  203. z = y - ldexpl( z, 4 );  /* y - 16 * (y/16) */
  204.  
  205. j = (int )z; /* convert to integer for tests on the phase angle */
  206. /* map zeros to origin */
  207. if( j & 1 )
  208.     {
  209.     j += 1;
  210.     y += One;
  211.     }
  212. j = j & 07; /* octant modulo 360 degrees */
  213. /* reflect in x axis */
  214. if( j > 3)
  215.     {
  216.     sign = -sign;
  217.     j -= 4;
  218.     }
  219.  
  220. /* Extended precision modular arithmetic */
  221. z = ((x - y * sinDP1) - y * sinDP2) - y * sinDP3;
  222.  
  223. zz = z * z;
  224. if( (j==1) || (j==2) )
  225.     {
  226.     y = One - ldexpl(zz,-1) + zz * zz * polevll( zz, coscof, 6 );
  227.     }
  228. else
  229.     {
  230.     y = z  +  z * (zz * polevll( zz, sincof, 6 ));
  231.     }
  232.  
  233. if(sign < 0)
  234.     y = -y;
  235.  
  236. return(y);
  237. }
  238.  
  239.  
  240.  
  241.  
  242.  
  243. long double cosl(x)
  244. long double x;
  245. {
  246. long double y, z, zz;
  247. long i;
  248. int j, sign;
  249. long double polevll(), floorl(), ldexpl();
  250.  
  251.  
  252. /* make argument positive */
  253. sign = 1;
  254. if( x < 0 )
  255.     x = -x;
  256.  
  257. if( x > sinlossth )
  258.     {
  259.     mtherr( "cosl", TLOSS );
  260.     return(Zero);
  261.     }
  262.  
  263. y = floorl( x/PIO4L );
  264. z = ldexpl( y, -4 );
  265. z = floorl(z);        /* integer part of y/8 */
  266. z = y - ldexpl( z, 4 );  /* y - 16 * (y/16) */
  267.  
  268. /* integer and fractional part modulo one octant */
  269. i = (int )z;
  270. if( i & 1 )    /* map zeros to origin */
  271.     {
  272.     i += 1;
  273.     y += One;
  274.     }
  275. j = i & 07;
  276. if( j > 3)
  277.     {
  278.     j -=4;
  279.     sign = -sign;
  280.     }
  281.  
  282. if( j > 1 )
  283.     sign = -sign;
  284.  
  285. /* Extended precision modular arithmetic */
  286. z = ((x - y * sinDP1) - y * sinDP2) - y * sinDP3;
  287.  
  288. zz = z * z;
  289. if( (j==1) || (j==2) )
  290.     {
  291.     y = z  +  z * (zz * polevll( zz, sincof, 6 ));
  292.     }
  293. else
  294.     {
  295.     y = One - ldexpl(zz,-1) + zz * zz * polevll( zz, coscof, 6 );
  296.     }
  297.  
  298. if(sign < 0)
  299.     y = -y;
  300.  
  301. return(y);
  302. }
  303.