trig.h revision 1.6 1 1.6 agc /* $NetBSD: trig.h,v 1.6 2003/08/07 16:44:53 agc Exp $ */
2 1.1 ragge /*
3 1.1 ragge * Copyright (c) 1987, 1993
4 1.1 ragge * The Regents of the University of California. All rights reserved.
5 1.1 ragge *
6 1.1 ragge * Redistribution and use in source and binary forms, with or without
7 1.1 ragge * modification, are permitted provided that the following conditions
8 1.1 ragge * are met:
9 1.1 ragge * 1. Redistributions of source code must retain the above copyright
10 1.1 ragge * notice, this list of conditions and the following disclaimer.
11 1.1 ragge * 2. Redistributions in binary form must reproduce the above copyright
12 1.1 ragge * notice, this list of conditions and the following disclaimer in the
13 1.1 ragge * documentation and/or other materials provided with the distribution.
14 1.6 agc * 3. Neither the name of the University nor the names of its contributors
15 1.1 ragge * may be used to endorse or promote products derived from this software
16 1.1 ragge * without specific prior written permission.
17 1.1 ragge *
18 1.1 ragge * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19 1.1 ragge * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 1.1 ragge * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 1.1 ragge * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22 1.1 ragge * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 1.1 ragge * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24 1.1 ragge * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25 1.1 ragge * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 1.1 ragge * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27 1.1 ragge * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28 1.1 ragge * SUCH DAMAGE.
29 1.1 ragge *
30 1.1 ragge * @(#)trig.h 8.1 (Berkeley) 6/4/93
31 1.1 ragge */
32 1.1 ragge
33 1.1 ragge vc(thresh, 2.6117239648121182150E-1 ,b863,3f85,6ea0,6b02, -1, .85B8636B026EA0)
34 1.1 ragge vc(PIo4, 7.8539816339744830676E-1 ,0fda,4049,68c2,a221, 0, .C90FDAA22168C2)
35 1.1 ragge vc(PIo2, 1.5707963267948966135E0 ,0fda,40c9,68c2,a221, 1, .C90FDAA22168C2)
36 1.1 ragge vc(PI3o4, 2.3561944901923449203E0 ,cbe3,4116,0e92,f999, 2, .96CBE3F9990E92)
37 1.1 ragge vc(PI, 3.1415926535897932270E0 ,0fda,4149,68c2,a221, 2, .C90FDAA22168C2)
38 1.1 ragge vc(PI2, 6.2831853071795864540E0 ,0fda,41c9,68c2,a221, 3, .C90FDAA22168C2)
39 1.1 ragge
40 1.1 ragge ic(thresh, 2.6117239648121182150E-1 , -2, 1.0B70C6D604DD4)
41 1.1 ragge ic(PIo4, 7.8539816339744827900E-1 , -1, 1.921FB54442D18)
42 1.1 ragge ic(PIo2, 1.5707963267948965580E0 , 0, 1.921FB54442D18)
43 1.1 ragge ic(PI3o4, 2.3561944901923448370E0 , 1, 1.2D97C7F3321D2)
44 1.1 ragge ic(PI, 3.1415926535897931160E0 , 1, 1.921FB54442D18)
45 1.1 ragge ic(PI2, 6.2831853071795862320E0 , 2, 1.921FB54442D18)
46 1.1 ragge
47 1.1 ragge #ifdef vccast
48 1.1 ragge #define thresh vccast(thresh)
49 1.1 ragge #define PIo4 vccast(PIo4)
50 1.1 ragge #define PIo2 vccast(PIo2)
51 1.1 ragge #define PI3o4 vccast(PI3o4)
52 1.1 ragge #define PI vccast(PI)
53 1.1 ragge #define PI2 vccast(PI2)
54 1.1 ragge #endif
55 1.1 ragge
56 1.1 ragge #ifdef national
57 1.1 ragge static long fmaxx[] = { 0xffffffff, 0x7fefffff};
58 1.1 ragge #define fmax (*(double*)fmaxx)
59 1.1 ragge #endif /* national */
60 1.1 ragge
61 1.5 matt #ifdef _LIBM_DECLARE
62 1.5 matt const double
63 1.5 matt __zero = 0,
64 1.5 matt __one = 1,
65 1.5 matt __negone = -1,
66 1.5 matt __half = 1.0/2.0,
67 1.5 matt #ifdef __vax__
68 1.5 matt __small = 1E-9, /* 1+small**2 == 1; better values for small:
69 1.5 matt * small = 1.5E-9 for VAX D
70 1.5 matt * = 1.2E-8 for IEEE Double
71 1.5 matt * = 2.8E-10 for IEEE Extended
72 1.5 matt */
73 1.5 matt __big = 1E18; /* big := 1/(small**2) */
74 1.5 matt #else
75 1.5 matt __small = 1E-10, /* 1+small**2 == 1; better values for small:
76 1.5 matt * small = 1.5E-9 for VAX D
77 1.5 matt * = 1.2E-8 for IEEE Double
78 1.5 matt * = 2.8E-10 for IEEE Extended
79 1.5 matt */
80 1.5 matt __big = 1E20; /* big := 1/(small**2) */
81 1.5 matt #endif
82 1.5 matt #else
83 1.5 matt extern const double __zero, __one, __negone, __half, __small, __big;
84 1.5 matt #endif
85 1.1 ragge
86 1.1 ragge /* sin__S(x*x) ... re-implemented as a macro
87 1.1 ragge * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
88 1.4 simonb * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X)
89 1.4 simonb * CODED IN C BY K.C. NG, 1/21/85;
90 1.1 ragge * REVISED BY K.C. NG on 8/13/85.
91 1.1 ragge *
92 1.1 ragge * sin(x*k) - x
93 1.1 ragge * RETURN --------------- on [-PI/4,PI/4] , where k=pi/PI, PI is the rounded
94 1.4 simonb * x
95 1.1 ragge * value of pi in machine precision:
96 1.1 ragge *
97 1.1 ragge * Decimal:
98 1.4 simonb * pi = 3.141592653589793 23846264338327 .....
99 1.1 ragge * 53 bits PI = 3.141592653589793 115997963 ..... ,
100 1.4 simonb * 56 bits PI = 3.141592653589793 227020265 ..... ,
101 1.1 ragge *
102 1.1 ragge * Hexadecimal:
103 1.1 ragge * pi = 3.243F6A8885A308D313198A2E....
104 1.1 ragge * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18
105 1.4 simonb * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2
106 1.1 ragge *
107 1.1 ragge * Method:
108 1.4 simonb * 1. Let z=x*x. Create a polynomial approximation to
109 1.1 ragge * (sin(k*x)-x)/x = z*(S0 + S1*z^1 + ... + S5*z^5).
110 1.1 ragge * Then
111 1.1 ragge * sin__S(x*x) = z*(S0 + S1*z^1 + ... + S5*z^5)
112 1.1 ragge *
113 1.1 ragge * The coefficient S's are obtained by a special Remez algorithm.
114 1.1 ragge *
115 1.1 ragge * Accuracy:
116 1.4 simonb * In the absence of rounding error, the approximation has absolute error
117 1.4 simonb * less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE.
118 1.1 ragge *
119 1.1 ragge * Constants:
120 1.1 ragge * The hexadecimal values are the intended ones for the following constants.
121 1.1 ragge * The decimal values may be used, provided that the compiler will convert
122 1.1 ragge * from decimal to binary accurately enough to produce the hexadecimal values
123 1.1 ragge * shown.
124 1.1 ragge *
125 1.1 ragge */
126 1.1 ragge
127 1.1 ragge vc(S0, -1.6666666666666646660E-1 ,aaaa,bf2a,aa71,aaaa, -2, -.AAAAAAAAAAAA71)
128 1.1 ragge vc(S1, 8.3333333333297230413E-3 ,8888,3d08,477f,8888, -6, .8888888888477F)
129 1.1 ragge vc(S2, -1.9841269838362403710E-4 ,0d00,ba50,1057,cf8a, -12, -.D00D00CF8A1057)
130 1.1 ragge vc(S3, 2.7557318019967078930E-6 ,ef1c,3738,bedc,a326, -18, .B8EF1CA326BEDC)
131 1.1 ragge vc(S4, -2.5051841873876551398E-8 ,3195,b3d7,e1d3,374c, -25, -.D73195374CE1D3)
132 1.1 ragge vc(S5, 1.6028995389845827653E-10 ,3d9c,3030,cccc,6d26, -32, .B03D9C6D26CCCC)
133 1.1 ragge vc(S6, -6.2723499671769283121E-13 ,8d0b,ac30,ea82,7561, -40, -.B08D0B7561EA82)
134 1.1 ragge
135 1.1 ragge ic(S0, -1.6666666666666463126E-1 , -3, -1.555555555550C)
136 1.1 ragge ic(S1, 8.3333333332992771264E-3 , -7, 1.111111110C461)
137 1.1 ragge ic(S2, -1.9841269816180999116E-4 , -13, -1.A01A019746345)
138 1.1 ragge ic(S3, 2.7557309793219876880E-6 , -19, 1.71DE3209CDCD9)
139 1.1 ragge ic(S4, -2.5050225177523807003E-8 , -26, -1.AE5C0E319A4EF)
140 1.1 ragge ic(S5, 1.5868926979889205164E-10 , -33, 1.5CF61DF672B13)
141 1.1 ragge
142 1.1 ragge #ifdef vccast
143 1.1 ragge #define S0 vccast(S0)
144 1.1 ragge #define S1 vccast(S1)
145 1.1 ragge #define S2 vccast(S2)
146 1.1 ragge #define S3 vccast(S3)
147 1.1 ragge #define S4 vccast(S4)
148 1.1 ragge #define S5 vccast(S5)
149 1.1 ragge #define S6 vccast(S6)
150 1.1 ragge #endif
151 1.1 ragge
152 1.2 matt #if defined(__vax__)||defined(tahoe)
153 1.1 ragge # define sin__S(z) (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*(S5+z*S6)))))))
154 1.2 matt #else /* defined(__vax__)||defined(tahoe) */
155 1.1 ragge # define sin__S(z) (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))))
156 1.2 matt #endif /* defined(__vax__)||defined(tahoe) */
157 1.1 ragge
158 1.1 ragge /* cos__C(x*x) ... re-implemented as a macro
159 1.1 ragge * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
160 1.4 simonb * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X)
161 1.4 simonb * CODED IN C BY K.C. NG, 1/21/85;
162 1.1 ragge * REVISED BY K.C. NG on 8/13/85.
163 1.1 ragge *
164 1.4 simonb * x*x
165 1.1 ragge * RETURN cos(k*x) - 1 + ----- on [-PI/4,PI/4], where k = pi/PI,
166 1.4 simonb * 2
167 1.1 ragge * PI is the rounded value of pi in machine precision :
168 1.1 ragge *
169 1.1 ragge * Decimal:
170 1.4 simonb * pi = 3.141592653589793 23846264338327 .....
171 1.1 ragge * 53 bits PI = 3.141592653589793 115997963 ..... ,
172 1.4 simonb * 56 bits PI = 3.141592653589793 227020265 ..... ,
173 1.1 ragge *
174 1.1 ragge * Hexadecimal:
175 1.1 ragge * pi = 3.243F6A8885A308D313198A2E....
176 1.1 ragge * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18
177 1.4 simonb * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2
178 1.1 ragge *
179 1.1 ragge *
180 1.1 ragge * Method:
181 1.4 simonb * 1. Let z=x*x. Create a polynomial approximation to
182 1.1 ragge * cos(k*x)-1+z/2 = z*z*(C0 + C1*z^1 + ... + C5*z^5)
183 1.1 ragge * then
184 1.1 ragge * cos__C(z) = z*z*(C0 + C1*z^1 + ... + C5*z^5)
185 1.1 ragge *
186 1.1 ragge * The coefficient C's are obtained by a special Remez algorithm.
187 1.1 ragge *
188 1.1 ragge * Accuracy:
189 1.4 simonb * In the absence of rounding error, the approximation has absolute error
190 1.4 simonb * less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE.
191 1.4 simonb *
192 1.1 ragge *
193 1.1 ragge * Constants:
194 1.1 ragge * The hexadecimal values are the intended ones for the following constants.
195 1.1 ragge * The decimal values may be used, provided that the compiler will convert
196 1.1 ragge * from decimal to binary accurately enough to produce the hexadecimal values
197 1.1 ragge * shown.
198 1.1 ragge */
199 1.1 ragge
200 1.1 ragge vc(C0, 4.1666666666666504759E-2 ,aaaa,3e2a,a9f0,aaaa, -4, .AAAAAAAAAAA9F0)
201 1.1 ragge vc(C1, -1.3888888888865302059E-3 ,0b60,bbb6,0cca,b60a, -9, -.B60B60B60A0CCA)
202 1.1 ragge vc(C2, 2.4801587285601038265E-5 ,0d00,38d0,098f,cdcd, -15, .D00D00CDCD098F)
203 1.1 ragge vc(C3, -2.7557313470902390219E-7 ,f27b,b593,e805,b593, -21, -.93F27BB593E805)
204 1.1 ragge vc(C4, 2.0875623401082232009E-9 ,74c8,320f,3ff0,fa1e, -28, .8F74C8FA1E3FF0)
205 1.1 ragge vc(C5, -1.1355178117642986178E-11 ,c32d,ae47,5a63,0a5c, -36, -.C7C32D0A5C5A63)
206 1.1 ragge
207 1.1 ragge ic(C0, 4.1666666666666504759E-2 , -5, 1.555555555553E)
208 1.1 ragge ic(C1, -1.3888888888865301516E-3 , -10, -1.6C16C16C14199)
209 1.1 ragge ic(C2, 2.4801587269650015769E-5 , -16, 1.A01A01971CAEB)
210 1.1 ragge ic(C3, -2.7557304623183959811E-7 , -22, -1.27E4F1314AD1A)
211 1.1 ragge ic(C4, 2.0873958177697780076E-9 , -29, 1.1EE3B60DDDC8C)
212 1.1 ragge ic(C5, -1.1250289076471311557E-11 , -37, -1.8BD5986B2A52E)
213 1.1 ragge
214 1.1 ragge #ifdef vccast
215 1.1 ragge #define C0 vccast(C0)
216 1.1 ragge #define C1 vccast(C1)
217 1.1 ragge #define C2 vccast(C2)
218 1.1 ragge #define C3 vccast(C3)
219 1.1 ragge #define C4 vccast(C4)
220 1.1 ragge #define C5 vccast(C5)
221 1.1 ragge #endif
222 1.1 ragge
223 1.1 ragge #define cos__C(z) (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))))
224