|
|
1.1 ! root 1: /* @(#)sin.c 4.1 12/25/82 */ ! 2: ! 3: /* ! 4: C program for floating point sin/cos. ! 5: Calls modf. ! 6: There are no error exits. ! 7: Coefficients are #3370 from Hart & Cheney (18.80D). ! 8: */ ! 9: ! 10: static double twoopi = 0.63661977236758134308; ! 11: static double p0 = .1357884097877375669092680e8; ! 12: static double p1 = -.4942908100902844161158627e7; ! 13: static double p2 = .4401030535375266501944918e6; ! 14: static double p3 = -.1384727249982452873054457e5; ! 15: static double p4 = .1459688406665768722226959e3; ! 16: static double q0 = .8644558652922534429915149e7; ! 17: static double q1 = .4081792252343299749395779e6; ! 18: static double q2 = .9463096101538208180571257e4; ! 19: static double q3 = .1326534908786136358911494e3; ! 20: ! 21: double ! 22: cos(arg) ! 23: double arg; ! 24: { ! 25: double sinus(); ! 26: if(arg<0) ! 27: arg = -arg; ! 28: return(sinus(arg, 1)); ! 29: } ! 30: ! 31: double ! 32: sin(arg) ! 33: double arg; ! 34: { ! 35: double sinus(); ! 36: return(sinus(arg, 0)); ! 37: } ! 38: ! 39: static double ! 40: sinus(arg, quad) ! 41: double arg; ! 42: int quad; ! 43: { ! 44: double modf(); ! 45: double e, f; ! 46: double ysq; ! 47: double x,y; ! 48: int k; ! 49: double temp1, temp2; ! 50: ! 51: x = arg; ! 52: if(x<0) { ! 53: x = -x; ! 54: quad = quad + 2; ! 55: } ! 56: x = x*twoopi; /*underflow?*/ ! 57: if(x>32764){ ! 58: y = modf(x,&e); ! 59: e = e + quad; ! 60: modf(0.25*e,&f); ! 61: quad = e - 4*f; ! 62: }else{ ! 63: k = x; ! 64: y = x - k; ! 65: quad = (quad + k) & 03; ! 66: } ! 67: if (quad & 01) ! 68: y = 1-y; ! 69: if(quad > 1) ! 70: y = -y; ! 71: ! 72: ysq = y*y; ! 73: temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; ! 74: temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); ! 75: return(temp1/temp2); ! 76: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.