← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cephes/sincos.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cephes/sincos.c
 
1
/*                                                      sincos.c
 
2
 *
 
3
 *      Circular sine and cosine of argument in degrees
 
4
 *      Table lookup and interpolation algorithm
 
5
 *
 
6
 *
 
7
 *
 
8
 * SYNOPSIS:
 
9
 *
 
10
 * double x, sine, cosine, flg, sincos();
 
11
 *
 
12
 * sincos( x, &sine, &cosine, flg );
 
13
 *
 
14
 *
 
15
 *
 
16
 * DESCRIPTION:
 
17
 *
 
18
 * Returns both the sine and the cosine of the argument x.
 
19
 * Several different compile time options and minimax
 
20
 * approximations are supplied to permit tailoring the
 
21
 * tradeoff between computation speed and accuracy.
 
22
 * 
 
23
 * Since range reduction is time consuming, the reduction
 
24
 * of x modulo 360 degrees is also made optional.
 
25
 *
 
26
 * sin(i) is internally tabulated for 0 <= i <= 90 degrees.
 
27
 * Approximation polynomials, ranging from linear interpolation
 
28
 * to cubics in (x-i)**2, compute the sine and cosine
 
29
 * of the residual x-i which is between -0.5 and +0.5 degree.
 
30
 * In the case of the high accuracy options, the residual
 
31
 * and the tabulated values are combined using the trigonometry
 
32
 * formulas for sin(A+B) and cos(A+B).
 
33
 *
 
34
 * Compile time options are supplied for 5, 11, or 17 decimal
 
35
 * relative accuracy (ACC5, ACC11, ACC17 respectively).
 
36
 * A subroutine flag argument "flg" chooses betwen this
 
37
 * accuracy and table lookup only (peak absolute error
 
38
 * = 0.0087).
 
39
 *
 
40
 * If the argument flg = 1, then the tabulated value is
 
41
 * returned for the nearest whole number of degrees. The
 
42
 * approximation polynomials are not computed.  At
 
43
 * x = 0.5 deg, the absolute error is then sin(0.5) = 0.0087.
 
44
 *
 
45
 * An intermediate speed and precision can be obtained using
 
46
 * the compile time option LINTERP and flg = 1.  This yields
 
47
 * a linear interpolation using a slope estimated from the sine
 
48
 * or cosine at the nearest integer argument.  The peak absolute
 
49
 * error with this option is 3.8e-5.  Relative error at small
 
50
 * angles is about 1e-5.
 
51
 *
 
52
 * If flg = 0, then the approximation polynomials are computed
 
53
 * and applied.
 
54
 *
 
55
 *
 
56
 *
 
57
 * SPEED:
 
58
 *
 
59
 * Relative speed comparisons follow for 6MHz IBM AT clone
 
60
 * and Microsoft C version 4.0.  These figures include
 
61
 * software overhead of do loop and function calls.
 
62
 * Since system hardware and software vary widely, the
 
63
 * numbers should be taken as representative only.
 
64
 *
 
65
 *                      flg=0   flg=0   flg=1   flg=1
 
66
 *                      ACC11   ACC5    LINTERP Lookup only
 
67
 * In-line 8087 (/FPi)
 
68
 * sin(), cos()         1.0     1.0     1.0     1.0
 
69
 *
 
70
 * In-line 8087 (/FPi)
 
71
 * sincos()             1.1     1.4     1.9     3.0
 
72
 *
 
73
 * Software (/FPa)
 
74
 * sin(), cos()         0.19    0.19    0.19    0.19
 
75
 *
 
76
 * Software (/FPa)
 
77
 * sincos()             0.39    0.50    0.73    1.7
 
78
 *
 
79
 *
 
80
 *
 
81
 * ACCURACY:
 
82
 *
 
83
 * The accurate approximations are designed with a relative error
 
84
 * criterion.  The absolute error is greatest at x = 0.5 degree.
 
85
 * It decreases from a local maximum at i+0.5 degrees to full
 
86
 * machine precision at each integer i degrees.  With the
 
87
 * ACC5 option, the relative error of 6.3e-6 is equivalent to
 
88
 * an absolute angular error of 0.01 arc second in the argument
 
89
 * at x = i+0.5 degrees.  For small angles < 0.5 deg, the ACC5
 
90
 * accuracy is 6.3e-6 (.00063%) of reading; i.e., the absolute
 
91
 * error decreases in proportion to the argument.  This is true
 
92
 * for both the sine and cosine approximations, since the latter
 
93
 * is for the function 1 - cos(x).
 
94
 *
 
95
 * If absolute error is of most concern, use the compile time
 
96
 * option ABSERR to obtain an absolute error of 2.7e-8 for ACC5
 
97
 * precision.  This is about half the absolute error of the
 
98
 * relative precision option.  In this case the relative error
 
99
 * for small angles will increase to 9.5e-6 -- a reasonable
 
100
 * tradeoff.
 
101
 */
 
102
 
 
103
 
 
104
#include "mconf.h"
 
