source: trunk/minix/lib/math/sin.c@ 11

Last change on this file since 11 was 9, checked in by Mattia Monga, 14 years ago

Minix 3.1.2a

File size: 1.8 KB
Line 
1/*
2 * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
3 * See the copyright notice in the ACK home directory, in the file "Copyright".
4 *
5 * Author: Ceriel J.H. Jacobs
6 */
7/* $Header: /cvsup/minix/src/lib/math/sin.c,v 1.1.1.1 2005/04/21 14:56:26 beng Exp $ */
8
9#include <math.h>
10#include <float.h>
11#include <errno.h>
12#include "localmath.h"
13
14static double
15sinus(double x, int cos_flag)
16{
17 /* Algorithm and coefficients from:
18 "Software manual for the elementary functions"
19 by W.J. Cody and W. Waite, Prentice-Hall, 1980
20 */
21
22 static double r[] = {
23 -0.16666666666666665052e+0,
24 0.83333333333331650314e-2,
25 -0.19841269841201840457e-3,
26 0.27557319210152756119e-5,
27 -0.25052106798274584544e-7,
28 0.16058936490371589114e-9,
29 -0.76429178068910467734e-12,
30 0.27204790957888846175e-14
31 };
32
33 double y;
34 int neg = 1;
35
36 if (__IsNan(x)) {
37 errno = EDOM;
38 return x;
39 }
40 if (x < 0) {
41 x = -x;
42 neg = -1;
43 }
44 if (cos_flag) {
45 neg = 1;
46 y = M_PI_2 + x;
47 }
48 else y = x;
49
50 /* ??? avoid loss of significance, if y is too large, error ??? */
51
52 y = y * M_1_PI + 0.5;
53
54 if (y >= DBL_MAX/M_PI) return 0.0;
55
56 /* Use extended precision to calculate reduced argument.
57 Here we used 12 bits of the mantissa for a1.
58 Also split x in integer part x1 and fraction part x2.
59 */
60#define A1 3.1416015625
61#define A2 -8.908910206761537356617e-6
62 {
63 double x1, x2;
64
65 modf(y, &y);
66 if (modf(0.5*y, &x1)) neg = -neg;
67 if (cos_flag) y -= 0.5;
68 x2 = modf(x, &x1);
69 x = x1 - y * A1;
70 x += x2;
71 x -= y * A2;
72#undef A1
73#undef A2
74 }
75
76 if (x < 0) {
77 neg = -neg;
78 x = -x;
79 }
80
81 /* ??? avoid underflow ??? */
82
83 y = x * x;
84 x += x * y * POLYNOM7(y, r);
85 return neg==-1 ? -x : x;
86}
87
88double
89sin(double x)
90{
91 return sinus(x, 0);
92}
93
94double
95cos(double x)
96{
97 if (x < 0) x = -x;
98 return sinus(x, 1);
99}
Note: See TracBrowser for help on using the repository browser.