home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ldouble / asinl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  4.1 KB  |  232 lines
  1. /*                            asinl.c
  2.  *
  3.  *    Inverse circular sine, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, asinl();
  10.  *
  11.  * y = asinl( 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   domain     # trials      peak         rms
  30.  *    IEEE     -1, 1        30000       2.7e-19     4.8e-20
  31.  *
  32.  *
  33.  * ERROR MESSAGES:
  34.  *
  35.  *   message         condition      value returned
  36.  * asin domain        |x| > 1           0.0
  37.  *
  38.  */
  39. /*                            acosl()
  40.  *
  41.  *    Inverse circular cosine, long double precision
  42.  *
  43.  *
  44.  *
  45.  * SYNOPSIS:
  46.  *
  47.  * double x, y, acosl();
  48.  *
  49.  * y = acosl( x );
  50.  *
  51.  *
  52.  *
  53.  * DESCRIPTION:
  54.  *
  55.  * Returns radian angle between -pi/2 and +pi/2 whose cosine
  56.  * is x.
  57.  *
  58.  * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
  59.  * near 1, there is cancellation error in subtracting asin(x)
  60.  * from pi/2.  Hence if x < -0.5,
  61.  *
  62.  *    acos(x) =     pi - 2.0 * asin( sqrt((1+x)/2) );
  63.  *
  64.  * or if x > +0.5,
  65.  *
  66.  *    acos(x) =     2.0 * asin(  sqrt((1-x)/2) ).
  67.  *
  68.  *
  69.  * ACCURACY:
  70.  *
  71.  *                      Relative error:
  72.  * arithmetic   domain     # trials      peak         rms
  73.  *    IEEE      -1, 1       30000       1.4e-19     3.5e-20
  74.  *
  75.  *
  76.  * ERROR MESSAGES:
  77.  *
  78.  *   message         condition      value returned
  79.  * asin domain        |x| > 1           0.0
  80.  */
  81.  
  82. /*                            asin.c    */
  83.  
  84. /*
  85. Cephes Math Library Release 2.2:  December, 1990
  86. Copyright 1984, 1990 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 long double P[] = {
  94.  3.7769340062433674871612E-3L,
  95. -6.1212919176969202969441E-1L,
  96.  5.9303993515791417710775E0L,
  97. -1.8631697621590161441592E1L,
  98.  2.3314603132141795720634E1L,
  99. -1.0087146579384916260197E1L,
  100. };
  101. static long double Q[] = {
  102. /* 1.0000000000000000000000E0L,*/
  103. -1.5684335624873146511217E1L,
  104.  7.8702951549021104258866E1L,
  105. -1.7078401170625864261444E2L,
  106.  1.6712291455718995937376E2L,
  107. -6.0522879476309497128868E1L,
  108. };
  109. #endif
  110.  
  111. #ifdef IBMPC
  112. static short P[] = {
  113. 0x59d1,0x3509,0x7009,0xf786,0x3ff6,
  114. 0xbe97,0x93e6,0x7fab,0x9cb4,0xbffe,
  115. 0x8bf5,0x6810,0xd4dc,0xbdc5,0x4001,
  116. 0x9bd4,0x8d86,0xb77b,0x950d,0xc003,
  117. 0x3b0f,0x9e25,0x4ea5,0xba84,0x4003,
  118. 0xea38,0xc6a9,0xf3cf,0xa164,0xc002,
  119. };
  120. static short Q[] = {
  121. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  122. 0x1229,0x8516,0x09e9,0xfaf3,0xc002,
  123. 0xb5c3,0xf36f,0xe943,0x9d67,0x4005,
  124. 0xe11a,0xbe0f,0xb4fd,0xaac8,0xc006,
  125. 0x4c69,0x1355,0x7754,0xa71f,0x4006,
  126. 0xded7,0xa9fe,0x6db7,0xf217,0xc004,
  127. };
  128. #endif
  129.  
  130. #ifdef MIEEE
  131. static long P[] = {
  132. 0x3ff60000,0xf7867009,0x350959d1,
  133. 0xbffe0000,0x9cb47fab,0x93e6be97,
  134. 0x40010000,0xbdc5d4dc,0x68108bf5,
  135. 0xc0030000,0x950db77b,0x8d869bd4,
  136. 0x40030000,0xba844ea5,0x9e253b0f,
  137. 0xc0020000,0xa164f3cf,0xc6a9ea38,
  138. };
  139. static long Q[] = {
  140. /*0x3fff0000,0x80000000,0x00000000,*/
  141. 0xc0020000,0xfaf309e9,0x85161229,
  142. 0x40050000,0x9d67e943,0xf36fb5c3,
  143. 0xc0060000,0xaac8b4fd,0xbe0fe11a,
  144. 0x40060000,0xa71f7754,0x13554c69,
  145. 0xc0040000,0xf2176db7,0xa9feded7,
  146. };
  147. #endif
  148.  
  149. long double asinl(x)
  150. long double x;
  151. {
  152. long double a, p, z, zz;
  153. long double ldexpl(), sqrtl(), polevll(), p1evll();
  154. short sign, flag;
  155. extern long double PIO2L;
  156.  
  157. if( x > 0 )
  158.     {
  159.     sign = 1;
  160.     a = x;
  161.     }
  162. else
  163.     {
  164.     sign = -1;
  165.     a = -x;
  166.     }
  167.  
  168. if( a > 1.0L )
  169.     {
  170.     mtherr( "asinl", DOMAIN );
  171.     return( 0.0L );
  172.     }
  173.  
  174. if( a < 1.0e-8L )
  175.     {
  176.     z = a;
  177.     goto done;
  178.     }
  179.  
  180. if( a > 0.5L )
  181.     {
  182.     zz = 0.5L -a;
  183.     zz = ldexpl( zz + 0.5L, -1 );
  184.     z = sqrtl( zz );
  185.     flag = 1;
  186.     }
  187. else
  188.     {
  189.     z = a;
  190.     zz = z * z;
  191.     flag = 0;
  192.     }
  193.  
  194. p = zz * polevll( zz, P, 5)/p1evll( zz, Q, 5);
  195. z = z * p + z;
  196. if( flag != 0 )
  197.     {
  198.     z = z + z;
  199.     z = PIO2L - z;
  200.     }
  201. done:
  202. if( sign < 0 )
  203.     z = -z;
  204. return(z);
  205. }
  206.  
  207.  
  208. extern long double PIO2L, PIL;
  209.  
  210. long double acosl(x)
  211. long double x;
  212. {
  213. long double asinl(), sqrtl();
  214.  
  215. if( x < -1.0L )
  216.     goto domerr;
  217.  
  218. if( x < -0.5L) 
  219.     return( PIL - 2.0L * asinl( sqrtl(0.5L*(1.0L+x)) ) );
  220.  
  221. if( x > 1.0L )
  222.     {
  223. domerr:    mtherr( "acosl", DOMAIN );
  224.     return( 0.0L );
  225.     }
  226.  
  227. if( x > 0.5L )
  228.     return( 2.0L * asinl(  sqrtl(0.5L*(1.0L-x) ) ) );
  229.  
  230. return( PIO2L - asinl(x) );
  231. }
  232.