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