home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / ASIN.C next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  4.5 KB  |  248 lines
  1. /*                            asin.c
  2.  *
  3.  *    Inverse circular sine
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, asin();
  10.  *
  11.  * y = asin( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
  18.  *
  19.  * A rational function of the form x + x**3 P(x**2)/Q(x**2)
  20.  * is used for |x| in the interval [0, 0.5].  If |x| > 0.5 it is
  21.  * transformed by the identity
  22.  *
  23.  *    asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *                      Relative error:
  29.  * arithmetic   range      # trials      peak         rms
  30.  *    DEC      -1, 1        20000       5.7e-17     1.2e-17
  31.  *    IEEE     -1, 1        10000       4.6e-16     8.4e-17
  32.  *
  33.  *
  34.  * ERROR MESSAGES:
  35.  *
  36.  *   message         condition      value returned
  37.  * asin domain        |x| > 1           0.0
  38.  *
  39.  */
  40. /*                            acos()
  41.  *
  42.  *    Inverse circular cosine
  43.  *
  44.  *
  45.  *
  46.  * SYNOPSIS:
  47.  *
  48.  * double x, y, acos();
  49.  *
  50.  * y = acos( x );
  51.  *
  52.  *
  53.  *
  54.  * DESCRIPTION:
  55.  *
  56.  * Returns radian angle between -pi/2 and +pi/2 whose cosine
  57.  * is x.
  58.  *
  59.  * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
  60.  * near 1, there is cancellation error in subtracting asin(x)
  61.  * from pi/2.  Hence if x < -0.5,
  62.  *
  63.  *    acos(x) =     pi - 2.0 * asin( sqrt((1+x)/2) );
  64.  *
  65.  * or if x > +0.5,
  66.  *
  67.  *    acos(x) =     2.0 * asin(  sqrt((1-x)/2) ).
  68.  *
  69.  *
  70.  * ACCURACY:
  71.  *
  72.  *                      Relative error:
  73.  * arithmetic   range      # trials      peak         rms
  74.  *    DEC       -1, 1       13000       3.2e-17     7.8e-18
  75.  *    IEEE      -1, 1       10000       2.1e-16     6.7e-17
  76.  *
  77.  *
  78.  * ERROR MESSAGES:
  79.  *
  80.  *   message         condition      value returned
  81.  * asin domain        |x| > 1           0.0
  82.  */
  83.  
  84. /*                            asin.c    */
  85.  
  86. /*
  87. Cephes Math Library Release 2.0:  April, 1987
  88. Copyright 1984, 1987 by Stephen L. Moshier
  89. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  90. */
  91.  
  92. #include "mconf.h"
  93.  
  94. #ifdef UNK
  95. static double P[] = {
  96. -6.96822599948686174217E-1,
  97.  1.01598286089872099722E1,
  98. -3.97340771391578757294E1,
  99.  5.72912144709846496134E1,
  100. -2.74148200465925708020E1
  101. };
  102. static double Q[] = {
  103. /* 1.00000000000000000000E0,*/
  104. -2.38368245005177488242E1,
  105.  1.51095072703128995631E2,
  106. -3.82340216045978957023E2,
  107.  4.17767300951716199422E2,
  108. -1.64488920279555473283E2
  109. };
  110. #endif
  111.  
  112. #ifdef DEC
  113. static short P[] = {
  114. 0140062,0061367,0042744,0127464,
  115. 0041042,0107250,0070611,0005470,
  116. 0141436,0167661,0165345,0131200,
  117. 0041545,0025064,0020124,0000411,
  118. 0141333,0050615,0026056,0134351
  119. };
  120. static short Q[] = {
  121. /*0040200,0000000,0000000,0000000,*/
  122. 0141276,0130721,0005461,0134336,
  123. 0042027,0014126,0127506,0127155,
  124. 0142277,0025614,0031413,0103353,
  125. 0042320,0161066,0165346,0163707,
  126. 0142044,0076451,0160443,0005274
  127. };
  128. #endif
  129.  
  130. #ifdef IBMPC
  131. static short P[] = {
  132. 0x95e7,0xe8bc,0x4c5e,0xbfe6,
  133. 0x2167,0x0e31,0x51d5,0x4024,
  134. 0xb650,0x3d5c,0xddf6,0xc043,
  135. 0x8021,0x840a,0xa546,0x404c,
  136. 0xd71d,0xa585,0x6a31,0xc03b
  137. };
  138. static short Q[] = {
  139. /*0x0000,0x0000,0x0000,0x3ff0,*/
  140. 0x371c,0x2166,0xd63a,0xc037,
  141. 0xd5ce,0xd5e8,0xe30a,0x4062,
  142. 0x70dd,0x8661,0xe571,0xc077,
  143. 0xdcf9,0xdd5c,0x1c46,0x407a,
  144. 0x6158,0x3c24,0x8fa5,0xc064
  145. };
  146. #endif
  147.  
  148. #ifdef MIEEE
  149. static short P[] = {
  150. 0xbfe6,0x4c5e,0xe8bc,0x95e7,
  151. 0x4024,0x51d5,0x0e31,0x2167,
  152. 0xc043,0xddf6,0x3d5c,0xb650,
  153. 0x404c,0xa546,0x840a,0x8021,
  154. 0xc03b,0x6a31,0xa585,0xd71d
  155. };
  156. static short Q[] = {
  157. 0xc037,0xd63a,0x2166,0x371c,
  158. 0x4062,0xe30a,0xd5e8,0xd5ce,
  159. 0xc077,0xe571,0x8661,0x70dd,
  160. 0x407a,0x1c46,0xdd5c,0xdcf9,
  161. 0xc0640x8fa5,0x3c24,0x6158
  162. };
  163. #endif
  164.  
  165. double asin(x)
  166. double x;
  167. {
  168. double a, p, z, zz;
  169. double sqrt(), polevl(), p1evl();
  170. short sign, flag;
  171. extern double PIO2;
  172.  
  173. if( x > 0 )
  174.     {
  175.     sign = 1;
  176.     a = x;
  177.     }
  178. else
  179.     {
  180.     sign = -1;
  181.     a = -x;
  182.     }
  183.  
  184. if( a > 1.0 )
  185.     {
  186.     mtherr( "asin", DOMAIN );
  187.     return( 0.0 );
  188.     }
  189.  
  190. if( a < 1.0e-7 )
  191.     {
  192.     z = a;
  193.     goto done;
  194.     }
  195.  
  196. if( a > 0.5 )
  197.     {
  198.     zz = 0.5 -a;
  199.     zz = (zz + 0.5)/2.0;
  200.     z = sqrt( zz );
  201.     flag = 1;
  202.     }
  203. else
  204.     {
  205.     z = a;
  206.     zz = z * z;
  207.     flag = 0;
  208.     }
  209.  
  210. p = zz * polevl( zz, P, 4)/p1evl( zz, Q, 5);
  211. z = z * p + z;
  212. if( flag != 0 )
  213.     {
  214.     z = z + z;
  215.     z = PIO2 - z;
  216.     }
  217. done:
  218. if( sign < 0 )
  219.     z = -z;
  220. return(z);
  221. }
  222.  
  223.  
  224. extern double PIO2, PI;
  225.  
  226. double acos(x)
  227. double x;
  228. {
  229. double asin(), sqrt();
  230.  
  231. if( x < -1.0 )
  232.     goto domerr;
  233.  
  234. if( x < -0.5) 
  235.     return( PI - 2.0 * asin( sqrt(0.5*(1.0+x)) ) );
  236.  
  237. if( x > 1.0 )
  238.     {
  239. domerr:    mtherr( "acos", DOMAIN );
  240.     return( 0.0 );
  241.     }
  242.  
  243. if( x > 0.5 )
  244.     return( 2.0 * asin(  sqrt(0.5*(1.0-x) ) ) );
  245.  
  246. return( PIO2 - asin(x) );
  247. }
  248.