Home My Page Projects Code Snippets Project Openings SML/NJ
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

[smlnj] Annotation of /sml/trunk/src/compiler/PervEnv/Basis/math64.sml
ViewVC logotype

Annotation of /sml/trunk/src/compiler/PervEnv/Basis/math64.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 89 - (view) (download)

1 : monnier 89 (* math64.sml
2 :     *
3 :     * COPYRIGHT (c) 1995 AT&T Bell Laboratories.
4 :     *
5 :     * The following functions were adapted from the 4.3BSD math library.
6 :     * Eventually, each machine supported should have a hand-coded math
7 :     * functor with more efficient versions of these functions.
8 :     *
9 :     ***************************************************************************
10 :     * *
11 :     * Copyright (c) 1985 Regents of the University of California. *
12 :     * *
13 :     * Use and reproduction of this software are granted in accordance with *
14 :     * the terms and conditions specified in the Berkeley Software License *
15 :     * Agreement (in particular, this entails acknowledgement of the programs' *
16 :     * source, and inclusion of this notice) with the additional understanding *
17 :     * that all recipients should regard themselves as participants in an *
18 :     * ongoing research project and hence should feel obligated to report *
19 :     * their experiences (good or bad) with these elementary function codes, *
20 :     * using "sendbug 4bsd-bugs@BERKELEY", to the authors. *
21 :     * *
22 :     * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. *
23 :     * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. *
24 :     * *
25 :     ***************************************************************************
26 :     *
27 :     *)
28 :    
29 :     structure Math64 : MATH =
30 :     struct
31 :    
32 :     type real = real
33 :    
34 :     infix 4 ==
35 :     val op + = InlineT.Real64.+
36 :     val op - = InlineT.Real64.-
37 :     val op * = InlineT.Real64.*
38 :     val op / = InlineT.Real64./
39 :     val op ~ = InlineT.Real64.~
40 :     val op < = InlineT.Real64.<
41 :     val op <= = InlineT.Real64.<=
42 :     val op > = InlineT.Real64.>
43 :     val op >= = InlineT.Real64.>=
44 :     val op == = InlineT.Real64.==
45 :    
46 :    
47 :     structure I = InlineT.DfltInt
48 :     val lessu : int * int -> bool = I.ltu
49 :    
50 :     val pi = 3.14159265358979323846
51 :     val e = 2.7182818284590452354
52 :    
53 :     fun isNan x = Bool.not(x==x)
54 :     val plusInfinity = 1E300 * 1E300
55 :     val minusInfinity = ~plusInfinity
56 :     val NaN = 0.0 / 0.0
57 :    
58 :     val two_to_the_54 = 18014398509481984.0
59 :     val two_to_the_minus_54 = 1.0 / 18014398509481984.0
60 :    
61 :     (* This function is IEEE double-precision specific;
62 :     it works correctly on subnormal inputs and outputs;
63 :     we do not apply it to inf's and nan's *)
64 :     fun scalb (x, k) =
65 :     let val j = Assembly.A.logb x
66 :     val k' = I.+(k,j)
67 :     in if j = ~1023
68 :     then scalb(x*two_to_the_54, I.-(k,54)) (*2*)
69 :     else if lessu(I.+(k',1022),2046)
70 :     then Assembly.A.scalb(x,k) (*1*)
71 :     else if I.<(k',0)
72 :     then if I.<(k',I.-(~1022,54))
73 :     then 0.0 (*3*)
74 :     else scalb(x,I.+(k,54)) *
75 :     two_to_the_minus_54 (*4*)
76 :     else x * plusInfinity (*5*)
77 :     end
78 :     (* Proof of correctness of scalb: (Appel)
79 :     1. if x is normal and x*2^k is normal
80 :     then case (*1*) applies, computes right answer
81 :     2. if x is subnormal and x*2^k is normal
82 :     then case (*2*) reduces problem to case 1.
83 :     3. if x*2^k is sub-subnormal (i.e. 0)
84 :     then case (*3*) applies, returns 0.0
85 :     4. if x*2^k is subnormal
86 :     then ~1076 < k' <= ~1023, case (*4*) applies,
87 :     computes right answer
88 :     5. if x*2^k is supernormal (i.e. infinity)
89 :     then case (*5*) computes right answer
90 :     *)
91 :    
92 :    
93 :    
94 :    
95 :    
96 :     (* This function is IEEE double-precision specific;
97 :     it works correctly on subnormal inputs;
98 :     must not be applied to inf's and nan's *)
99 :     fun logb(x) = (case Assembly.A.logb x
100 :     of ~1023 => (* denormalized number *)
101 :     I.-(Assembly.A.logb(x*two_to_the_54), 54)
102 :     | i => i
103 :     (* end case *))
104 :    
105 :     val negone = ~1.0
106 :     val zero = 0.0
107 :     val half = 0.5
108 :     val one = 1.0
109 :     val two = 2.0
110 :     val four = 4.0
111 :    
112 :     (** SHOULD BE INLINE OP **)
113 :     (* may be applied to inf's and nan's
114 :     GETS MINUS-ZERO WRONG!
115 :     *)
116 :     fun copysign(a,b) = (case (a<zero, b<zero)
117 :     of (true,true) => a
118 :     | (false,false) => a
119 :     | _ => ~a
120 :     (* end case *))
121 :    
122 :     (* may be applied to inf's and nan's *)
123 :     fun abs x = if x < zero then ~x else x
124 :    
125 :     fun op mod (a, b) = I.-(a, I.*(I.div(a, b), b))
126 :    
127 :     (* we will never call floor with an inf or nan *)
128 :     fun floor x = if x < 1073741824.0 andalso x >= ~1073741824.0
129 :     then Assembly.A.floor x
130 :     else if isNan x then raise General.Domain
131 :     else raise General.Overflow
132 :     val real = InlineT.real
133 :    
134 :     (* This is the IEEE double-precision maxint; won't work accurately on VAX *)
135 :     val maxint = 4503599627370496.0
136 :    
137 :     (* realround(x) returns x rounded to some nearby integer, almost always
138 :     * the nearest integer.
139 :     * May be applied to inf's and nan's.
140 :     *)
141 :     fun realround x = if x>=0.0 then x+maxint-maxint else x-maxint+maxint
142 :    
143 :     (* sin/cos *)
144 :     local
145 :     val S0 = ~1.6666666666666463126E~1
146 :     val S1 = 8.3333333332992771264E~3
147 :     val S2 = ~1.9841269816180999116E~4
148 :     val S3 = 2.7557309793219876880E~6
149 :     val S4 = ~2.5050225177523807003E~8
150 :     val S5 = 1.5868926979889205164E~10
151 :     in
152 :     fun sin__S z = (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))))
153 :     end
154 :    
155 :     local
156 :     val C0 = 4.1666666666666504759E~2
157 :     val C1 = ~1.3888888888865301516E~3
158 :     val C2 = 2.4801587269650015769E~5
159 :     val C3 = ~2.7557304623183959811E~7
160 :     val C4 = 2.0873958177697780076E~9
161 :     val C5 = ~1.1250289076471311557E~11
162 :     in
163 :     fun cos__C z = (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))))
164 :     end
165 :    
166 :     val PIo4 = 7.853981633974483096E~1
167 :     val PIo2 = 1.5707963267948966192E0
168 :     val PI3o4 = 2.3561944901923449288E0
169 :     val PI = pi
170 :     val PI2 = 6.2831853071795864769E0
171 :     val oneOver2Pi = 0.1591549430918953357688837633725143620345
172 :    
173 :     local
174 :     val thresh = 2.6117239648121182150E~1
175 :     in fun S y = y + y * sin__S(y*y)
176 :     fun C y =
177 :     let val yy = y*y
178 :     val c = cos__C yy
179 :     val Y = yy/two
180 :     in if Y < thresh then one - (Y - c)
181 :     else half - (Y - half - c)
182 :     end
183 :     end
184 :     fun sin x = (* This function propagages Inf's and Nan's correctly. *)
185 :     let (* x may be finite, inf, or nan at this point. *)
186 :     val xover2pi = x * oneOver2Pi
187 :     val x = PI2*(xover2pi - realround(xover2pi))
188 :     (* now, probably, ~pi <= x <= pi, except on vaxes.
189 :     x may be a nan, but cannot now be an inf. *)
190 :     fun lessPIo2 x = if x>PIo4 then C(PIo2-x) else S x
191 :     fun lessPI x = if x>PIo2 then lessPIo2(PI-x) else lessPIo2 x
192 :     fun positive x = if x>PI then sin(x-PI2) (* exceedingly rare *)
193 :     else lessPI x
194 :     in if x>=0.0
195 :     then positive x
196 :     else ~(positive(~x))
197 :     end
198 :    
199 :     fun cos x = sin(PIo2-x)
200 :    
201 :     fun tan x = sin x / cos x
202 :    
203 :     local
204 :     val p1 = 1.3887401997267371720E~2
205 :     val p2 = 3.3044019718331897649E~5
206 :     val q1 = 1.1110813732786649355E~1
207 :     val q2 = 9.9176615021572857300E~4
208 :     in fun exp__E(x:real,c:real) =
209 :     let val z = x*x
210 :     val p = z*(p1+z*p2)
211 :     val q = z*(q1+z*q2)
212 :     val xp= x*p
213 :     val xh= x*half
214 :     val w = xh-(q-xp)
215 :     val c = c+x*((xh*w-(q-(p+p+xp)))/(one-w)+c)
216 :     in z*half+c
217 :     end
218 :     end
219 :    
220 :     (* for exp and ln *)
221 :     val ln2hi = 6.9314718036912381649E~1
222 :     val ln2lo = 1.9082149292705877000E~10
223 :     val sqrt2 = 1.4142135623730951455E0
224 :     val lnhuge = 7.1602103751842355450E2
225 :     val lntiny = ~7.5137154372698068983E2
226 :     val invln2 = 1.4426950408889633870E0
227 :    
228 :     fun exp(x:real) = (* propagates and generates inf's and nan's correctly *)
229 :     let fun exp_norm x =
230 :     let (* argument reduction : x --> x - k*ln2 *)
231 :     val k = floor(invln2*x+copysign(half,x)) (* k=NINT(x/ln2) *)
232 :     val K = real k
233 :     (* express x-k*ln2 as z+c *)
234 :     val hi = x-K*ln2hi
235 :     val lo = K*ln2lo
236 :     val z = hi - lo
237 :     val c = (hi-z)-lo
238 :     (* return 2^k*[expm1(x) + 1] *)
239 :     val z = z + exp__E(z,c)
240 :     in scalb(z+one,k)
241 :     end
242 :     in if x <= lnhuge
243 :     then if x >= lntiny
244 :     then exp_norm x
245 :     else zero
246 :     else if isNan x then x else plusInfinity
247 :     end
248 :    
249 :     local
250 :     val L1 = 6.6666666666667340202E~1
251 :     val L2 = 3.9999999999416702146E~1
252 :     val L3 = 2.8571428742008753154E~1
253 :     val L4 = 2.2222198607186277597E~1
254 :     val L5 = 1.8183562745289935658E~1
255 :     val L6 = 1.5314087275331442206E~1
256 :     val L7 = 1.4795612545334174692E~1
257 :     in fun log__L(z) = z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))
258 :     end
259 :    
260 :     fun ln(x:real) = (* handles inf's and nan's correctly *)
261 :     if x>0.0
262 :     then if x < plusInfinity
263 :     then let
264 :     val k = logb(x)
265 :     val x = scalb(x, I.~ k)
266 :     val (k,x) = if x >= sqrt2 then (I.+(k, 1), x*half) else (k,x)
267 :     val K = real k
268 :     val x = x - one
269 :     (* compute log(1+x) *)
270 :     val s = x/(two+x)
271 :     val t = x*x*half
272 :     val z = K*ln2lo+s*(t+log__L(s*s))
273 :     val x = x + (z - t)
274 :     in
275 :     K*ln2hi+x
276 :     end
277 :     else x
278 :     else if (x == 0.0)
279 :     then minusInfinity
280 :     else if isNan x then x else NaN
281 :    
282 :     val oneOverln10 = 1.0 / ln 10.0
283 :    
284 :     fun log10 x = ln x * oneOverln10
285 :    
286 :     fun isInt y = realround(y)-y == 0.0
287 :     fun isOddInt(y) = isInt((y-1.0)*0.5)
288 :    
289 :     fun intpow(x,0) = 1.0
290 :     | intpow(x,y) = let val h = I.rshift(y,1)
291 :     val z = intpow(x,h)
292 :     val zz = z*z
293 :     in if y=I.+(h,h) then zz else x*zz
294 :     end
295 :    
296 :     (* SML/NJ doesn't properly handle negative zeros.
297 :     Also, the copysign function works incorrectly on negative zero.
298 :     The code for "pow" below should work correctly when these other
299 :     bugs are fixed. A. Appel, 5/8/97 *)
300 :     fun pow(x,y) = if y>0.0
301 :     then if y<plusInfinity
302 :     then if x > minusInfinity
303 :     then if x > 0.0
304 :     then exp(y*ln(x))
305 :     else if x == 0.0
306 :     then if isOddInt(y)
307 :     then x
308 :     else 0.0
309 :     else if isInt(y)
310 :     then intpow(x,floor(y+0.5))
311 :     else NaN
312 :     else if isNan x
313 :     then x
314 :     else if isOddInt(y)
315 :     then x
316 :     else plusInfinity
317 :     else let val ax = abs(x)
318 :     in if ax>1.0 then plusInfinity
319 :     else if ax<1.0 then 0.0
320 :     else NaN
321 :     end
322 :     else if y < 0.0
323 :     then if y>minusInfinity
324 :     then if x > minusInfinity
325 :     then if x > 0.0
326 :     then exp(y*ln(x))
327 :     else if x==0.0
328 :     then if isOddInt(y)
329 :     then copysign(plusInfinity,x)
330 :     else plusInfinity
331 :     else if isInt(y)
332 :     then 1.0 / intpow(x, floor(~y+0.5))
333 :     else NaN
334 :     else if isNan x
335 :     then x
336 :     else if isOddInt(y)
337 :     then ~0.0
338 :     else 0.0
339 :     else let val ax = abs(x)
340 :     in if ax>1.0 then 0.0
341 :     else if ax<1.0 then plusInfinity
342 :     else NaN
343 :     end
344 :     else if isNan y
345 :     then y
346 :     else 1.0
347 :     local
348 :     val athfhi = 4.6364760900080611433E~1
349 :     val athflo = 1.0147340032515978826E~18
350 :     val at1hi = 0.78539816339744830676
351 :     val at1lo = 1.11258708870781088040E~18
352 :     val a1 = 3.3333333333333942106E~1
353 :     val a2 = ~1.9999999999979536924E~1
354 :     val a3 = 1.4285714278004377209E~1
355 :     val a4 = ~1.1111110579344973814E~1
356 :     val a5 = 9.0908906105474668324E~2
357 :     val a6 = ~7.6919217767468239799E~2
358 :     val a7 = 6.6614695906082474486E~2
359 :     val a8 = ~5.8358371008508623523E~2
360 :     val a9 = 4.9850617156082015213E~2
361 :     val a10 = ~3.6700606902093604877E~2
362 :     val a11 = 1.6438029044759730479E~2
363 :    
364 :     fun atn(t,hi,lo) = (* for ~0.4375 <= t <= 0.4375 *)
365 :     let val z = t*t
366 :     in hi+(t+(lo-t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+
367 :     z*(a8+z*(a9+z*(a10+z*a11)))))))))))))
368 :     end
369 :    
370 :     fun atan(t) = (* 0 <= t <= 1 *)
371 :     if t <= 0.4375 then atn(t,zero,zero)
372 :     else if t <= 0.6875 then atn((t-half)/(one+half*t), athfhi, athflo)
373 :     else atn((t-one)/(one+t), at1hi,at1lo)
374 :    
375 :     fun atanpy y = (* y>=0 *)
376 :     if y>one then PIo2 - atan(one/y) else atan(y)
377 :    
378 :     fun atan2pypx(x,y) =
379 :     if y>x then PIo2 - atan(x/y) else atan(y/x)
380 :    
381 :     fun atan2py(x,y) =
382 :     if x > 0.0 then atan2pypx(x,y)
383 :     else if x == 0.0 andalso y == 0.0 then 0.0
384 :     else PI - atan2pypx(~x,y)
385 :    
386 :     in fun atan y = (* miraculously handles inf's and nan's correctly *)
387 :     if y<=0.0 then ~(atanpy(~y)) else atanpy y
388 :     fun atan2(y,x) = (* miraculously handles inf's and nan's correctly *)
389 :     if y>=0.0 then atan2py(x,y) else ~(atan2py(x,~y))
390 :     end
391 :    
392 :     fun sqrt(x: real) = (* handles inf's and nan's correctly *)
393 :     if x>zero
394 :     then if x < plusInfinity
395 :     then let
396 :     val k = 6 (* log base 2 of the precision *)
397 :     val n = I.rshift(logb x, 1)
398 :     val x = scalb(x, I.~(I.+(n, n)))
399 :     fun iter(0, g) = g
400 :     | iter(i, g) = iter(I.-(i, 1), half * (g + x/g))
401 :     in
402 :     scalb(iter(k,one),n)
403 :     end
404 :     else x
405 :     else if x<zero then NaN else x
406 :    
407 :    
408 :     fun asin x = atan2(x, sqrt(1.0-x*x))
409 :     fun acos x = 2.0 * atan2(sqrt((1.0-x)/(1.0+x)),1.0)
410 :    
411 :     fun cosh u = let val a = exp u in if a==0.0
412 :     then plusInfinity
413 :     else 0.5 * (a + 1.0 / a)
414 :     end
415 :     fun sinh u = let val a = exp u
416 :     in if a==0.0
417 :     then copysign(plusInfinity,u)
418 :     else 0.5 * (a - 1.0 / a)
419 :     end
420 :     fun tanh u = let val a = exp u
421 :     val b = 1.0 / a
422 :     in if a==0.0 then copysign(1.0,u)
423 :     else (a-b) / (a+b)
424 :     end
425 :     end
426 :    
427 :    
428 :     (*
429 :     * $Log: math64.sml,v $
430 :     * Revision 1.1.1.1 1998/04/08 18:40:04 george
431 :     * Version 110.5
432 :     *
433 :     *)

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0