![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to math.sml CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / sml / src / boot|
Return to math.sml CVS log | Up to [Research Unix] / researchv10no / cmd / sml / src / boot |
1.1 root 1: (* Copyright 1989 by AT&T Bell Laboratories *)
2: (* The following functions were adapted from the 4.3BSD math library.
3: Eventually, each machine supported should have a hand-coded math
4: functor with more efficient versions of these functions.
5:
6: ***************************************************************************
7: * *
8: * Copyright (c) 1985 Regents of the University of California. *
9: * *
10: * Use and reproduction of this software are granted in accordance with *
11: * the terms and conditions specified in the Berkeley Software License *
12: * Agreement (in particular, this entails acknowledgement of the programs' *
13: * source, and inclusion of this notice) with the additional understanding *
14: * that all recipients should regard themselves as participants in an *
15: * ongoing research project and hence should feel obligated to report *
16: * their experiences (good or bad) with these elementary function codes, *
17: * using "sendbug 4bsd-bugs@BERKELEY", to the authors. *
18: * *
19: * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. *
20: * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. *
21: * *
22: ***************************************************************************
23:
24: *)
25:
26: signature MATH =
27: sig
28: exception Sqrt and Exp and Ln
29: val sqrt : real -> real
30: val sin : real -> real
31: val cos : real -> real
32: val arctan : real -> real
33: val exp : real -> real
34: val ln : real -> real
35: end
36:
37: structure Math : MATH = struct
38:
39: structure Assembly : ASSEMBLY = Core.Assembly
40:
41: infix 7 * /
42: infix 6 + -
43: infix 4 > < >= <=
44:
45: exception Real = Assembly.Real
46: exception Sqrt
47: exception Exp
48: exception Ln
49:
50: val scalb = Assembly.A.scalb
51: val logb = Assembly.A.logb
52:
53: val ~ = InLine.fneg
54: val op + = InLine.fadd
55: val op - = InLine.fsub
56: val op * = InLine.fmul
57: val op / = InLine.fdiv
58: val op > = InLine.fgt
59: val op < = InLine.flt
60: val op >= = InLine.fge
61: val op <= = InLine.fle
62:
63: val ieql = InLine.ieql
64: val ineq = InLine.ineq
65: val feql = InLine.feql
66:
67: val negone = ~1.0
68: val zero = 0.0
69: val half = 0.5
70: val one = 1.0
71: val two = 2.0
72: val four = 4.0
73:
74: fun copysign(a,b) =
75: case (a<zero,b<zero) of
76: (true,true) => a
77: | (false,false) => a
78: | _ => ~a
79:
80: fun abs x = if x < zero then ~x else x
81: fun mod(a,b) = InLine.-(a,InLine.*(InLine.div(a,b),b))
82: local fun rtoi (pred,init) =
83: let fun loop n =
84: if pred n
85: then init
86: else let val (i,r) = loop(n * half)
87: val i' = InLine.lshift(i,1)
88: val r' = r+r
89: in if n-r' < one then (i',r') else (InLine.+(i',1),r'+one)
90: end
91: in loop
92: end
93: val pton = rtoi(fn x => x < one, (0,zero))
94: val nton = rtoi(fn x => x >= negone, (~1,negone))
95: fun fl n = if n < zero then nton n else pton n
96: fun tr n = #1(pton n)
97: in
98: fun truncate n = if n < zero then InLine.~(tr(~n)) else tr n
99: val floor = Assembly.A.floor
100: fun realfloor n = #2(fl n)
101: fun ceiling n = InLine.~(floor(~n))
102: end
103: local fun loop 0 = zero
104: | loop n = let val x = two * loop(InLine.rshift(n,1))
105: in case InLine.andb(n,1) of 0 => x | 1 => x + one
106: end
107: in fun real n = if InLine.<(n,0) then ~(loop(InLine.~ n)) else loop n
108: end
109:
110: (* Not a true drem - rounding on .5 should go to the even case. *)
111: fun drem(x,y) =
112: let val n = x/y+half
113: val N = if abs n < 9.007199254741E15 (* about 2^53 *)
114: then realfloor n
115: handle Real _ =>
116: scalb(realfloor(scalb(n,InLine.~(logb n))),logb n)
117: else n
118: in x - N * y
119: end
120:
121: (* sin/cos *)
122: local
123: val S0 = ~1.6666666666666463126E~1
124: val S1 = 8.3333333332992771264E~3
125: val S2 = ~1.9841269816180999116E~4
126: val S3 = 2.7557309793219876880E~6
127: val S4 = ~2.5050225177523807003E~8
128: val S5 = 1.5868926979889205164E~10
129: in fun sin__S z = (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))))
130: end
131:
132: local
133: val C0 = 4.1666666666666504759E~2
134: val C1 = ~1.3888888888865301516E~3
135: val C2 = 2.4801587269650015769E~5
136: val C3 = ~2.7557304623183959811E~7
137: val C4 = 2.0873958177697780076E~9
138: val C5 = ~1.1250289076471311557E~11
139: in fun cos__C z = (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))))
140: end
141:
142: val PIo4 = 7.8539816339744827900E~1
143: val PIo2 = 1.5707963267948965580E0
144: val PI3o4 = 2.3561944901923448370E0
145: val PI = 3.1415926535897931160E0
146: val PI2 = 6.2831853071795862320E0
147:
148: local
149: val thresh = 2.6117239648121182150E~1
150: in fun S y = y + y * sin__S(y*y)
151: fun C y =
152: let val yy = y*y
153: val c = cos__C yy
154: val Y = yy/two
155: in if Y < thresh then one - (Y - c)
156: else half - (Y - half - c)
157: end
158: end
159: fun sin x =
160: let val x = drem(x,PI2)
161: val a = abs x
162: val y = if a >= PI3o4 then copysign(PI-a,x) else a
163: in if y >= PIo4 then copysign(C(PIo2-a),x)
164: else S y
165: end
166:
167: fun cos x =
168: let val x = drem(x,PI2)
169: val a = abs x
170: val (s,y) = if a >= PI3o4
171: then (fn x => ~x,PI-a)
172: else (fn x => x,a)
173: in if y < PIo4 then s(C y)
174: else S(PIo2-a)
175: end
176:
177: local
178: val p1 = 1.3887401997267371720E~2
179: val p2 = 3.3044019718331897649E~5
180: val q1 = 1.1110813732786649355E~1
181: val q2 = 9.9176615021572857300E~4
182: in fun exp__E(x:real,c:real) =
183: let val small=1.0E~19
184: val x = abs x
185: val z = x*x
186: val p = z*(p1+z*p2)
187: val q = z*(q1+z*q2)
188: val xp= x*p
189: val xh= x*half
190: val w = xh-(q-xp)
191: val p = p+p
192: val c = c+x*((xh*w-(q-(p+xp)))/(one-w)+c)
193: in if(copysign(x,one)>small) then z*half+c
194: else copysign(zero,x) (* could be inexact *)
195: end
196: end
197:
198: (* for exp and ln *)
199: val ln2hi = 6.9314718036912381649E~1
200: val ln2lo = 1.9082149292705877000E~10
201: val sqrt2 = 1.4142135623730951455E0
202: val lnhuge = 7.1602103751842355450E2
203: val lntiny = ~7.5137154372698068983E2
204: val invln2 = 1.4426950408889633870E0
205:
206: fun exp(x:real) =
207: let fun exp_norm() =
208: let (* argument reduction : x --> x - k*ln2 *)
209: val k = floor(invln2*x+copysign(half,x)) (* k=NINT(x/ln2) *)
210: val K = real k
211: (* express x-k*ln2 as z+c *)
212: val hi = x-K*ln2hi
213: val lo = K*ln2lo
214: val z = hi - lo
215: val c = (hi-z)-lo
216: (* return 2^k*[expm1(x) + 1] *)
217: val z = z + exp__E(z,c)
218: in scalb(z+one,k)
219: end
220: in if x <= lnhuge then if x >= lntiny
221: then exp_norm() handle Real _ => raise Exp
222: else zero (* exp(-big#) underflows to zero *)
223: (* exp(-INF) is zero *)
224: else raise Exp (* exp(INF) is INF, exp(+big#) overflows to INF *)
225: end
226:
227: local
228: val L1 = 6.6666666666667340202E~1
229: val L2 = 3.9999999999416702146E~1
230: val L3 = 2.8571428742008753154E~1
231: val L4 = 2.2222198607186277597E~1
232: val L5 = 1.8183562745289935658E~1
233: val L6 = 1.5314087275331442206E~1
234: val L7 = 1.4795612545334174692E~1
235: in fun log__L(z) = z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))
236: end
237:
238: fun ln(x:real) =
239: let fun findln() =
240: let (* argument reduction *)
241: val k = logb(x)
242: val x = scalb(x,InLine.~ k)
243: val (x,k) = if ieql(k,~1022) (* subnormal no. *)
244: then let val n = logb(x)
245: in (scalb(x,InLine.~ n),InLine.+(k,n)) end
246: else (x,k)
247: val (k,x) = if x >= sqrt2 then (InLine.+(k,1),x*half) else (k,x)
248: val K = real k
249: val x = x - one
250: (* compute log(1+x) *)
251: val s = x/(two+x)
252: val t = x*x*half
253: val z = K*ln2lo+s*(t+log__L(s*s))
254: val x = x + (z - t)
255: in K*ln2hi+x end
256: in if x <= zero then raise Ln
257: else findln() handle Real _ => raise Ln
258: end
259:
260: local
261: val small=1.0E~9
262: val big=1.0E18
263: val athfhi = 4.6364760900080609352E~1
264: val athflo = 4.6249969567426939759E~18
265: val at1fhi = 9.8279372324732905408E~1
266: val at1flo = ~2.4407677060164810007E~17
267: val a1 = 3.3333333333333942106E~1
268: val a2 = ~1.9999999999979536924E~1
269: val a3 = 1.4285714278004377209E~1
270: val a4 = ~1.1111110579344973814E~1
271: val a5 = 9.0908906105474668324E~2
272: val a6 = ~7.6919217767468239799E~2
273: val a7 = 6.6614695906082474486E~2
274: val a8 = ~5.8358371008508623523E~2
275: val a9 = 4.9850617156082015213E~2
276: val a10 = ~3.6700606902093604877E~2
277: val a11 = 1.6438029044759730479E~2
278: in
279: fun atan2(y:real,x:real) =
280: let (* copy down the sign of y and x *)
281: val signy = copysign(one,y)
282: val signx = copysign(one,x)
283: exception Done of real
284: fun findatan2(x,y,t) =
285: let (* begin argument reduction *)
286: val k = floor(four * (t+0.0625))
287: val (hi,lo,t) =
288: if t < 2.4375 then
289: (* truncate 4(t+1/16) to integer for branching *)
290: case k of
291: (* t is in [0,7/16] *)
292: 0 => if t < small
293: then raise Done
294: (copysign(if signx>zero then t
295: else PI-t,
296: signy))
297: else (zero,zero,t)
298: | 1 => if t < small
299: then raise Done
300: (copysign(if signx>zero then t
301: else PI-t,
302: signy))
303: else (zero,zero,t)
304: (* t is in [7/16,11/16] *)
305: | 2 => (athfhi,athflo,((y+y)-x)/(x+x+y))
306: (* t is in [11/16,19/16] *)
307: | 3 => (PIo4,zero,(y-x)/(x+y))
308: | 4 => (PIo4,zero,(y-x)/(x+y))
309: (* t is in [19/16,39/16] *)
310: | _ => let val z = y-x
311: val y = y+y+y
312: val t = x+x
313: in (at1fhi,at1flo,( (z+z)-x ) / ( t + y ))
314: end
315: (* end of if (t < 2.4375) *)
316: else let val t = if t <= big
317: (* t is in [2.4375, big] *)
318: then ~x / y
319: (* t is in [big, INF] *)
320: else zero
321: in (PIo2,zero,t) end
322: (* end of argument reduction *)
323: (* compute atan(t) for t in [~.4375, .4375] *)
324: val z = t*t
325: val z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+
326: z*(a9+z*(a10+z*a11)))))))))))
327: val z = lo - z
328: val z = z + t
329: val z = z + hi
330: in copysign(if signx>zero then z else PI-z,signy)
331: end
332: in
333: (* if x is 1.0, compute *)
334: if feql(x,one) then findatan2(x,copysign(y,one),y) handle Done r => r
335: (* when y = 0 *)
336: else if feql(y,zero) then if feql(signx,one) then y else copysign(PI,signy)
337: (* when x = 0 *)
338: else if feql(x,zero) then copysign(PIo2,signy)
339: else (* compute y/x *)
340: let val x = copysign(x,one)
341: val y = copysign(y,one)
342: val k = logb(y)
343: val m = InLine.-(k,logb(x))
344: val (t,y,x) = if InLine.>(m,60) then (big+big,y,x)
345: else if InLine.<(m,~80) then (y/x,y,x)
346: else (y/x,scalb(y,InLine.~ k),
347: scalb(x,InLine.~ k))
348: in findatan2(x,y,t) handle Done r => r end
349: end
350: end
351:
352: fun arctan(x:real) = atan2(x,one)
353:
354: fun sqrt(x: real) =
355: if x<=zero then if x<zero then raise Sqrt else x
356: else
357: let val k = 6 (* log base 2 of the precision *)
358: val n = InLine.rshift(logb(x),1)
359: val x = scalb(x, InLine.~(InLine.lshift(n,1)))
360: fun iter(0,g) = g
361: | iter(i,g) = iter(InLine.-(i,1), half * (g + x/g))
362: in scalb(iter(k,one),n)
363: end
364:
365: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.