Subversion Repositories or1k

/* @(#)z_asine.c 1.0 98/08/13 */

/******************************************************************

 * The following routines are coded directly from the algorithms

 * and coefficients given in "Software Manual for the Elementary

 * Functions" by William J. Cody, Jr. and William Waite, Prentice

 * Hall, 1980.

 ******************************************************************/

/*

FUNCTION

        <<asin>>, <<asinf>>, <<acos>>, <<acosf>>, <<asine>>, <<asinef>>---arc sine or cosine

INDEX

   asin

INDEX

   asinf

INDEX

   acos

INDEX

   acosf

INDEX

   asine

INDEX

   asinef

ANSI_SYNOPSIS

        #include <math.h>

        double asine(double <[x]>);

        float asinef(float <[x]>);

        double asin(double <[x]>);

        float asinf(float <[x]>);

        double acos(double <[x]>);

        float acosf(float <[x]>);

TRAD_SYNOPSIS

        #include <math.h>

        double asine(<[x]>);

        double <[x]>;

        float asinef(<[x]>);

        float <[x]>;

        double asin(<[x]>)

        double <[x]>;

        float asinf(<[x]>)

        float <[x]>;

        double acos(<[x]>)

        double <[x]>;

        float acosf(<[x]>)

        float <[x]>;

DESCRIPTION

<<asin>> computes the inverse sine or cosine of the argument <[x]>.

Arguments to <<asin>> and <<acos>> must be in the range @minus{}1 to 1.

<<asinf>> and <<acosf>> are identical to <<asin>> and <<acos>>, other

than taking and returning floats.

RETURNS

@ifinfo

<<asin>> and <<acos>> return values in radians, in the range of -pi/2 to pi/2.

@end ifinfo

@tex

<<asin>> and <<acos>> return values in radians, in the range of $-\pi/2$ to $\pi/2$.

@end tex

If <[x]> is not in the range @minus{}1 to 1, <<asin>> and <<asinf>>

return NaN (not a number), set the global variable <<errno>> to

<<EDOM>>, and issue a <<DOMAIN error>> message.

*/

/******************************************************************

 * Arcsine

 *

 * Input:

 *   x - floating point value

 *   acosine - indicates acos calculation

 *

 * Output:

 *   Arcsine of x.

 *

 * Description:

 *   This routine calculates arcsine / arccosine.

 *

 *****************************************************************/

#include "fdlibm.h"

#include "zmath.h"

#ifndef _DOUBLE_IS_32BITS

static const double p[] = { -0.27368494524164255994e+2,

                             0.57208227877891731407e+2,

                            -0.39688862997404877339e+2,

                             0.10152522233806463645e+2,

                            -0.69674573447350646411 };

static const double q[] = { -0.16421096714498560795e+3,

                             0.41714430248260412556e+3,

                            -0.38186303361750149284e+3,

                             0.15095270841030604719e+3,

                            -0.23823859153670238830e+2 };

static const double a[] = { 0.0, 0.78539816339744830962 };

static const double b[] = { 1.57079632679489661923, 0.78539816339744830962 };

double

_DEFUN (asine, (double, int),

        double x _AND

        int acosine)

{

  int flag, i;

  int branch = 0;

  double g, res, R, P, Q, y;

  /* Check for special values. */

  i = numtest (x);

  if (i == NAN || i == INF)

    {

      errno = EDOM;

      if (i == NAN)

        return (x);

      else

        return (z_infinity.d);

    }

  y = fabs (x);

  flag = acosine;

  if (y > 0.5)

    {

      i = 1 - flag;

      /* Check for range error. */

      if (y > 1.0)

        {

          errno = ERANGE;

          return (z_notanum.d);

        }

      g = (1 - y) / 2.0;

      y = -2 * sqrt (g);

      branch = 1;

    }

  else

    {

      i = flag;

      if (y < z_rooteps)

        res = y;

      else

        g = y * y;

    }

  if (y >= z_rooteps || branch == 1)

    {

      /* Calculate the Taylor series. */

      P = ((((p[4] * g + p[3]) * g + p[2]) * g + p[1]) * g + p[0]) * g;

      Q = ((((g + q[4]) * g + q[3]) * g + q[2]) * g + q[1]) * g + q[0];

      R = P / Q;

      res = y + y * R;

    }

  /* Calculate asine or acose. */

  if (flag == 0)

    {

      res = (a[i] + res) + a[i];

      if (x < 0.0)

        res = -res;

    }

  else

    {

      if (x < 0.0)

        res = (b[i] + res) + b[i];

      else

        res = (a[i] - res) + a[i];

    }

  return (res);

}

#endif /* _DOUBLE_IS_32BITS */