105
 
 
106
/* Define one of the following to be 1:
 
107
 */
 
108
#define ACC5 1
 
109
#define ACC11 0
 
110
#define ACC17 0
 
111
 
 
112
/* Option for linear interpolation when flg = 1
 
113
 */
 
114
#define LINTERP 1
 
115
 
 
116
/* Option for absolute error criterion
 
117
 */
 
118
#define ABSERR 1
 
119
 
 
120
/* Option to include modulo 360 function:
 
121
 */
 
122
#define MOD360 0
 
123
 
 
124
/*
 
125
Cephes Math Library Release 2.1
 
126
Copyright 1987 by Stephen L. Moshier
 
127
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
128
*/
 
129
 
 
130
 
 
131
/* Table of sin(i degrees)
 
132
 * for 0 <= i <= 90
 
133
 */
 
134
static double sintbl[92] = {
 
135
  0.00000000000000000000E0,
 
136
  1.74524064372835128194E-2,
 
137
  3.48994967025009716460E-2,
 
138
  5.23359562429438327221E-2,
 
139
  6.97564737441253007760E-2,
 
140
  8.71557427476581735581E-2,
 
141
  1.04528463267653471400E-1,
 
142
  1.21869343405147481113E-1,
 
143
  1.39173100960065444112E-1,
 
144
  1.56434465040230869010E-1,
 
145
  1.73648177666930348852E-1,
 
146
  1.90808995376544812405E-1,
 
147
  2.07911690817759337102E-1,
 
148
  2.24951054343864998051E-1,
 
149
  2.41921895599667722560E-1,
 
150
  2.58819045102520762349E-1,
 
151
  2.75637355816999185650E-1,
 
152
  2.92371704722736728097E-1,
 
153
  3.09016994374947424102E-1,
 
154
  3.25568154457156668714E-1,
 
155
  3.42020143325668733044E-1,
 
156
  3.58367949545300273484E-1,
 
157
  3.74606593415912035415E-1,
 
158
  3.90731128489273755062E-1,
 
159
  4.06736643075800207754E-1,
 
160
  4.22618261740699436187E-1,
 
161
  4.38371146789077417453E-1,
 
162
  4.53990499739546791560E-1,
 
163
  4.69471562785890775959E-1,
 
164
  4.84809620246337029075E-1,
 
165
  5.00000000000000000000E-1,
 
166
  5.15038074910054210082E-1,
 
167
  5.29919264233204954047E-1,
 
168
  5.44639035015027082224E-1,
 
169
  5.59192903470746830160E-1,
 
170
  5.73576436351046096108E-1,
 
171
  5.87785252292473129169E-1,
 
172
  6.01815023152048279918E-1,
 
173
  6.15661475325658279669E-1,
 
174
  6.29320391049837452706E-1,
 
175
  6.42787609686539326323E-1,
 
176
  6.56059028990507284782E-1,
 
177
  6.69130606358858213826E-1,
 
178
  6.81998360062498500442E-1,
 
179
  6.94658370458997286656E-1,
 
180
  7.07106781186547524401E-1,
 
181
  7.19339800338651139356E-1,
 
182
  7.31353701619170483288E-1,
 
183
  7.43144825477394235015E-1,
 
184
  7.54709580222771997943E-1,
 
185
  7.66044443118978035202E-1,
 
186
  7.77145961456970879980E-1,
 
187
  7.88010753606721956694E-1,
 
188
  7.98635510047292846284E-1,
 
189
  8.09016994374947424102E-1,
 
190
  8.19152044288991789684E-1,
 
191
  8.29037572555041692006E-1,
 
192
  8.38670567945424029638E-1,
 
193
  8.48048096156425970386E-1,
 
194
  8.57167300702112287465E-1,
 
195
  8.66025403784438646764E-1,
 
196
  8.74619707139395800285E-1,
 
197
  8.82947592858926942032E-1,
 
198
  8.91006524188367862360E-1,
 
199
  8.98794046299166992782E-1,
 
200
  9.06307787036649963243E-1,
 
201
  9.13545457642600895502E-1,
 
202
  9.20504853452440327397E-1,
 
203
  9.27183854566787400806E-1,
 
204
  9.33580426497201748990E-1,
 
205
  9.39692620785908384054E-1,
 
206
  9.45518575599316810348E-1,
 
207
  9.51056516295153572116E-1,
 
208
  9.56304755963035481339E-1,
 
209
  9.61261695938318861916E-1,
 
210
  9.65925826289068286750E-1,
 
211
  9.70295726275996472306E-1,
 
212
  9.74370064785235228540E-1,
 
213
  9.78147600733805637929E-1,
 
214
  9.81627183447663953497E-1,
 
215
  9.84807753012208059367E-1,
 
216
  9.87688340595137726190E-1,
 
217
  9.90268068741570315084E-1,
 
218
  9.92546151641322034980E-1,
 
219
  9.94521895368273336923E-1,
 
220
  9.96194698091745532295E-1,
 
221
  9.97564050259824247613E-1,
 
222
  9.98629534754573873784E-1,
 
223
  9.99390827019095730006E-1,
 
224
  9.99847695156391239157E-1,
 
225
  1.00000000000000000000E0,
 
226
  9.99847695156391239157E-1,
 
227
};
 
