Commit | Line | Data |
---|---|---|
7f918cf1 CE |
1 | (* Copyright (C) 2011-2014,2017 Matthew Fluet. |
2 | * Copyright (C) 2003-2007 Henry Cejtin, Matthew Fluet, Suresh | |
3 | * Jagannathan, and Stephen Weeks. | |
4 | * | |
5 | * MLton is released under a BSD-style license. | |
6 | * See the file MLton-LICENSE for details. | |
7 | *) | |
8 | ||
9 | functor Real (structure W: WORD_EXTRA | |
10 | structure R: | |
11 | sig | |
12 | include PRE_REAL | |
13 | val castToWord: real -> W.word | |
14 | val castFromWord: W.word -> real | |
15 | end): REAL_EXTRA = | |
16 | struct | |
17 | structure MLton = Primitive.MLton | |
18 | structure Prim = R | |
19 | local | |
20 | open IEEEReal | |
21 | in | |
22 | datatype float_class = datatype float_class | |
23 | datatype rounding_mode = datatype rounding_mode | |
24 | end | |
25 | infix 4 == != ?= | |
26 | type real = R.real | |
27 | ||
28 | local | |
29 | open Prim | |
30 | in | |
31 | val realSize = Int32.toInt realSize | |
32 | val exponentBias = Int32.toInt exponentBias | |
33 | val precision = Int32.toInt precision | |
34 | val radix = Int32.toInt radix | |
35 | end | |
36 | ||
37 | val signBits = Word.one | |
38 | val exponentSignificandBits = Word.- (Word.fromInt realSize, signBits) | |
39 | val significandBits = Word.- (Word.fromInt precision, Word.one) | |
40 | val exponentBits = Word.- (exponentSignificandBits, significandBits) | |
41 | ||
42 | local | |
43 | val mkMask : Word.word -> W.word = | |
44 | fn b => W.notb (W.<< (W.notb W.zero, b)) | |
45 | in | |
46 | val signMask = | |
47 | W.<< (mkMask signBits, exponentSignificandBits) | |
48 | val exponentMask = | |
49 | W.<< (mkMask exponentBits, significandBits) | |
50 | val significandMask = | |
51 | mkMask significandBits | |
52 | end | |
53 | ||
54 | val class : real -> float_class = | |
55 | fn r => | |
56 | let | |
57 | val w = R.castToWord r | |
58 | in | |
59 | if W.andb (w, exponentMask) = exponentMask | |
60 | then if W.andb (w, significandMask) = W.zero | |
61 | then IEEEReal.INF | |
62 | else IEEEReal.NAN | |
63 | else if W.andb (w, exponentMask) = W.zero | |
64 | then if W.andb (w, significandMask) = W.zero | |
65 | then IEEEReal.ZERO | |
66 | else IEEEReal.SUBNORMAL | |
67 | else IEEEReal.NORMAL | |
68 | end | |
69 | ||
70 | val toBits : real -> {sign: bool, exponent: W.word, significand: W.word} = | |
71 | fn r => | |
72 | let | |
73 | val w = R.castToWord r | |
74 | val significand = | |
75 | W.andb (w, significandMask) | |
76 | val exponent = | |
77 | W.>> (W.andb (w, exponentMask), significandBits) | |
78 | val sign = | |
79 | W.andb (w, signMask) = signMask | |
80 | in | |
81 | {sign = sign, | |
82 | exponent = exponent, | |
83 | significand = significand} | |
84 | end | |
85 | ||
86 | val fromBits : {sign: bool, exponent: W.word, significand: W.word} -> real = | |
87 | fn {sign, exponent, significand} => | |
88 | let | |
89 | val w = | |
90 | W.orb (if sign then W.<< (W.one, exponentSignificandBits) else W.zero, | |
91 | W.orb (W.andb (W.<< (exponent, significandBits), exponentMask), | |
92 | W.andb (significand, significandMask))) | |
93 | val r = R.castFromWord w | |
94 | in | |
95 | r | |
96 | end | |
97 | ||
98 | local | |
99 | open Prim | |
100 | in | |
101 | val op *+ = op *+ | |
102 | val op *- = op *- | |
103 | val op * = op * | |
104 | val op + = op + | |
105 | val op - = op - | |
106 | val op / = op / | |
107 | val op < = op < | |
108 | val op <= = op <= | |
109 | val op > = op > | |
110 | val op >= = op >= | |
111 | val ~ = ~ | |
112 | val abs = abs | |
113 | end | |
114 | ||
115 | local | |
116 | fun 'a make {fromRealUnsafe: 'a -> real, | |
117 | toRealUnsafe: real -> 'a, | |
118 | other : {precision: Primitive.Int32.int}} = | |
119 | if R.precision = #precision other | |
120 | then (fn (_: rounding_mode) => fromRealUnsafe, | |
121 | toRealUnsafe) | |
122 | else (fn (m: rounding_mode) => fn r => | |
123 | IEEEReal.withRoundingMode (m, fn () => fromRealUnsafe r), | |
124 | toRealUnsafe) | |
125 | in | |
126 | val (fromReal32,toReal32) = | |
127 | make {fromRealUnsafe = R.fromReal32Unsafe, | |
128 | toRealUnsafe = R.toReal32Unsafe, | |
129 | other = {precision = Primitive.Real32.precision}} | |
130 | val (fromReal64,toReal64) = | |
131 | make {fromRealUnsafe = R.fromReal64Unsafe, | |
132 | toRealUnsafe = R.toReal64Unsafe, | |
133 | other = {precision = Primitive.Real64.precision}} | |
134 | end | |
135 | local | |
136 | structure S = | |
137 | LargeReal_ChooseRealN | |
138 | (type 'a t = real -> 'a | |
139 | val fReal32 = toReal32 | |
140 | val fReal64 = toReal64) | |
141 | in | |
142 | val toLarge = S.f | |
143 | end | |
144 | local | |
145 | structure S = | |
146 | LargeReal_ChooseRealN | |
147 | (type 'a t = rounding_mode -> 'a -> real | |
148 | val fReal32 = fromReal32 | |
149 | val fReal64 = fromReal64) | |
150 | in | |
151 | val fromLarge = S.f | |
152 | end | |
153 | ||
154 | val negInf = R.castFromWord (W.orb (signMask, exponentMask)) | |
155 | val negOne = R.castFromWord (W.orb (signMask, W.<< (W.fromInt exponentBias, significandBits))) | |
156 | val negZero = R.castFromWord signMask | |
157 | val zero = R.castFromWord W.zero | |
158 | val minPos = R.castFromWord W.one | |
159 | val minNormalPos = R.castFromWord (W.<< (W.one, significandBits)) | |
160 | val half = R.castFromWord (W.<< (W.- (W.fromInt exponentBias, W.one), significandBits)) | |
161 | val one = R.castFromWord (W.<< (W.fromInt exponentBias, significandBits)) | |
162 | val two = R.castFromWord (W.<< (W.+ (W.fromInt exponentBias, W.one), significandBits)) | |
163 | val maxFinite = R.castFromWord (W.- (exponentMask, W.one)) | |
164 | val posInf = R.castFromWord exponentMask | |
165 | ||
166 | val nan = posInf + negInf | |
167 | val posNan = R.castFromWord (W.andb (R.castToWord nan, W.notb signMask)) | |
168 | val negNan = R.castFromWord (W.orb (R.castToWord nan, signMask)) | |
169 | ||
170 | ||
171 | fun isFinite r = | |
172 | abs r <= maxFinite | |
173 | ||
174 | val op == = Prim.== | |
175 | ||
176 | val op != = not o op == | |
177 | ||
178 | fun isNan r = r != r | |
179 | ||
180 | fun isNormal r = class r = NORMAL | |
181 | ||
182 | val op ?= = | |
183 | if MLton.Codegen.isAMD64 orelse MLton.Codegen.isLLVM orelse MLton.Codegen.isX86 | |
184 | then R.?= | |
185 | else | |
186 | fn (x, y) => | |
187 | case (class x, class y) of | |
188 | (NAN, _) => true | |
189 | | (_, NAN) => true | |
190 | | (ZERO, ZERO) => true | |
191 | | _ => R.== (x, y) | |
192 | ||
193 | fun min (x, y) = | |
194 | if x <= y then x | |
195 | else if x > y then y | |
196 | else if isNan y then x | |
197 | else y | |
198 | ||
199 | fun max (x, y) = | |
200 | if x >= y then x | |
201 | else if x < y then y | |
202 | else if isNan y then x | |
203 | else y | |
204 | ||
205 | fun sign (x: real): int = | |
206 | if x > zero then 1 | |
207 | else if x < zero then ~1 | |
208 | else if x == zero then 0 | |
209 | else raise Domain | |
210 | ||
211 | val signBit = #sign o toBits | |
212 | ||
213 | fun sameSign (x, y) = signBit x = signBit y | |
214 | ||
215 | fun copySign (x, y) = | |
216 | if sameSign (x, y) | |
217 | then x | |
218 | else ~ x | |
219 | ||
220 | local | |
221 | structure I = IEEEReal | |
222 | in | |
223 | fun compareReal (x, y) = | |
224 | if x < y then I.LESS | |
225 | else if x > y then I.GREATER | |
226 | else if x == y then I.EQUAL | |
227 | else I.UNORDERED | |
228 | end | |
229 | ||
230 | local | |
231 | structure I = IEEEReal | |
232 | structure G = General | |
233 | in | |
234 | fun compare (x, y) = | |
235 | case compareReal (x, y) of | |
236 | I.EQUAL => G.EQUAL | |
237 | | I.GREATER => G.GREATER | |
238 | | I.LESS => G.LESS | |
239 | | I.UNORDERED => raise IEEEReal.Unordered | |
240 | end | |
241 | ||
242 | fun unordered (x, y) = isNan x orelse isNan y | |
243 | ||
244 | (* nextAfter for subnormal and normal values works by converting | |
245 | * the real to a word of equivalent size and doing an increment | |
246 | * or decrement on the word. Because of the way IEEE floating | |
247 | * point numbers are represented, word {de,in}crement | |
248 | * automatically does the right thing at the boundary between | |
249 | * normals and denormals. Also, convienently, | |
250 | * maxFinite+1 = posInf and minFinite-1 = negInf. | |
251 | *) | |
252 | val nextAfter: real * real -> real = | |
253 | fn (r, t) => | |
254 | case (class r, class t) of | |
255 | (NAN, _) => nan | |
256 | | (_, NAN) => nan | |
257 | | (INF, _) => r | |
258 | | (ZERO, ZERO) => t (* want "t", not "r", to get the sign right *) | |
259 | | (ZERO, _) => if t > zero then minPos else ~minPos | |
260 | | _ => | |
261 | if r == t then | |
262 | r | |
263 | else if (r > t) = (r > zero) then | |
264 | R.castFromWord (W.- (R.castToWord r, W.one)) | |
265 | else | |
266 | R.castFromWord (W.+ (R.castToWord r, W.one)) | |
267 | ||
268 | local | |
269 | val one = One.make (fn () => ref (0 : C_Int.t)) | |
270 | in | |
271 | fun toManExp x = | |
272 | case class x of | |
273 | INF => {exp = 0, man = x} | |
274 | | NAN => {exp = 0, man = nan} | |
275 | | ZERO => {exp = 0, man = x} | |
276 | | _ => One.use (one, fn r => | |
277 | let | |
278 | val man = R.frexp (x, r) | |
279 | in | |
280 | {exp = C_Int.toInt (!r), man = man} | |
281 | end) | |
282 | end | |
283 | ||
284 | fun fromManExp {exp, man} = | |
285 | (R.ldexp (man, C_Int.fromInt exp)) | |
286 | handle Overflow => | |
287 | man * (if Int.< (exp, 0) then zero else posInf) | |
288 | ||
289 | val fromManExp = | |
290 | if MLton.Codegen.isX86 | |
291 | then fromManExp | |
292 | else | |
293 | fn {exp, man} => | |
294 | case class man of | |
295 | INF => man | |
296 | | NAN => man | |
297 | | ZERO => man | |
298 | | _ => fromManExp {exp = exp, man = man} | |
299 | ||
300 | local | |
301 | val one = One.make (fn () => ref zero) | |
302 | in | |
303 | fun split x = | |
304 | case class x of | |
305 | INF => {frac = if x > zero then zero else ~zero, | |
306 | whole = x} | |
307 | | NAN => {frac = nan, whole = nan} | |
308 | | _ => | |
309 | let | |
310 | val (frac, whole) = | |
311 | One.use (one, fn int => | |
312 | (R.modf (x, int), !int)) | |
313 | (* Some platforms' C libraries don't get sign of | |
314 | * zero right. | |
315 | *) | |
316 | fun fix y = | |
317 | if class y = ZERO andalso not (sameSign (x, y)) | |
318 | then ~ y | |
319 | else y | |
320 | in | |
321 | {frac = fix frac, | |
322 | whole = fix whole} | |
323 | end | |
324 | end | |
325 | ||
326 | val realMod = #frac o split | |
327 | ||
328 | fun checkFloat x = | |
329 | if isFinite x then x | |
330 | else if isNan x then raise Div | |
331 | else raise Overflow | |
332 | ||
333 | val realCeil = R.realCeil | |
334 | val realFloor = R.realFloor | |
335 | val realTrunc = R.realTrunc | |
336 | ||
337 | (* Unfortunately, libc round ties to zero instead of even values. *) | |
338 | (* Fortunately, if any rounding mode is supported, it's TO_NEAREST. *) | |
339 | val realRound = fn r => IEEEReal.withRoundingMode (TO_NEAREST, fn () => R.round r) | |
340 | ||
341 | fun rem (x, y) = | |
342 | (case class x of | |
343 | INF => nan | |
344 | | NAN => nan | |
345 | | ZERO => zero | |
346 | | _ => (case class y of | |
347 | INF => x | |
348 | | NAN => nan | |
349 | | ZERO => nan | |
350 | | _ => x - realTrunc (x/y) * y)) | |
351 | ||
352 | (* fromDecimal, scan, fromString: decimal -> binary conversions *) | |
353 | fun strtor (str: NullString.t, | |
354 | rounding_mode: IEEEReal.rounding_mode) = | |
355 | let | |
356 | val rounding : C_Int.int = | |
357 | case rounding_mode of | |
358 | TO_NEAREST => 1 | |
359 | | TO_NEGINF => 3 | |
360 | | TO_POSINF => 2 | |
361 | | TO_ZERO => 0 | |
362 | in | |
363 | Prim.strtor (str, rounding) | |
364 | end | |
365 | exception Bad | |
366 | fun fromDecimalWithRoundingMode | |
367 | ({class, digits, exp, sign}: IEEEReal.decimal_approx, | |
368 | rounding_mode: IEEEReal.rounding_mode) = | |
369 | let | |
370 | fun doit () = | |
371 | let | |
372 | val exp = | |
373 | if Int.< (exp, 0) | |
374 | then concat ["-", Int.toString (Int.~ exp)] | |
375 | else Int.toString exp | |
376 | (* | |
377 | val str = concat [if sign then "-" else "", | |
378 | "0.", digits, | |
379 | "E", exp, "\000"] | |
380 | *) | |
381 | val n = Int.+ (if sign then 1 else 0, | |
382 | Int.+ (4 (* "0." + "E" + "\000" *), | |
383 | Int.+ (List.length digits, | |
384 | String.size exp))) | |
385 | val a = Array.alloc n | |
386 | fun upd (i, c) = (Array.update (a, i, c); Int.+ (i, 1)) | |
387 | val i = 0 | |
388 | val i = if sign then upd (i, #"-") else i | |
389 | val i = upd (i, #"0") | |
390 | val i = upd (i, #".") | |
391 | val i = | |
392 | List.foldl | |
393 | (fn (d, i) => | |
394 | if Int.< (d, 0) orelse Int.> (d, 9) | |
395 | then raise Bad | |
396 | else upd (i, Char.chr (Int.+ (d, Char.ord #"0")))) | |
397 | i digits | |
398 | val i = upd (i, #"E") | |
399 | val i = CharVector.foldl (fn (c, i) => upd (i, c)) i exp | |
400 | val _ = upd (i, #"\000") | |
401 | val str = Vector.unsafeFromArray a | |
402 | val x = strtor (NullString.fromString str, rounding_mode) | |
403 | in | |
404 | x | |
405 | end | |
406 | in | |
407 | SOME (case class of | |
408 | INF => if sign then negInf else posInf | |
409 | | NAN => if sign then negNan else posNan | |
410 | | NORMAL => doit () | |
411 | | SUBNORMAL => doit () | |
412 | | ZERO => if sign then negZero else zero) | |
413 | handle Bad => NONE | |
414 | end | |
415 | ||
416 | fun fromDecimal da = fromDecimalWithRoundingMode (da, TO_NEAREST) | |
417 | ||
418 | fun scan reader state = | |
419 | case IEEEReal.scan reader state of | |
420 | NONE => NONE | |
421 | | SOME (da, state) => | |
422 | SOME (valOf (fromDecimalWithRoundingMode | |
423 | (da, IEEEReal.getRoundingMode ())), | |
424 | state) | |
425 | ||
426 | val fromString = StringCvt.scanString scan | |
427 | ||
428 | (* toDecimal, fmt, toString: binary -> decimal conversions. *) | |
429 | datatype mode = Fix | Gen | Sci | |
430 | local | |
431 | val one = One.make (fn () => ref (0: C_Int.int)) | |
432 | in | |
433 | fun gdtoa (x: real, mode: mode, ndig: int, | |
434 | rounding_mode: IEEEReal.rounding_mode) = | |
435 | let | |
436 | val mode : C_Int.int = | |
437 | case mode of | |
438 | Fix => 3 | |
439 | | Gen => 0 | |
440 | | Sci => 2 | |
441 | val ndig : C_Int.int = C_Int.fromInt ndig | |
442 | val rounding : C_Int.int = | |
443 | case rounding_mode of | |
444 | TO_NEAREST => 1 | |
445 | | TO_NEGINF => 3 | |
446 | | TO_POSINF => 2 | |
447 | | TO_ZERO => 0 | |
448 | in | |
449 | One.use (one, fn decpt => | |
450 | (Prim.gdtoa (x, mode, ndig, rounding, decpt), | |
451 | C_Int.toInt (!decpt))) | |
452 | end | |
453 | end | |
454 | ||
455 | fun toDecimal (x: real): IEEEReal.decimal_approx = | |
456 | case class x of | |
457 | INF => {class = INF, | |
458 | digits = [], | |
459 | exp = 0, | |
460 | sign = x < zero} | |
461 | | NAN => {class = NAN, | |
462 | digits = [], | |
463 | exp = 0, | |
464 | sign = signBit x} | |
465 | | ZERO => {class = ZERO, | |
466 | digits = [], | |
467 | exp = 0, | |
468 | sign = signBit x} | |
469 | | c => | |
470 | let | |
471 | val (cs, exp) = gdtoa (x, Gen, 0, TO_NEAREST) | |
472 | fun loop (i, ac) = | |
473 | if Int.< (i, 0) | |
474 | then ac | |
475 | else loop (Int.- (i, 1), | |
476 | (Int.- (Char.ord (CUtil.C_String.sub (cs, i)), | |
477 | Char.ord #"0")) | |
478 | :: ac) | |
479 | val digits = loop (Int.- (CUtil.C_String.length cs, 1), []) | |
480 | in | |
481 | {class = c, | |
482 | digits = digits, | |
483 | exp = exp, | |
484 | sign = x < zero} | |
485 | end | |
486 | ||
487 | datatype realfmt = datatype StringCvt.realfmt | |
488 | ||
489 | local | |
490 | fun fix (sign: string, cs: CUtil.C_String.t, decpt: int, ndig: int): string = | |
491 | let | |
492 | val length = CUtil.C_String.length cs | |
493 | in | |
494 | if Int.< (decpt, 0) | |
495 | then | |
496 | concat [sign, | |
497 | "0.", | |
498 | String.new (Int.~ decpt, #"0"), | |
499 | CUtil.C_String.toString cs, | |
500 | String.new (Int.+ (Int.- (ndig, length), | |
501 | decpt), | |
502 | #"0")] | |
503 | else | |
504 | let | |
505 | val whole = | |
506 | if decpt = 0 | |
507 | then "0" | |
508 | else | |
509 | String.tabulate (decpt, fn i => | |
510 | if Int.< (i, length) | |
511 | then CUtil.C_String.sub (cs, i) | |
512 | else #"0") | |
513 | in | |
514 | if 0 = ndig | |
515 | then concat [sign, whole] | |
516 | else | |
517 | let | |
518 | val frac = | |
519 | String.tabulate | |
520 | (ndig, fn i => | |
521 | let | |
522 | val j = Int.+ (i, decpt) | |
523 | in | |
524 | if Int.< (j, length) | |
525 | then CUtil.C_String.sub (cs, j) | |
526 | else #"0" | |
527 | end) | |
528 | in | |
529 | concat [sign, whole, ".", frac] | |
530 | end | |
531 | end | |
532 | end | |
533 | fun sci (x: real, ndig: int): string = | |
534 | let | |
535 | val sign = if x < zero then "~" else "" | |
536 | val (cs, decpt) = | |
537 | gdtoa (x, Sci, Int.+ (1, ndig), IEEEReal.getRoundingMode ()) | |
538 | val length = CUtil.C_String.length cs | |
539 | val whole = String.tabulate (1, fn _ => CUtil.C_String.sub (cs, 0)) | |
540 | val frac = | |
541 | if 0 = ndig | |
542 | then "" | |
543 | else concat [".", | |
544 | String.tabulate | |
545 | (ndig, fn i => | |
546 | let | |
547 | val j = Int.+ (i, 1) | |
548 | in | |
549 | if Int.< (j, length) | |
550 | then CUtil.C_String.sub (cs, j) | |
551 | else #"0" | |
552 | end)] | |
553 | val exp = Int.- (decpt, 1) | |
554 | val exp = | |
555 | let | |
556 | val (exp, sign) = | |
557 | if Int.< (exp, 0) | |
558 | then (Int.~ exp, "~") | |
559 | else (exp, "") | |
560 | in | |
561 | concat [sign, Int.toString exp] | |
562 | end | |
563 | in | |
564 | concat [sign, whole, frac, "E", exp] | |
565 | end | |
566 | fun gen (x: real, n: int): string = | |
567 | let | |
568 | val (prefix, x) = | |
569 | if x < zero | |
570 | then ("~", ~ x) | |
571 | else ("", x) | |
572 | val ss = Substring.full (sci (x, Int.- (n, 1))) | |
573 | fun isE c = c = #"E" | |
574 | fun isZero c = c = #"0" | |
575 | val expS = | |
576 | Substring.string (Substring.taker (not o isE) ss) | |
577 | val exp = valOf (Int.fromString expS) | |
578 | val man = | |
579 | String.translate | |
580 | (fn #"." => "" | c => str c) | |
581 | (Substring.string (Substring.dropr isZero | |
582 | (Substring.takel (not o isE) ss))) | |
583 | val manSize = String.size man | |
584 | fun zeros i = CharVector.tabulate (i, fn _ => #"0") | |
585 | fun dotAt i = | |
586 | concat [String.substring (man, 0, i), | |
587 | ".", String.extract (man, i, NONE)] | |
588 | fun sci () = concat [prefix, | |
589 | if manSize = 1 then man else dotAt 1, | |
590 | "E", expS] | |
591 | val op - = Int.- | |
592 | val op + = Int.+ | |
593 | val ~ = Int.~ | |
594 | val op >= = Int.>= | |
595 | in | |
596 | if exp >= (if manSize = 1 then 3 else manSize + 3) | |
597 | then sci () | |
598 | else if exp >= manSize - 1 | |
599 | then concat [prefix, man, zeros (exp - (manSize - 1))] | |
600 | else if exp >= 0 | |
601 | then concat [prefix, dotAt (exp + 1)] | |
602 | else if exp >= (if manSize = 1 then ~2 else ~3) | |
603 | then concat [prefix, "0.", zeros (~exp - 1), man] | |
604 | else sci () | |
605 | end | |
606 | in | |
607 | fun fmt spec = | |
608 | let | |
609 | val doit = | |
610 | case spec of | |
611 | EXACT => IEEEReal.toString o toDecimal | |
612 | | FIX opt => | |
613 | let | |
614 | val n = | |
615 | case opt of | |
616 | NONE => 6 | |
617 | | SOME n => | |
618 | if Primitive.Controls.safe andalso Int.< (n, 0) | |
619 | then raise Size | |
620 | else n | |
621 | in | |
622 | fn x => | |
623 | let | |
624 | val sign = if x < zero then "~" else "" | |
625 | val (cs, decpt) = | |
626 | gdtoa (x, Fix, n, IEEEReal.getRoundingMode ()) | |
627 | in | |
628 | fix (sign, cs, decpt, n) | |
629 | end | |
630 | end | |
631 | | GEN opt => | |
632 | let | |
633 | val n = | |
634 | case opt of | |
635 | NONE => 12 | |
636 | | SOME n => | |
637 | if Primitive.Controls.safe andalso Int.< (n, 1) | |
638 | then raise Size | |
639 | else n | |
640 | in | |
641 | fn x => gen (x, n) | |
642 | end | |
643 | | SCI opt => | |
644 | let | |
645 | val n = | |
646 | case opt of | |
647 | NONE => 6 | |
648 | | SOME n => | |
649 | if Primitive.Controls.safe andalso Int.< (n, 0) | |
650 | then raise Size | |
651 | else n | |
652 | in | |
653 | fn x => sci (x, n) | |
654 | end | |
655 | in | |
656 | fn x => | |
657 | case class x of | |
658 | NAN => (* if signBit x then "~nan" else *) "nan" | |
659 | | INF => if x > zero then "inf" else "~inf" | |
660 | | _ => doit x | |
661 | end | |
662 | end | |
663 | ||
664 | val toString = fmt (StringCvt.GEN NONE) | |
665 | ||
666 | (* Not all devices support all rounding modes. | |
667 | * However, every device has ceil/floor/round/trunc. | |
668 | *) | |
669 | fun safeConvert (m, cvt, x) = | |
670 | case m of | |
671 | TO_POSINF => cvt (realCeil x) | |
672 | | TO_NEGINF => cvt (realFloor x) | |
673 | | TO_NEAREST => cvt (realRound x) | |
674 | | TO_ZERO => cvt (realTrunc x) | |
675 | ||
676 | local | |
677 | fun 'a make {fromIntUnsafe: 'a -> real, | |
678 | toIntUnsafe: real -> 'a, | |
679 | other : {maxInt': Word.word -> 'a, | |
680 | minInt': 'a, | |
681 | precision': int}} = | |
682 | (fromIntUnsafe, | |
683 | if Int.< (precision, #precision' other) then | |
684 | let | |
685 | val trim = Int.- (Int.- (#precision' other, precision), 1) | |
686 | val maxInt' = (#maxInt' other) (Word.fromInt trim) | |
687 | val minInt' = #minInt' other | |
688 | val maxInt = fromIntUnsafe maxInt' | |
689 | val minInt = fromIntUnsafe minInt' | |
690 | in | |
691 | fn (m: rounding_mode) => fn x => | |
692 | if minInt <= x then | |
693 | if x <= maxInt then | |
694 | safeConvert (m, toIntUnsafe, x) | |
695 | else | |
696 | raise Overflow | |
697 | else | |
698 | if x < minInt then | |
699 | raise Overflow | |
700 | else | |
701 | raise Domain (* NaN *) | |
702 | end | |
703 | else | |
704 | let | |
705 | val maxInt' = (#maxInt' other) 0w0 | |
706 | val minInt' = #minInt' other | |
707 | val maxInt = fromIntUnsafe maxInt' | |
708 | val minInt = fromIntUnsafe minInt' | |
709 | in | |
710 | fn (m: rounding_mode) => fn x => | |
711 | if minInt <= x then | |
712 | if x <= maxInt then | |
713 | safeConvert (m, toIntUnsafe, x) | |
714 | else | |
715 | if x < maxInt + one then | |
716 | (case m of | |
717 | TO_NEGINF => maxInt' | |
718 | | TO_POSINF => raise Overflow | |
719 | | TO_ZERO => maxInt' | |
720 | | TO_NEAREST => | |
721 | (* Depends on maxInt being odd. *) | |
722 | if x - maxInt >= half then | |
723 | raise Overflow | |
724 | else | |
725 | maxInt') | |
726 | else | |
727 | raise Overflow | |
728 | else | |
729 | if x < minInt then | |
730 | if minInt - one < x then | |
731 | (case m of | |
732 | TO_NEGINF => raise Overflow | |
733 | | TO_POSINF => minInt' | |
734 | | TO_ZERO => minInt' | |
735 | | TO_NEAREST => | |
736 | (* Depends on minInt being even. *) | |
737 | if x - minInt < ~half then | |
738 | raise Overflow | |
739 | else | |
740 | minInt') | |
741 | else | |
742 | raise Overflow | |
743 | else | |
744 | raise Domain (* NaN *) | |
745 | end) | |
746 | in | |
747 | val (fromInt8,toInt8) = | |
748 | make {fromIntUnsafe = R.fromInt8Unsafe, | |
749 | toIntUnsafe = R.toInt8Unsafe, | |
750 | other = {maxInt' = fn w => Int8.<< (Int8.>> (Int8.maxInt', w), w), | |
751 | minInt' = Int8.minInt', | |
752 | precision' = Int8.precision'}} | |
753 | val (fromInt16,toInt16) = | |
754 | make {fromIntUnsafe = R.fromInt16Unsafe, | |
755 | toIntUnsafe = R.toInt16Unsafe, | |
756 | other = {maxInt' = fn w => Int16.<< (Int16.>> (Int16.maxInt', w), w), | |
757 | minInt' = Int16.minInt', | |
758 | precision' = Int16.precision'}} | |
759 | val (fromInt32,toInt32) = | |
760 | make {fromIntUnsafe = R.fromInt32Unsafe, | |
761 | toIntUnsafe = R.toInt32Unsafe, | |
762 | other = {maxInt' = fn w => Int32.<< (Int32.>> (Int32.maxInt', w), w), | |
763 | minInt' = Int32.minInt', | |
764 | precision' = Int32.precision'}} | |
765 | val (fromInt64,toInt64) = | |
766 | make {fromIntUnsafe = R.fromInt64Unsafe, | |
767 | toIntUnsafe = R.toInt64Unsafe, | |
768 | other = {maxInt' = fn w => Int64.<< (Int64.>> (Int64.maxInt', w), w), | |
769 | minInt' = Int64.minInt', | |
770 | precision' = Int64.precision'}} | |
771 | end | |
772 | ||
773 | val fromIntInf: IntInf.int -> real = | |
774 | fn i => | |
775 | let | |
776 | val str = | |
777 | if IntInf.< (i, 0) | |
778 | then "-" ^ (IntInf.toString (IntInf.~ i)) | |
779 | else IntInf.toString i | |
780 | val x = strtor (NullString.nullTerm str, | |
781 | IEEEReal.getRoundingMode ()) | |
782 | in | |
783 | x | |
784 | end | |
785 | ||
786 | val toIntInf: rounding_mode -> real -> LargeInt.int = | |
787 | fn mode => fn x => | |
788 | case class x of | |
789 | INF => raise Overflow | |
790 | | NAN => raise Domain | |
791 | | ZERO => (0 : LargeInt.int) | |
792 | | _ => | |
793 | let | |
794 | (* This round may turn x into an INF, so we need to check the | |
795 | * class again. | |
796 | *) | |
797 | val x = | |
798 | case mode of | |
799 | TO_POSINF => realCeil x | |
800 | | TO_NEGINF => realFloor x | |
801 | | TO_NEAREST => realRound x | |
802 | | TO_ZERO => realTrunc x | |
803 | in | |
804 | case class x of | |
805 | INF => raise Overflow | |
806 | | _ => valOf (IntInf.fromString (fmt (StringCvt.FIX (SOME 0)) x)) | |
807 | end | |
808 | ||
809 | local | |
810 | structure S = | |
811 | Int_ChooseInt | |
812 | (type 'a t = 'a -> real | |
813 | val fInt8 = fromInt8 | |
814 | val fInt16 = fromInt16 | |
815 | val fInt32 = fromInt32 | |
816 | val fInt64 = fromInt64 | |
817 | val fIntInf = fromIntInf) | |
818 | in | |
819 | val fromInt = S.f | |
820 | end | |
821 | local | |
822 | structure S = | |
823 | LargeInt_ChooseInt | |
824 | (type 'a t = 'a -> real | |
825 | val fInt8 = fromInt8 | |
826 | val fInt16 = fromInt16 | |
827 | val fInt32 = fromInt32 | |
828 | val fInt64 = fromInt64 | |
829 | val fIntInf = fromIntInf) | |
830 | in | |
831 | val fromLargeInt = S.f | |
832 | end | |
833 | local | |
834 | structure S = | |
835 | Int_ChooseInt | |
836 | (type 'a t = rounding_mode -> real -> 'a | |
837 | val fInt8 = toInt8 | |
838 | val fInt16 = toInt16 | |
839 | val fInt32 = toInt32 | |
840 | val fInt64 = toInt64 | |
841 | val fIntInf = toIntInf) | |
842 | in | |
843 | val toInt = S.f | |
844 | end | |
845 | local | |
846 | structure S = | |
847 | LargeInt_ChooseInt | |
848 | (type 'a t = rounding_mode -> real -> 'a | |
849 | val fInt8 = toInt8 | |
850 | val fInt16 = toInt16 | |
851 | val fInt32 = toInt32 | |
852 | val fInt64 = toInt64 | |
853 | val fIntInf = toIntInf) | |
854 | in | |
855 | val toLargeInt = S.f | |
856 | end | |
857 | ||
858 | val floor = toInt TO_NEGINF | |
859 | val ceil = toInt TO_POSINF | |
860 | val trunc = toInt TO_ZERO | |
861 | val round = toInt TO_NEAREST | |
862 | ||
863 | local | |
864 | fun 'a make {fromWordUnsafe: 'a -> real, | |
865 | toWordUnsafe: real -> 'a, | |
866 | other : {maxWord': Word.word -> 'a, | |
867 | wordSize: int, | |
868 | zeroWord: 'a}} = | |
869 | (fromWordUnsafe, | |
870 | if Int.<= (precision, #wordSize other) | |
871 | then let | |
872 | val trim = Int.- (#wordSize other, precision) | |
873 | val maxWord' = (#maxWord' other) (Word.fromInt trim) | |
874 | val maxWord = fromWordUnsafe maxWord' | |
875 | val zeroWord = #zeroWord other | |
876 | in | |
877 | fn (m: rounding_mode) => fn x => | |
878 | case class x of | |
879 | INF => raise Overflow | |
880 | | NAN => raise Domain | |
881 | | _ => if zero <= x | |
882 | then if x <= maxWord | |
883 | then safeConvert (m, toWordUnsafe, x) | |
884 | else raise Overflow | |
885 | else if x > ~one | |
886 | then (case m of | |
887 | TO_NEGINF => raise Overflow | |
888 | | TO_POSINF => zeroWord | |
889 | | TO_ZERO => zeroWord | |
890 | | TO_NEAREST => | |
891 | if x < ~half | |
892 | then raise Overflow | |
893 | else zeroWord) | |
894 | else raise Overflow | |
895 | end | |
896 | else let | |
897 | val maxWord' = (#maxWord' other) 0w0 | |
898 | val maxWord = fromWordUnsafe maxWord' | |
899 | val zeroWord = #zeroWord other | |
900 | in | |
901 | fn (m: rounding_mode) => fn x => | |
902 | case class x of | |
903 | INF => raise Overflow | |
904 | | NAN => raise Domain | |
905 | | _ => if zero <= x | |
906 | then if x <= maxWord | |
907 | then safeConvert (m, toWordUnsafe, x) | |
908 | else if x < maxWord + one | |
909 | then (case m of | |
910 | TO_NEGINF => maxWord' | |
911 | | TO_POSINF => raise Overflow | |
912 | | TO_ZERO => maxWord' | |
913 | | TO_NEAREST => | |
914 | (* Depends on maxWord being odd. *) | |
915 | if x - maxWord >= half | |
916 | then raise Overflow | |
917 | else maxWord') | |
918 | else raise Overflow | |
919 | else if x > ~one | |
920 | then (case m of | |
921 | TO_NEGINF => raise Overflow | |
922 | | TO_POSINF => zeroWord | |
923 | | TO_ZERO => zeroWord | |
924 | | TO_NEAREST => | |
925 | if x < ~half | |
926 | then raise Overflow | |
927 | else zeroWord) | |
928 | else raise Overflow | |
929 | end) | |
930 | in | |
931 | val (fromWord8,toWord8) = | |
932 | make {fromWordUnsafe = R.fromWord8Unsafe, | |
933 | toWordUnsafe = R.toWord8Unsafe, | |
934 | other = {maxWord' = fn w => Word8.<< (Word8.>> (Word8.maxWord', w), w), | |
935 | wordSize = Word8.wordSize, | |
936 | zeroWord = Word8.zero}} | |
937 | val (fromWord16,toWord16) = | |
938 | make {fromWordUnsafe = R.fromWord16Unsafe, | |
939 | toWordUnsafe = R.toWord16Unsafe, | |
940 | other = {maxWord' = fn w => Word16.<< (Word16.>> (Word16.maxWord', w), w), | |
941 | wordSize = Word16.wordSize, | |
942 | zeroWord = Word16.zero}} | |
943 | val (fromWord32,toWord32) = | |
944 | make {fromWordUnsafe = R.fromWord32Unsafe, | |
945 | toWordUnsafe = R.toWord32Unsafe, | |
946 | other = {maxWord' = fn w => Word32.<< (Word32.>> (Word32.maxWord', w), w), | |
947 | wordSize = Word32.wordSize, | |
948 | zeroWord = Word32.zero}} | |
949 | val (fromWord64,toWord64) = | |
950 | make {fromWordUnsafe = R.fromWord64Unsafe, | |
951 | toWordUnsafe = R.toWord64Unsafe, | |
952 | other = {maxWord' = fn w => Word64.<< (Word64.>> (Word64.maxWord', w), w), | |
953 | wordSize = Word64.wordSize, | |
954 | zeroWord = Word64.zero}} | |
955 | end | |
956 | ||
957 | local | |
958 | structure S = | |
959 | Word_ChooseWordN | |
960 | (type 'a t = 'a -> real | |
961 | val fWord8 = fromWord8 | |
962 | val fWord16 = fromWord16 | |
963 | val fWord32 = fromWord32 | |
964 | val fWord64 = fromWord64) | |
965 | in | |
966 | val fromWord = S.f | |
967 | end | |
968 | local | |
969 | structure S = | |
970 | LargeWord_ChooseWordN | |
971 | (type 'a t = 'a -> real | |
972 | val fWord8 = fromWord8 | |
973 | val fWord16 = fromWord16 | |
974 | val fWord32 = fromWord32 | |
975 | val fWord64 = fromWord64) | |
976 | in | |
977 | val fromLargeWord = S.f | |
978 | end | |
979 | local | |
980 | structure S = | |
981 | Word_ChooseWordN | |
982 | (type 'a t = rounding_mode -> real -> 'a | |
983 | val fWord8 = toWord8 | |
984 | val fWord16 = toWord16 | |
985 | val fWord32 = toWord32 | |
986 | val fWord64 = toWord64) | |
987 | in | |
988 | val toWord = S.f | |
989 | end | |
990 | local | |
991 | structure S = | |
992 | LargeWord_ChooseWordN | |
993 | (type 'a t = rounding_mode -> real -> 'a | |
994 | val fWord8 = toWord8 | |
995 | val fWord16 = toWord16 | |
996 | val fWord32 = toWord32 | |
997 | val fWord64 = toWord64) | |
998 | in | |
999 | val toLargeWord = S.f | |
1000 | end | |
1001 | ||
1002 | structure Math = | |
1003 | struct | |
1004 | open Prim.Math | |
1005 | ||
1006 | (* Patch functions to handle out-of-range args. Many C math | |
1007 | * libraries do not do what the SML Basis Spec requires. | |
1008 | *) | |
1009 | ||
1010 | local | |
1011 | fun patch f x = | |
1012 | if x < ~one orelse x > one | |
1013 | then nan | |
1014 | else f x | |
1015 | in | |
1016 | val acos = patch acos | |
1017 | val asin = patch asin | |
1018 | end | |
1019 | ||
1020 | local | |
1021 | fun patch f x = if x < zero then nan else f x | |
1022 | in | |
1023 | val ln = patch ln | |
1024 | val log10 = patch log10 | |
1025 | end | |
1026 | ||
1027 | (* The x86 doesn't get exp right on infs. *) | |
1028 | val exp = | |
1029 | if MLton.Codegen.isX86 | |
1030 | andalso let open MLton.Platform.Arch in host = X86 end | |
1031 | then (fn x => | |
1032 | case class x of | |
1033 | INF => if x > zero then posInf else zero | |
1034 | | _ => exp x) | |
1035 | else exp | |
1036 | ||
1037 | (* The Cygwin math library doesn't get pow right on some exceptional | |
1038 | * cases. | |
1039 | * | |
1040 | * The Linux math library doesn't get pow (x, y) right when x < 0 | |
1041 | * and y is large (but finite). | |
1042 | * | |
1043 | * So, we define a pow function that gives the correct result on | |
1044 | * exceptional cases, and only calls the C pow with x > 0. | |
1045 | *) | |
1046 | fun isInt (x: real): bool = x == realFloor x | |
1047 | ||
1048 | (* isEven x assumes isInt x. *) | |
1049 | fun isEven (x: real): bool = isInt (x / two) | |
1050 | ||
1051 | fun isOddInt x = isInt x andalso not (isEven x) | |
1052 | ||
1053 | fun isNeg x = x < zero | |
1054 | ||
1055 | fun pow (x, y) = | |
1056 | case class y of | |
1057 | INF => | |
1058 | if class x = NAN | |
1059 | then nan | |
1060 | else if x < negOne orelse x > one | |
1061 | then if isNeg y then zero else posInf | |
1062 | else if negOne < x andalso x < one | |
1063 | then if isNeg y then posInf else zero | |
1064 | else (* x = 1 orelse x = ~1 *) | |
1065 | nan | |
1066 | | NAN => nan | |
1067 | | ZERO => one | |
1068 | | _ => | |
1069 | (case class x of | |
1070 | INF => | |
1071 | if isNeg x | |
1072 | then if isNeg y | |
1073 | then if isOddInt y | |
1074 | then negZero | |
1075 | else zero | |
1076 | else if isOddInt y | |
1077 | then negInf | |
1078 | else posInf | |
1079 | else (* x = posInf *) | |
1080 | if isNeg y then zero else posInf | |
1081 | | NAN => nan | |
1082 | | ZERO => | |
1083 | if isNeg y | |
1084 | then if isOddInt y | |
1085 | then copySign (posInf, x) | |
1086 | else posInf | |
1087 | else if isOddInt y | |
1088 | then x | |
1089 | else zero | |
1090 | | _ => | |
1091 | if isNeg x | |
1092 | then if isInt y | |
1093 | then if isEven y | |
1094 | then Prim.Math.pow (~ x, y) | |
1095 | else negOne * Prim.Math.pow (~ x, y) | |
1096 | else nan | |
1097 | else Prim.Math.pow (x, y)) | |
1098 | ||
1099 | fun cosh x = | |
1100 | case class x of | |
1101 | INF => x | |
1102 | | ZERO => one | |
1103 | | _ => R.Math.cosh x | |
1104 | ||
1105 | fun sinh x = | |
1106 | case class x of | |
1107 | INF => x | |
1108 | | ZERO => x | |
1109 | | _ => R.Math.sinh x | |
1110 | ||
1111 | fun tanh x = | |
1112 | case class x of | |
1113 | INF => if x > zero then one else negOne | |
1114 | | ZERO => x | |
1115 | | _ => R.Math.tanh x | |
1116 | end | |
1117 | end | |
1118 | ||
1119 | structure Real32 = Real (structure W = Word32 | |
1120 | structure R = | |
1121 | struct | |
1122 | open Primitive.Real32 | |
1123 | local open Primitive.PackReal32 in | |
1124 | val castToWord = castToWord | |
1125 | val castFromWord = castFromWord | |
1126 | end | |
1127 | end) | |
1128 | structure Real64 = Real (structure W = Word64 | |
1129 | structure R = | |
1130 | struct | |
1131 | open Primitive.Real64 | |
1132 | local open Primitive.PackReal64 in | |
1133 | val castToWord = castToWord | |
1134 | val castFromWord = castFromWord | |
1135 | end | |
1136 | end) |