![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to math.icn CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10dc / cmd / icon / ipl / source|
Return to math.icn CVS log | Up to [Research Unix] / researchv10dc / cmd / icon / ipl / source |
1.1 ! root 1: # MATH ! 2: # ! 3: # Mathematical procedures ! 4: # ! 5: # George D. Yee ! 6: # 1847 N. Frances Blvd. ! 7: # Tucson, AZ 85712 ! 8: # ! 9: # Last modified 5/6/86 by Ralph E. Griswold ! 10: # ! 11: # ! 12: # Free distribution and use of this material is granted provided the ! 13: # above credit is left intact on all source copies. No warranties ! 14: # are made as to the correctness or suitability of these procedures ! 15: # for any purpose. Please send any suggestions to the author at the above ! 16: # address. ! 17: # ! 18: ! 19: procedure sin(x) ! 20: return _sinus(numeric(x),0) ! 21: end ! 22: ! 23: procedure cos(x) ! 24: return _sinus(abs(numeric(x)),1) ! 25: end ! 26: ! 27: procedure tan(x) ! 28: return sin(x) / (0.0 ~= cos(x)) ! 29: end ! 30: ! 31: # atan returns the value of the arctangent of its ! 32: # argument in the range [-pi/2,pi/2]. ! 33: procedure atan(x) ! 34: if numeric(x) then ! 35: return if x > 0.0 then _satan(x) else -_satan(-x) ! 36: end ! 37: ! 38: # atan2 returns the arctangent of y/x ! 39: # in the range [-pi,pi]. ! 40: procedure atan2(y,x) ! 41: local r ! 42: static pi ! 43: initial pi := 3.141592653589793238462643 ! 44: return if numeric(y) & numeric(x) then { ! 45: if x > 0.0 then ! 46: atan(y/x) ! 47: else if x < 0.0 then { ! 48: r := pi - atan(abs(y/x)) ! 49: if y >= 0.0 then r else -r ! 50: } ! 51: else if x = y = 0.0 then ! 52: 0.0 # special value if both x and y are zero ! 53: else ! 54: if y >= 0.0 then pi/2.0 else -pi/2.0 ! 55: } ! 56: end ! 57: ! 58: procedure asin(x) ! 59: if abs(numeric(x)) <= 1.0 then ! 60: return atan2(x, (1.0-(x^2))^0.5) ! 61: end ! 62: ! 63: procedure acos(x) ! 64: return 1.570796326794896619231e0 - asin(x) ! 65: end ! 66: ! 67: procedure dtor(deg) ! 68: return numeric(deg)/57.29577951308232 ! 69: end ! 70: ! 71: procedure rtod(rad) ! 72: return numeric(rad)*57.29577951308232 ! 73: end ! 74: ! 75: procedure sqrt(x) ! 76: return (0.0 <= numeric(x)) ^ 0.5 ! 77: end ! 78: ! 79: procedure floor(x) ! 80: return if numeric(x) then ! 81: if x>=0.0 | real(x)=integer(x) then integer(x) else -integer(-x+1) ! 82: end ! 83: ! 84: procedure ceil(x) ! 85: return -floor(-numeric(x)) ! 86: end ! 87: ! 88: procedure log(x) ! 89: local z, zsq, ex ! 90: static log2, sqrto2, p0, p1, p2, p3, q0, q1, q2 ! 91: initial { ! 92: # The coefficients are #2705 from Hart & Cheney. (19.38D) ! 93: log2 := 0.693147180559945309e0 ! 94: sqrto2 := 0.707106781186547524e0 ! 95: p0 := -0.240139179559210510e2 ! 96: p1 := 0.309572928215376501e2 ! 97: p2 := -0.963769093368686593e1 ! 98: p3 := 0.421087371217979714e0 ! 99: q0 := -0.120069589779605255e2 ! 100: q1 := 0.194809660700889731e2 ! 101: q2 := -0.891110902798312337e1 ! 102: } ! 103: if numeric(x) > 0.0 then { ! 104: ex := 0 ! 105: while x >= 1.0 do { ! 106: x /:= 2.0 ! 107: ex +:= 1 ! 108: } ! 109: while x < 0.5 do { ! 110: x *:= 2.0 ! 111: ex -:= 1 ! 112: } ! 113: if x < sqrto2 then { ! 114: x *:= 2.0 ! 115: ex -:= 1 ! 116: } ! 117: return ((((p3*(zsq:=(z:=(x-1.0)/(x+1.0))^2)+p2)*zsq+p1)*zsq+p0)/ ! 118: (((1.0*zsq+q2)*zsq+q1)*zsq+q0))*z+ex*log2 ! 119: } ! 120: end ! 121: ! 122: procedure exp(x) ! 123: return 2.718281828459045235360287 ^ numeric(x) ! 124: end ! 125: ! 126: procedure log10(x) ! 127: return log(x)/2.30258509299404568402 ! 128: end ! 129: ! 130: procedure _sinus(x,quad) ! 131: local ysq, y, k ! 132: static twoopi, p0, p1, p2, p3, p4, q0, q1, q2, q3 ! 133: initial { ! 134: # Coefficients are #3370 from Hart & Cheney (18.80D). ! 135: twoopi := 0.63661977236758134308 ! 136: p0 := 0.1357884097877375669092680e8 ! 137: p1 := -0.4942908100902844161158627e7 ! 138: p2 := 0.4401030535375266501944918e6 ! 139: p3 := -0.1384727249982452873054457e5 ! 140: p4 := 0.1459688406665768722226959e3 ! 141: q0 := 0.8644558652922534429915149e7 ! 142: q1 := 0.4081792252343299749395779e6 ! 143: q2 := 0.9463096101538208180571257e4 ! 144: q3 := 0.1326534908786136358911494e3 ! 145: } ! 146: if x < 0.0 then { ! 147: x := -x ! 148: quad +:= 2 ! 149: } ! 150: y := (x *:= twoopi) - (k := integer(x)) ! 151: if (quad := (quad + k) % 4) = (1|3) then ! 152: y := 1.0 - y ! 153: if quad > 1 then ! 154: y := -y ! 155: return (((((p4*(ysq:=y^2)+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y) / ! 156: ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0) ! 157: end ! 158: ! 159: procedure _satan(x) ! 160: static sq2p1,sq2m1,pio2,pio4 ! 161: initial { ! 162: sq2p1 := 2.414213562373095048802e0 ! 163: sq2m1 := 0.414213562373095048802e0 ! 164: pio2 := 1.570796326794896619231e0 ! 165: pio4 := 0.785398163397448309615e0 ! 166: } ! 167: return if x < sq2m1 then ! 168: _xatan(x) ! 169: else if x > sq2p1 then ! 170: pio2 - _xatan(1.0/x) ! 171: else ! 172: pio4 + _xatan((x-1.0)/(x+1.0)) ! 173: end ! 174: ! 175: procedure _xatan(x) ! 176: local xsq ! 177: static p4,p3,p2,p1,p0,q4,q3,q2,q1,q0 ! 178: initial { ! 179: # coefficients are #5077 from Hart & Cheney. (19.56D) ! 180: p4 := 0.161536412982230228262e2 ! 181: p3 := 0.26842548195503973794141e3 ! 182: p2 := 0.11530293515404850115428136e4 ! 183: p1 := 0.178040631643319697105464587e4 ! 184: p0 := 0.89678597403663861959987488e3 ! 185: q4 := 0.5895697050844462222791e2 ! 186: q3 := 0.536265374031215315104235e3 ! 187: q2 := 0.16667838148816337184521798e4 ! 188: q1 := 0.207933497444540981287275926e4 ! 189: q0 := 0.89678597403663861962481162e3 ! 190: } ! 191: return x * ((((p4*(xsq:=x^2)+p3)*xsq+p2)*xsq+p1)*xsq+p0) / ! 192: (((((xsq+q4)*xsq+q3)*xsq+q2)*xsq+q1)*xsq+q0) ! 193: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.