228
 
 
229
#ifndef ANSIPROT
 
230
double floor();
 
231
#else
 
232
extern void sincos ( double x, double *s, double *c, int flg );
 
233
#endif
 
234
 
 
235
void
 
236
sincos(x, s, c, flg)
 
237
double x;
 
238
double *s, *c;
 
239
int flg;
 
240
{
 
241
int ix, ssign, csign, xsign;
 
242
double y, z, sx, sz, cx, cz;
 
243
 
 
244
/* Make argument nonnegative.
 
245
 */
 
246
xsign = 1;
 
247
if( x < 0.0 )
 
248
        {
 
249
        xsign = -1;
 
250
        x = -x;
 
251
        }
 
252
 
 
253
 
 
254
#if MOD360
 
255
x = x  -  360.0 * floor( x/360.0 );
 
256
#endif
 
257
 
 
258
/* Find nearest integer to x.
 
259
 * Note there should be a domain error test here,
 
260
 * but this is omitted to gain speed.
 
261
 */
 
262
ix = x + 0.5;
 
263
z = x - ix;             /* the residual */
 
264
 
 
265
/* Look up the sine and cosine of the integer.
 
266
 */
 
267
if( ix <= 180 )
 
268
        {
 
269
        ssign = 1;
 
270
        csign = 1;
 
271
        }
 
272
else
 
273
        {
 
274
        ssign = -1;
 
275
        csign = -1;
 
276
        ix -= 180;
 
277
        }
 
278
 
 
279
if( ix > 90 )
 
280
        {
 
281
        csign = -csign;
 
282
        ix = 180 - ix;
 
283
        }
 
284
 
 
285
sx = sintbl[ix];
 
286
if( ssign < 0 )
 
287
        sx = -sx;
 
288
cx = sintbl[ 90-ix ];
 
289
if( csign < 0 )
 
290
        cx = -cx;
 
291
 
 
292
/* If the flag argument is set, then just return
 
293
 * the tabulated values for arg to the nearest whole degree.
 
294
 */
 
295
if( flg )
 
296
        {
 
297
#if LINTERP
 
298
        y = sx + 1.74531263774940077459e-2 * z * cx;
 
299
        cx -= 1.74531263774940077459e-2 * z * sx;
 
300
        sx = y;
 
301
#endif
 
302
        if( xsign < 0 )
 
303
                sx = -sx;
 
304
        *s = sx;        /* sine */
 
305
        *c = cx;        /* cosine */
 
306
        return;
 
307
        }
 
308
 
 
309
 
 
310
if( ssign < 0 )
 
311
        sx = -sx;
 
312
if( csign < 0 )
 
313
        cx = -cx;
 
314
 
 
315
/* Find sine and cosine
 
316
 * of the residual angle between -0.5 and +0.5 degree.
 
317
 */
 
318
#if ACC5
 
319
#if ABSERR
 
320
/* absolute error = 2.769e-8: */
 
321
sz = 1.74531263774940077459e-2 * z;
 
322
/* absolute error = 4.146e-11: */
 
323
cz = 1.0 - 1.52307909153324666207e-4 * z * z;
 
324
#else
 
325
/* relative error = 6.346e-6: */
 
326
sz = 1.74531817576426662296e-2 * z;
 
327
/* relative error = 3.173e-6: */
 
328
cz = 1.0 - 1.52308226602566149927e-4 * z * z;
 
329
#endif
 
330
#else
 
331
y = z * z;
 
332
#endif
 
333
 
 
334
 
 
335
#if ACC11
 
336
sz = ( -8.86092781698004819918e-7 * y
 
337
      + 1.74532925198378577601e-2     ) * z;
 
338
 
 
339
cz = 1.0 - ( -3.86631403698859047896e-9 * y
 
340
            + 1.52308709893047593702e-4     ) * y;
 
341
#endif
 
342
 
 
343
 
 
344
#if ACC17
 
345
sz = ((  1.34959795251974073996e-11 * y
 
346
       - 8.86096155697856783296e-7     ) * y
 
347
       + 1.74532925199432957214e-2          ) * z;
 
348
 
 
349
cz = 1.0 - ((  3.92582397764340914444e-14 * y
 
350
             - 3.86632385155548605680e-9     ) * y
 
351
             + 1.52308709893354299569e-4          ) * y;
 
352
#endif
 
353
 
 
354
 
 
355
/* Combine the tabulated part and the calculated part
 
356
 * by trigonometry.
 
357
 */
 
358
y = sx * cz  +  cx * sz;
 
359
if( xsign < 0 )
 
360
        y = - y;
 
361
*s = y; /* sine */
 
362
 
 
363
*c = cx * cz  -  sx * sz; /* cosine */
 
364
}