newlib/libm/mathfp/sf_asine.c

   1
   2 /* @(#)z_asinef.c 1.0 98/08/13 */
   3 /******************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  ******************************************************************/
   9 /******************************************************************
  10  * Arcsine
  11  *
  12  * Input:
  13  *   x - floating point value
  14  *   acosine - indicates acos calculation
  15  *
  16  * Output:
  17  *   Arcsine of x.
  18  *
  19  * Description:
  20  *   This routine calculates arcsine / arccosine.
  21  *
  22  *****************************************************************/
  23
  24 #include "fdlibm.h"
  25 #include "zmath.h"
  26
  27 static const float p[] = { 0.933935835, -0.504400557 };
  28 static const float q[] = { 0.560363004e+1, -0.554846723e+1 };
  29 static const float a[] = { 0.0, 0.785398163 };
  30 static const float b[] = { 1.570796326, 0.785398163 };
  31
  32 float
  33 asinef (float x,
  34         int acosine)
  35 {
  36   int flag, i;
  37   int branch = 0;
  38   float g, res, R, P, Q, y;
  39
  40   /* Check for special values. */
  41   i = numtestf (x);
  42   if (i == NAN || i == INF)
  43     {
  44       errno = EDOM;
  45       if (i == NAN)
  46         return (x);
  47       else
  48         return (z_infinity_f.f);
  49     }
  50
  51   y = fabsf (x);
  52   flag = acosine;
  53
  54   if (y > 0.5)
  55     {
  56       i = 1 - flag;
  57
  58       /* Check for range error. */
  59       if (y > 1.0)
  60         {
  61           errno = ERANGE;
  62           return (z_notanum_f.f);
  63         }
  64
  65       g = (1 - y) / 2.0;
  66       y = -2 * sqrt (g);
  67       branch = 1;
  68     }
  69   else
  70     {
  71       i = flag;
  72       if (y < z_rooteps_f)
  73         res = y;
  74       else
  75         g = y * y;
  76     }
  77
  78   if (y >= z_rooteps_f || branch == 1)
  79     {
  80       /* Calculate the Taylor series. */
  81       P = (p[1] * g + p[0]) * g;
  82       Q = (g + q[1]) * g + q[0];
  83       R = P / Q;
  84
  85       res = y + y * R;
  86     }
  87
  88   /* Calculate asine or acose. */
  89   if (flag == 0)
  90     {
  91       res = (a[i] + res) + a[i];
  92       if (x < 0.0)
  93         res = -res;
  94     }
  95   else
  96     {
  97       if (x < 0.0)
  98         res = (b[i] + res) + b[i];
  99       else
 100         res = (a[i] - res) + a[i];
 101     }
 102
 103   return (res);
 104 }