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/real64.sml
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 89 - (view) (download)

1 : monnier 89 (* real64.sml
2 :     *
3 :     * COPYRIGHT (c) 1995 AT&T Bell Laboratories.
4 :     *
5 :     *)
6 :    
7 :     structure Real64 : REAL =
8 :     struct
9 :     structure I = InlineT.DfltInt
10 :    
11 :     structure Math = Math64
12 :    
13 :     infix 4 == !=
14 :     type real = real
15 :     val ~ = InlineT.Real64.~
16 :     val op + = InlineT.Real64.+
17 :     val op - = InlineT.Real64.-
18 :     val op * = InlineT.Real64.*
19 :     val op / = InlineT.Real64./
20 :     fun *+(a,b,c) = a*b+c
21 :     fun *-(a,b,c) = a*b-c
22 :    
23 :     val op > = InlineT.Real64.>
24 :     val op < = InlineT.Real64.<
25 :     val op >= = InlineT.Real64.>=
26 :     val op <= = InlineT.Real64.<=
27 :    
28 :     val op == = InlineT.Real64.==
29 :     val op != = InlineT.Real64.!=
30 :    
31 :     fun unordered(x,y) = Bool.not(x>y orelse x <= y)
32 :     fun ?= (x, y) = (x == y) orelse unordered(x, y)
33 :    
34 :     fun isNormal x = (case Assembly.A.logb x
35 :     of ~1023 => false (* 0.0 or subnormal *)
36 :     | 1024 => false (* inf or nan *)
37 :     | _ => true
38 :     (* end case *))
39 :    
40 :     (* The next three values are computed laboriously, partly to
41 :     * avoid problems with inaccurate string->float conversions
42 :     * in the compiler itself.
43 :     *)
44 :     val maxFinite = let
45 :     fun f(x,i) = if i=1023 then x else f(x*2.0, I.+(i, 1))
46 :     val y = f(1.0,0)
47 :     fun g(z, y, 0) = z
48 :     | g(z, y, i) = g(z+y, y*0.5, I.-(i, 1))
49 :     in
50 :     g(0.0,y,53)
51 :     end
52 :    
53 :     val minNormalPos = let
54 :     fun f(x) = let
55 :     val y = x * 0.5
56 :     in
57 :     if isNormal y then f y else x
58 :     end
59 :     in
60 :     f 1.0
61 :     end
62 :    
63 :     local
64 :     (* The x86 uses extended precision (80 bits) internally, therefore
65 :     * it is necessary to write out the result of r * 0.5 to get
66 :     * 64 bit precision.
67 :     *)
68 :     val mem = InlineT.PolyArray.array(1, minNormalPos)
69 :     val update = InlineT.PolyArray.update
70 :     val subscript = InlineT.PolyArray.chkSub
71 :     fun f () = let
72 :     val r = subscript(mem, 0)
73 :     val y = r * 0.5
74 :     in
75 :     update(mem, 0, y);
76 :     if subscript(mem, 0) == 0.0 then r else f ()
77 :     end
78 :     in
79 :     val minPos = f()
80 :     end
81 :    
82 :     val posInf = maxFinite * maxFinite
83 :     val negInf = ~posInf
84 :    
85 :     fun isFinite x = negInf < x andalso x < posInf
86 :     fun isNan x = Bool.not(x==x)
87 :     fun floor x = if x < 1073741824.0 andalso x >= ~1073741824.0
88 :     then Assembly.A.floor x
89 :     else if isNan x then raise General.Domain
90 :     else raise General.Overflow
91 :    
92 :     fun ceil n = I.-(~1,floor(~(n+1.0)))
93 :     fun trunc n = if n < 0.0 then ceil n else floor n
94 :     fun round x = floor(x+0.5) (* bug: does not do round-to-nearest *)
95 :    
96 :     (* This is the IEEE double-precision maxint *)
97 :     val maxInt = 4503599627370496.0
98 :    
99 :     local
100 :     (* realround mode x returns x rounded to the nearest integer using the
101 :     * given rounding mode.
102 :     * May be applied to inf's and nan's.
103 :     *)
104 :     fun realround mode x = let
105 :     val saveMode = IEEEReal.getRoundingMode ()
106 :     in
107 :     IEEEReal.setRoundingMode mode;
108 :     if x>=0.0 then x+maxInt-maxInt else x-maxInt+maxInt
109 :     before IEEEReal.setRoundingMode saveMode
110 :     end
111 :     in
112 :     val realFloor = realround IEEEReal.TO_NEGINF
113 :     val realCeil = realround IEEEReal.TO_POSINF
114 :     val realTrunc = realround IEEEReal.TO_ZERO
115 :     end
116 :    
117 :     val abs : real -> real = InlineT.Real64.abs
118 :     val fromInt : int -> real = InlineT.real
119 :    
120 :     (* bug: operates correctly but slowly *)
121 :     fun fromLargeInt(x : Int32.int) =
122 :     let val i = Int32.quot(x,2)
123 :     val j = Int32.-(x,Int32.+(i,i))
124 :     val i' = Int32.toInt i
125 :     val j' = Int32.toInt j
126 :     in fromInt(i')*2.0+fromInt(j')
127 :     end
128 :    
129 :     (* bug: only one rounding mode implemented *)
130 :     fun toInt IEEEReal.TO_NEGINF = floor
131 :     | toInt _ = raise Fail "toInt supports only NEGINF rounding mode now"
132 :    
133 :     (* bug: doesn't support full range of large ints *)
134 :     fun toLargeInt mode x = Int32.fromInt(toInt mode x)
135 :    
136 :     fun toLarge x = x
137 :     fun fromLarge _ x = x
138 :    
139 :     fun sign x = if (x < 0.0) then ~1 else if (x > 0.0) then 1
140 :     else if isNan x then raise Domain else 0
141 :     fun signBit x = (* Bug: negative zero not handled properly *)
142 :     Assembly.A.scalb(x, I.~(Assembly.A.logb x)) < 0.0
143 :    
144 :     fun sameSign (x, y) = signBit x = signBit y
145 :    
146 :     fun copySign(x,y) = (* may not work if x is Nan *)
147 :     if sameSign(x,y) then x else ~x
148 :    
149 :     fun compare(x,y) = if x<y then General.LESS else if x>y then General.GREATER
150 :     else if x == y then General.EQUAL
151 :     else raise IEEEReal.Unordered
152 :    
153 :     fun compareReal(x,y) =
154 :     if x<y then IEEEReal.LESS else if x>y then IEEEReal.GREATER
155 :     else if x == y then IEEEReal.EQUAL
156 :     else IEEEReal.UNORDERED
157 :    
158 :    
159 :     (** This proably needs to be reorganized **)
160 :     fun class x = (* does not distinguish between quiet and signalling NaN *)
161 :     if signBit x
162 :     then if x>negInf then if x == 0.0 then IEEEReal.ZERO
163 :     else if Assembly.A.logb x = ~1023
164 :     then IEEEReal.SUBNORMAL
165 :     else IEEEReal.NORMAL
166 :     else if x==x then IEEEReal.INF
167 :     else IEEEReal.NAN IEEEReal.QUIET
168 :     else if x<posInf then if x == 0.0 then IEEEReal.ZERO
169 :     else if Assembly.A.logb x = ~1023
170 :     then IEEEReal.SUBNORMAL
171 :     else IEEEReal.NORMAL
172 :     else if x==x then IEEEReal.INF
173 :     else IEEEReal.NAN IEEEReal.QUIET
174 :    
175 :     val radix = 2
176 :     val precision = 52
177 :    
178 :     val two_to_the_54 = 18014398509481984.0
179 :    
180 :     val two_to_the_neg_1000 =
181 :     let fun f(i,x) = if i=0 then x else f(I.-(i,1), x*0.5)
182 :     in f(1000, 1.0)
183 :     end
184 :    
185 :     (* AARGH! Our version of logb gives a value that's one less than the
186 :     rest of the world's logb functions.
187 :     We should fix this systematically some time. *)
188 :    
189 :     fun toManExp x =
190 :     case I.+(Assembly.A.logb x, 1)
191 :     of ~1023 => if x==0.0 then {man=x,exp=0}
192 :     else let val {man=m,exp=e} = toManExp(x*1048576.0)
193 :     in {man=m,exp=I.-(e,20)}
194 :     end
195 :     | 1024 => {man=x,exp=0}
196 :     | i => {man=Assembly.A.scalb(x,I.~ i),exp=i}
197 :    
198 :     fun fromManExp {man=m,exp=e:int} =
199 :     if (m >= 0.5 andalso m <= 1.0 orelse m <= ~0.5 andalso m >= ~1.0)
200 :     then if I.>(e, 1020)
201 :     then if I.>(e, 1050) then if m>0.0 then posInf else negInf
202 :     else let fun f(i,x) = if i=0 then x else f(I.-(i,1),x+x)
203 :     in f(I.-(e,1020), Assembly.A.scalb(m,1020))
204 :     end
205 :     else if I.<(e, I.~ 1020)
206 :     then if I.<(e, I.~ 1200) then 0.0
207 :     else let fun f(i,x) = if i=0 then x else f(I.-(i,1), x*0.5)
208 :     in f(I.-(1020,e), Assembly.A.scalb(m,I.~ 1020))
209 :     end
210 :     else Assembly.A.scalb(m,e) (* This is the common case! *)
211 :     else let val {man=m',exp=e'} = toManExp m
212 :     in fromManExp{man=m', exp=I.+(e',e)}
213 :     end
214 :    
215 :     (* whole and split could be implemented more efficiently if we had
216 :     * control over the rounding mode; but for now we don't.
217 :     *)
218 :     fun whole x = if x>0.0
219 :     then if x > 0.5
220 :     then x-0.5+maxInt-maxInt
221 :     else whole(x+1.0)-1.0
222 :     else if x<0.0
223 :     then if x < ~0.5
224 :     then x+0.5-maxInt+maxInt
225 :     else whole(x-1.0)+1.0
226 :     else x
227 :    
228 :     fun split x = let val w = whole x
229 :     val f = x-w
230 :     in if abs(f)==1.0
231 :     then {whole=w+f,frac=0.0}
232 :     else {whole=w, frac=f}
233 :     end
234 :    
235 :     fun realMod x = let
236 :     val f = x - whole x
237 :     in
238 :     if abs f == 1.0 then 0.0 else f
239 :     end
240 :    
241 :     fun rem(x,y) = y * #frac(split(x/y))
242 :    
243 :     fun checkFloat x = if x>negInf andalso x<posInf then x
244 :     else if isNan x then raise General.Div
245 :     else raise General.Overflow
246 :    
247 :     (** NOTE logb and scalb are also defined in math64.sml; do we need both??? **)
248 :     fun logb x = (case Assembly.A.logb x
249 :     of ~1023 => (* denormalized number *)
250 :     I.-(Assembly.A.logb(x * two_to_the_54), 54)
251 :     | i => i
252 :     (* end case *))
253 :    
254 :     fun scalb (x, k) = if I.ltu(I.+(k,1022),2046)
255 :     then Assembly.A.scalb(x,k)
256 :     else let val k1 = I.div(k, 2)
257 :     in
258 :     scalb(scalb(x, k1), I.-(k, k1))
259 :     end
260 :    
261 :     fun nextAfter _ = raise Fail "Real.nextAfter unimplemented"
262 :    
263 :     fun min(x,y) = if x<y orelse isNan y then x else y
264 :     fun max(x,y) = if x>y orelse isNan y then x else y
265 :    
266 :     fun toDecimal _ = raise Fail "Real.toDecimal unimplemented"
267 :     fun fromDecimal _ = raise Fail "Real.fromDecimal unimplemented"
268 :    
269 :     val fmt = RealFormat.fmtReal
270 :     val toString = fmt (StringCvt.GEN NONE)
271 :     val scan = NumScan.scanReal
272 :     val fromString = StringCvt.scanString scan
273 :    
274 :     end (* Real64 *)
275 :    
276 :    
277 :     (*
278 :     * $Log: real64.sml,v $
279 :     * Revision 1.1.1.1 1998/04/08 18:40:04 george
280 :     * Version 110.5
281 :     *
282 :     *)

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