source: trunk/libs/newlib/src/newlib/libm/mathfp/s_sine.c @ 452

Last change on this file since 452 was 444, checked in by satin@…, 6 years ago

add newlib,libalmos-mkh, restructure shared_syscalls.h and mini-libc
File size: 3.2 KB
Line 
1
2/* @(#)z_sine.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/*
11FUNCTION
12        <<sin>>, <<cos>>, <<sine>>, <<sinf>>, <<cosf>>, <<sinef>>---sine or cosine
13INDEX
14sin
15INDEX
16sinf
17INDEX
18cos
19INDEX
20cosf
21SYNOPSIS
22        #include <math.h>
23        double sin(double <[x]>);
24        float  sinf(float <[x]>);
25        double cos(double <[x]>);
26        float cosf(float <[x]>);
27
28DESCRIPTION
29        <<sin>> and <<cos>> compute (respectively) the sine and cosine
30        of the argument <[x]>.  Angles are specified in radians.
31RETURNS
32        The sine or cosine of <[x]> is returned.
33
34PORTABILITY
35        <<sin>> and <<cos>> are ANSI C.
36        <<sinf>> and <<cosf>> are extensions.
37
38QUICKREF
39        sin ansi pure
40        sinf - pure
41*/
42
43/******************************************************************
44 * sine
45 *
46 * Input:
47 *   x - floating point value
48 *   cosine - indicates cosine value
49 *
50 * Output:
51 *   Sine of x.
52 *
53 * Description:
54 *   This routine calculates sines and cosines.
55 *
56 *****************************************************************/
57
58#include "fdlibm.h"
59#include "zmath.h"
60
61#ifndef _DOUBLE_IS_32BITS
62
63static const double HALF_PI = 1.57079632679489661923;
64static const double ONE_OVER_PI = 0.31830988618379067154;
65static const double r[] = { -0.16666666666666665052,
66                             0.83333333333331650314e-02,
67                            -0.19841269841201840457e-03,
68                             0.27557319210152756119e-05,
69                            -0.25052106798274584544e-07,
70                             0.16058936490371589114e-09,
71                            -0.76429178068910467734e-12,
72                             0.27204790957888846175e-14 };
73
74double
75sine (double x,
76        int cosine)
77{
78  int sgn, N;
79  double y, XN, g, R, res;
80  double YMAX = 210828714.0;
81
82  switch (numtest (x))
83    {
84      case NAN:
85        errno = EDOM;
86        return (x);
87      case INF:
88        errno = EDOM;
89        return (z_notanum.d); 
90    }
91
92  /* Use sin and cos properties to ease computations. */
93  if (cosine)
94    {
95      sgn = 1;
96      y = fabs (x) + HALF_PI;
97    }
98  else
99    {
100      if (x < 0.0)
101        {
102          sgn = -1;
103          y = -x;
104        }
105      else
106        {
107          sgn = 1;
108          y = x;
109        }
110    }
111
112  /* Check for values of y that will overflow here. */
113  if (y > YMAX)
114    {
115      errno = ERANGE;
116      return (x);
117    }
118
119  /* Calculate the exponent. */
120  if (y < 0.0)
121    N = (int) (y * ONE_OVER_PI - 0.5);
122  else
123    N = (int) (y * ONE_OVER_PI + 0.5);
124  XN = (double) N;
125
126  if (N & 1)
127    sgn = -sgn;
128
129  if (cosine)
130    XN -= 0.5;
131
132  y = fabs (x) - XN * __PI;
133
134  if (-z_rooteps < y && y < z_rooteps)
135    res = y;
136
137  else
138    {
139      g = y * y;
140
141      /* Calculate the Taylor series. */
142      R = (((((((r[6] * g + r[5]) * g + r[4]) * g + r[3]) * g + r[2]) * g + r[1]) * g + r[0]) * g);
143
144      /* Finally, compute the result. */
145      res = y + y * R;
146    }
147 
148  res *= sgn;
149
150  return (res);
151}
152
153#endif /* _DOUBLE_IS_32BITS */
Note: See TracBrowser for help on using the repository browser.