![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to sind.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDTahoe / cci / libM|
Return to sind.c CVS log | Up to [CSRG BSD Unix] / 43BSDTahoe / cci / libM |
1.1 root 1: /* @(#)sin.c 4.1/4.2 10/31/84 CCI-CPG */
2:
3:
4: /*
5: * Double-precision floating point sin/cos function.
6: * Calls modf.
7: * Pn coefficients are #3345 from Hart & etc.(18.79D).
8: *
9: * New version by Les Powers (3/26/85).
10: */
11:
12:
13: static double twoopi = 0.63661977236758134308;
14: static double piotwo = 3.14159265358979323846264338327950288419716939937511/2.;
15: static double piotwols = .572118872610983179676266255859012e-17;
16: static double piofour = 3.14159265358979323846264338327950288419716939937511/4.;
17: static double threepio4 = 3.14159265358979323846264338327950288419716939937511*.75;
18:
19:
20: static double p0 = .15707963267948966188272e1;
21: static double p1 = -.64596409750624619108547e0;
22: static double p2 = .7969262624616543562977e-1;
23: static double p3 = -.468175413530264260121e-2;
24: static double p4 = .1604411847068220716e-3;
25: static double p5 = -.359884300720869272e-5;
26: static double p6 = .5692134872719023e-7;
27: static double p7 = -.66843217206396e-9;
28: static double p8 = .587061098171e-11;
29:
30: /*
31: * Sn coefficients are #3043 from Hart & etc.(17.48D).
32: * Cn coefficients are #3824 from Hart & etc.(19.45D).
33: * The coefficients are multiplied by the appropriate power
34: * of 4/pi using binary calculator with 40 digit arithmetic.
35: * The following statements in comments indicate the
36: * operations originally used to generate the coefficients.
37: *
38: * #define FP (4/3.14159265358979323846264338327950288419716939937511)
39: * #define FS FP*FP
40: *
41: * static double s0 = .785398163397448307014e0*FP;
42: * static double s1 = -.80745512188280530192e-1*FP*FS;
43: * static double s2 = .2490394570188736117e-2*FP*FS*FS;
44: * static double s3 = -.36576204158455695e-4*FP*FS*FS*FS;
45: * static double s4 = .313361621661904e-6*FP*FS*FS*FS*FS;
46: * static double s5 = -.1757149292755e-8*FP*FS*FS*FS*FS*FS;
47: * static double s6 = .6877100349e-11*FP*FS*FS*FS*FS*FS*FS;
48: *
49: * static double c0 = .99999999999999999996415e0;
50: * static double c1 = -.30842513753404245242414e0*FS;
51: * static double c2 = .1585434424381541089754e-1*FS*FS;
52: * static double c3 = -.32599188692668755044e-3*FS*FS*FS;
53: * static double c4 = .359086044588581953e-5*FS*FS*FS*FS;
54: * static double c5 = -.2461136382637005e-7*FS*FS*FS*FS*FS;
55: * static double c6 = .11500497024263e-9*FS*FS*FS*FS*FS*FS;
56: * static double c7 = -.38577620372e-12*FS*FS*FS*FS*FS*FS*FS;
57: */
58:
59: static double s0 = .9999999999999999966874625291130031549274;
60: static double s1 = -.1666666666666661475150697663344387409506;
61: static double s2 = .0083333333333200019685518812023464185712;
62: static double s3 = -.0001984126982840175394529286095367983768;
63: static double s4 = .0000027557313298885682473645051747116530;
64: static double s5 = -.0000000250507058263681718011320436448899;
65: static double s6 = .0000000001589413523004632819499530093565;
66:
67: static double c0 = 0.99999999999999999996415;
68: static double c1 = -.4999999999999999928435829208978050890709;
69: static double c2 = .0416666666666664302573692446884954307432;
70: static double c3 = -.0013888888888858960434395194724966502603;
71: static double c4 = .0000248015872828994630247806807330993132;
72: static double c5 = -.0000002755731286569608222434728722627332;
73: static double c6 = .0000000020875555145677882874779379760684;
74: static double c7 = -.0000000000113521232057839395845871741059;
75:
76:
77: double
78: cos(arg)
79: double arg;
80: {
81: double sinus();
82: double s;
83: union {
84: double d;
85: unsigned i;
86: } u;
87: if ( arg < 0 )
88: arg = -arg;
89: u.d = arg;
90: if ( u.i < 0x40490fdb ) { /* arg < pi * 1/4 */
91: s = arg * arg;
92: return(c0+s*(c1+s*(c2+s*(c3+s*(c4+s*(c5+s*(c6+s*c7)))))));
93: }
94: if ( u.i < 0x4116cbe3 ) { /* arg < pi * 3/4 */
95: arg = (piotwo - arg)+piotwols;
96: s = arg * arg;
97: return(arg*(s0+s*(s1+s*(s2+s*(s3+s*(s4+s*(s5+s*s6)))))));
98: }
99: return(sinus(arg, 1));
100: }
101:
102: double
103: sin(arg)
104: double arg;
105: {
106: double sinus();
107: double s;
108: union {
109: double d;
110: unsigned i;
111: } u;
112: if ( arg < 0 )
113: return( -sin(-arg) );
114: u.d = arg;
115: if ( u.i < 0x40490fdb ) { /* arg < pi * 1/4 */
116: s = arg * arg;
117: return(arg*(s0+s*(s1+s*(s2+s*(s3+s*(s4+s*(s5+s*s6)))))));
118: }
119: if ( u.i < 0x4116cbe3 ) { /* arg < pi * 3/4 */
120: arg = (piotwo - arg)+piotwols;
121: s = arg * arg;
122: return(c0+s*(c1+s*(c2+s*(c3+s*(c4+s*(c5+s*(c6+s*c7)))))));
123: }
124: return(sinus(arg, 0));
125: }
126:
127: static double
128: sinus(arg, quad)
129: double arg;
130: int quad;
131: {
132: double modf();
133: double e, f;
134: double ysq;
135: double x,y;
136: int k;
137: double temp1, temp2;
138:
139: x = arg;
140: if(x < 0) {
141: x = -x;
142: quad = quad + 2;
143: }
144: x = x * twoopi; /*underflow?*/
145: if (x > 32764){
146: y = modf(x,&e);
147: e = e + quad;
148: modf(0.25 * e,&f);
149: quad = e - 4 * f;
150: }else{
151: k = x;
152: y = x - k;
153: quad = (quad + k) & 03;
154: }
155: if (quad & 01)
156: y = 1-y;
157: if (quad > 1)
158: y = -y;
159:
160: ysq = y * y;
161: temp1 = ((((((((p8*ysq+p7)*ysq+p6)*ysq+p5)*ysq+p4)*ysq+
162: p3)*ysq+p2)*ysq+p1)*ysq+p0);
163: return(temp1 * y);
164: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.