![[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.