|
|
/* @(#)z_sine.c 1.0 98/08/13 */
|
/* @(#)z_sine.c 1.0 98/08/13 */
|
/******************************************************************
|
/******************************************************************
|
* The following routines are coded directly from the algorithms
|
* The following routines are coded directly from the algorithms
|
* and coefficients given in "Software Manual for the Elementary
|
* and coefficients given in "Software Manual for the Elementary
|
* Functions" by William J. Cody, Jr. and William Waite, Prentice
|
* Functions" by William J. Cody, Jr. and William Waite, Prentice
|
* Hall, 1980.
|
* Hall, 1980.
|
******************************************************************/
|
******************************************************************/
|
|
|
/*
|
/*
|
FUNCTION
|
FUNCTION
|
<<sin>>, <<cos>>, <<sine>>, <<sinf>>, <<cosf>>, <<sinef>>---sine or cosine
|
<<sin>>, <<cos>>, <<sine>>, <<sinf>>, <<cosf>>, <<sinef>>---sine or cosine
|
INDEX
|
INDEX
|
sin
|
sin
|
INDEX
|
INDEX
|
sinf
|
sinf
|
INDEX
|
INDEX
|
cos
|
cos
|
INDEX
|
INDEX
|
cosf
|
cosf
|
ANSI_SYNOPSIS
|
ANSI_SYNOPSIS
|
#include <math.h>
|
#include <math.h>
|
double sin(double <[x]>);
|
double sin(double <[x]>);
|
float sinf(float <[x]>);
|
float sinf(float <[x]>);
|
double cos(double <[x]>);
|
double cos(double <[x]>);
|
float cosf(float <[x]>);
|
float cosf(float <[x]>);
|
|
|
TRAD_SYNOPSIS
|
TRAD_SYNOPSIS
|
#include <math.h>
|
#include <math.h>
|
double sin(<[x]>)
|
double sin(<[x]>)
|
double <[x]>;
|
double <[x]>;
|
float sinf(<[x]>)
|
float sinf(<[x]>)
|
float <[x]>;
|
float <[x]>;
|
|
|
double cos(<[x]>)
|
double cos(<[x]>)
|
double <[x]>;
|
double <[x]>;
|
float cosf(<[x]>)
|
float cosf(<[x]>)
|
float <[x]>;
|
float <[x]>;
|
|
|
DESCRIPTION
|
DESCRIPTION
|
<<sin>> and <<cos>> compute (respectively) the sine and cosine
|
<<sin>> and <<cos>> compute (respectively) the sine and cosine
|
of the argument <[x]>. Angles are specified in radians.
|
of the argument <[x]>. Angles are specified in radians.
|
RETURNS
|
RETURNS
|
The sine or cosine of <[x]> is returned.
|
The sine or cosine of <[x]> is returned.
|
|
|
PORTABILITY
|
PORTABILITY
|
<<sin>> and <<cos>> are ANSI C.
|
<<sin>> and <<cos>> are ANSI C.
|
<<sinf>> and <<cosf>> are extensions.
|
<<sinf>> and <<cosf>> are extensions.
|
|
|
QUICKREF
|
QUICKREF
|
sin ansi pure
|
sin ansi pure
|
sinf - pure
|
sinf - pure
|
*/
|
*/
|
|
|
/******************************************************************
|
/******************************************************************
|
* sine
|
* sine
|
*
|
*
|
* Input:
|
* Input:
|
* x - floating point value
|
* x - floating point value
|
* cosine - indicates cosine value
|
* cosine - indicates cosine value
|
*
|
*
|
* Output:
|
* Output:
|
* Sine of x.
|
* Sine of x.
|
*
|
*
|
* Description:
|
* Description:
|
* This routine calculates sines and cosines.
|
* This routine calculates sines and cosines.
|
*
|
*
|
*****************************************************************/
|
*****************************************************************/
|
|
|
#include "fdlibm.h"
|
#include "fdlibm.h"
|
#include "zmath.h"
|
#include "zmath.h"
|
|
|
#ifndef _DOUBLE_IS_32BITS
|
#ifndef _DOUBLE_IS_32BITS
|
|
|
static const double HALF_PI = 1.57079632679489661923;
|
static const double HALF_PI = 1.57079632679489661923;
|
static const double ONE_OVER_PI = 0.31830988618379067154;
|
static const double ONE_OVER_PI = 0.31830988618379067154;
|
static const double r[] = { -0.16666666666666665052,
|
static const double r[] = { -0.16666666666666665052,
|
0.83333333333331650314e-02,
|
0.83333333333331650314e-02,
|
-0.19841269841201840457e-03,
|
-0.19841269841201840457e-03,
|
0.27557319210152756119e-05,
|
0.27557319210152756119e-05,
|
-0.25052106798274584544e-07,
|
-0.25052106798274584544e-07,
|
0.16058936490371589114e-09,
|
0.16058936490371589114e-09,
|
-0.76429178068910467734e-12,
|
-0.76429178068910467734e-12,
|
0.27204790957888846175e-14 };
|
0.27204790957888846175e-14 };
|
|
|
double
|
double
|
_DEFUN (sine, (double, int),
|
_DEFUN (sine, (double, int),
|
double x _AND
|
double x _AND
|
int cosine)
|
int cosine)
|
{
|
{
|
int sgn, N;
|
int sgn, N;
|
double y, XN, g, R, res;
|
double y, XN, g, R, res;
|
double YMAX = 210828714.0;
|
double YMAX = 210828714.0;
|
|
|
switch (numtest (x))
|
switch (numtest (x))
|
{
|
{
|
case NAN:
|
case NAN:
|
errno = EDOM;
|
errno = EDOM;
|
return (x);
|
return (x);
|
case INF:
|
case INF:
|
errno = EDOM;
|
errno = EDOM;
|
return (z_notanum.d);
|
return (z_notanum.d);
|
}
|
}
|
|
|
/* Use sin and cos properties to ease computations. */
|
/* Use sin and cos properties to ease computations. */
|
if (cosine)
|
if (cosine)
|
{
|
{
|
sgn = 1;
|
sgn = 1;
|
y = fabs (x) + HALF_PI;
|
y = fabs (x) + HALF_PI;
|
}
|
}
|
else
|
else
|
{
|
{
|
if (x < 0.0)
|
if (x < 0.0)
|
{
|
{
|
sgn = -1;
|
sgn = -1;
|
y = -x;
|
y = -x;
|
}
|
}
|
else
|
else
|
{
|
{
|
sgn = 1;
|
sgn = 1;
|
y = x;
|
y = x;
|
}
|
}
|
}
|
}
|
|
|
/* Check for values of y that will overflow here. */
|
/* Check for values of y that will overflow here. */
|
if (y > YMAX)
|
if (y > YMAX)
|
{
|
{
|
errno = ERANGE;
|
errno = ERANGE;
|
return (x);
|
return (x);
|
}
|
}
|
|
|
/* Calculate the exponent. */
|
/* Calculate the exponent. */
|
if (y < 0.0)
|
if (y < 0.0)
|
N = (int) (y * ONE_OVER_PI - 0.5);
|
N = (int) (y * ONE_OVER_PI - 0.5);
|
else
|
else
|
N = (int) (y * ONE_OVER_PI + 0.5);
|
N = (int) (y * ONE_OVER_PI + 0.5);
|
XN = (double) N;
|
XN = (double) N;
|
|
|
if (N & 1)
|
if (N & 1)
|
sgn = -sgn;
|
sgn = -sgn;
|
|
|
if (cosine)
|
if (cosine)
|
XN -= 0.5;
|
XN -= 0.5;
|
|
|
y = fabs (x) - XN * __PI;
|
y = fabs (x) - XN * __PI;
|
|
|
if (-z_rooteps < y && y < z_rooteps)
|
if (-z_rooteps < y && y < z_rooteps)
|
res = y;
|
res = y;
|
|
|
else
|
else
|
{
|
{
|
g = y * y;
|
g = y * y;
|
|
|
/* Calculate the Taylor series. */
|
/* Calculate the Taylor series. */
|
R = (((((((r[6] * g + r[5]) * g + r[4]) * g + r[3]) * g + r[2]) * g + r[1]) * g + r[0]) * g);
|
R = (((((((r[6] * g + r[5]) * g + r[4]) * g + r[3]) * g + r[2]) * g + r[1]) * g + r[0]) * g);
|
|
|
/* Finally, compute the result. */
|
/* Finally, compute the result. */
|
res = y + y * R;
|
res = y + y * R;
|
}
|
}
|
|
|
res *= sgn;
|
res *= sgn;
|
|
|
return (res);
|
return (res);
|
}
|
}
|
|
|
#endif /* _DOUBLE_IS_32BITS */
|
#endif /* _DOUBLE_IS_32BITS */
|
|
|
