diff --git a/basis/Real.sml b/basis/Real.sml index c8e25496..2fc356c6 100644 --- a/basis/Real.sml +++ b/basis/Real.sml @@ -1,959 +1,803 @@ (* Title: Standard Basis Library: Real and Real32 structures. Author: David Matthews Copyright David Matthews 2000, 2005, 2008, 2016-18, 2023 This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License version 2.1 as published by the Free Software Foundation. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA *) local (* Functions common to Real and Real32 *) local open StringCvt in - (* Zero padding. Handle some of the shorter case to avoid too much concatentation. *) - fun padZero 0 = "" - | padZero 1 = "0" - | padZero 2 = "00" - | padZero 3 = "000" - | padZero 4 = "0000" - | padZero 5 = "00000" - | padZero 6 = "000000" - | padZero 7 = "0000000" - | padZero 8 = "00000000" - | padZero n = if n < 0 then raise Size else "00000000" ^ padZero(n-8) - (* How many digits will the mantissa take? It's possible to unroll this loop since the maximum for float is 9 and for double is 17. We want to avoid long-format arbitrary precision arithmetic as much as possible so we stop at 10 digits which is the maximum short-format power of ten in 32-bit mode. *) fun ndigits (i: LargeInt.int) = if i >= 1000000000 then ndigits(i div 1000000000) + 9 else if i >= 100000000 then 9 else if i >= 10000000 then 8 else if i >= 1000000 then 7 else if i >= 100000 then 6 else if i >= 10000 then 5 else if i >= 1000 then 4 else if i >= 100 then 3 else if i >= 10 then 2 else 1 - (* Power of ten - unroll a few values. *) - fun powerTen 0 = 1 - | powerTen 1 = 10 - | powerTen 2 = 100 - | powerTen 3 = 1000 - | powerTen 4 = 10000 - | powerTen 5 = 100000 - | powerTen 6 = 1000000 - | powerTen n = if n < 0 then raise Size else 1000000 * powerTen(n-6) - local (* PolyRealDoubleToDecimal returns an arbitrary precision number for the mantissa because it could be greater than 31 bits. In 64-bit mode it will always be short and we want to use fixed int arithmetic if we can. *) fun fixedDigitList (0, r) = r | fixedDigitList (n, r) = fixedDigitList (FixedInt.quot(n, 10), FixedInt.toInt(FixedInt.rem(n, 10)) :: r) in fun digitList(n, r) = if LibrarySupport.largeIntIsSmall n then fixedDigitList (FixedInt.fromLarge n, r) else let val (qu, rm) = IntInf.quotRem(n, 10) in digitList (qu, (Int.fromLarge rm) :: r) end end - datatype realConv = RCSpecial of string | RCNormal of bool * int * LargeInt.int - (* Common functions to convert real/Real32.real to strings *) (* "ndigs" is the number of digits after the decimal point. For sciFmt that means that "ndigs+1" digits are produced. For fixFmt at least "ndigs+1" digits are produced. If "ndigs" is 0 the value is printed as an integer without any ".0". In both cases ndigs > 0. *) - (* These functions start with the exact representation and round if necessary by - adding 0.5 to the last digit. Since the exact representation is itself a - rounded value it's possible that this could result in double rounding. *) fun exactFmt convert r = - case convert r of - RCSpecial s => s - | RCNormal (sign, exponent, mantissa) => + let + val {sign, exponent, mantissa, class} = convert r + open IEEEReal + in + case class of + NAN => "nan" + | INF => if sign then "~inf" else "inf" + | _ (* NORMAL, ZERO, SUBNORMAL *) => let val s = if sign then "~" else "" + val exponent = exponent + ndigits mantissa val (e, exp) = if exponent = 0 then ("", "") else ("E", Int.toString exponent) in String.concat[s, "0.", LargeInt.toString mantissa, e, exp] end - - and fixFmt convert ndigs r = - case convert r of - RCSpecial s => s - | RCNormal (sign, expo, mant) => - let - val signString = if sign then "~" else "" - val digits = ndigits mant - - val (roundedMantissa, exp) = - if digits-expo <= ndigs - then (* No rounding necessary *) (mant, expo) - else - let - val tens = powerTen(digits-expo-ndigs) - val rounded = (mant + tens div 2) div tens - in - (* If we have rounded to zero the exponent is zero. - We may also have rounded up a value of 9.999 - to add an extra digit. *) - if rounded = 0 - then (0, 0) - else (rounded, ndigits rounded - ndigs) - end - - val mantissa = LargeInt.toString roundedMantissa - val mantLength = String.size mantissa - in - if ndigs = 0 - then (* No decimal point or anything after. *) - ( - if exp >= mantLength - then String.concat[signString, mantissa, padZero(exp-mantLength)] - else if exp <= 0 - then String.concat[signString, "0"] - else String.concat[signString, String.substring(mantissa, 0, exp)] - ) - else if exp >= mantLength - then String.concat[signString, mantissa, padZero(exp-mantLength), ".", padZero ndigs] - else if exp <= 0 - then String.concat[signString, "0.", padZero(~exp), mantissa, - padZero(ndigs-mantLength+exp)] - else String.concat[signString, String.substring(mantissa, 0, exp), ".", - String.substring(mantissa, exp, mantLength-exp), padZero(ndigs-mantLength+exp)] - end - - (* sciFmt - always produces ndigs+1 significant digits *) - and sciFmt convert ndigs r = - case convert r of - RCSpecial s => s - | RCNormal (sign, expo, mant) => - let - val signString = if sign then "~" else "" - val digits = ndigits mant - val (roundedMantissa, exp) = - if digits <= ndigs+1 - then (* No rounding necessary *) (mant, expo-1) - else - let - val tens = powerTen(digits-ndigs-1) - val rounded = mant + tens div 2 - (* It's possible that this could increase the number of - digits and hence the exponent. e.g. 9.9999 -> 10.0 *) - in - (rounded, expo + ndigits rounded - digits - 1) - end - val mantissa = LargeInt.toString roundedMantissa - val mantLength = String.size mantissa - in - if ndigs = 0 - then (* No decimal point or anything after. *) - String.concat[signString, String.substring(mantissa, 0, 1), "E", Int.toString exp] - else String.concat[signString, String.substring(mantissa, 0, 1), ".", - String.substring(mantissa, 1, Int.min(mantLength-1, ndigs)), - padZero(Int.max(0, ndigs-mantLength+1)), "E", Int.toString exp] - end - - (* General format - produces up to ndigs of output. No trailing zeros are - produced except for any zeros before the DP. We also produce one ".0" if - necessary so that the result looks like a real number rather than an int. *) - and genFmt convert ndigs r = - case convert r of - RCSpecial s => s - | RCNormal (sign, expo, mant) => - let - val signString = if sign then "~" else "" - val digits = ndigits mant - - val (mantissa, exp) = - if digits <= ndigs - then (* No rounding necessary *) (LargeInt.toString mant, expo-1) - else - let - val tens = powerTen(digits-ndigs) - val rounded = mant + tens div 2 - (* It's possible that this could increase the number of - digits and hence the exponent. e.g. 9.9999 -> 10.0 - We need to remove any trailing zeros produced. *) - val asString = LargeInt.toString rounded - fun stripZeros 0 = 1 - | stripZeros n = - if String.sub(asString, n-1) = #"0" - then stripZeros(n-1) else n - val sLength = stripZeros ndigs - in - (String.substring(asString, 0, sLength), expo + ndigits rounded - digits - 1) - end - val mantLength = String.size mantissa (* <= ndigs *) - in - if exp > ndigs orelse exp < ~5 (* Use E format. No zero padding. *) - then - ( - if mantLength = 1 - then String.concat[signString, String.substring(mantissa, 0, 1), "E", Int.toString exp] - else String.concat[signString, String.substring(mantissa, 0, 1), ".", - String.substring(mantissa, 1, mantLength-1), "E", Int.toString exp] - ) - else (* Fixed format *) - if exp >= mantLength - then String.concat[signString, mantissa, padZero(exp+1-mantLength), ".0"] - else if exp+1 <= 0 - then String.concat[signString, "0.", padZero(~exp-1), mantissa] - else String.concat[signString, String.substring(mantissa, 0, exp+1), ".", - if mantLength = exp+1 then "0" else String.substring(mantissa, exp+1, mantLength-exp-1)] - end + end (* Note: The definition says, reasonably, that negative values for the number of digits raises Size. The tests also check for a very large value for the number of digits and seem to expect Size to be raised in that case. Note that the exception is raised when fmt spec is evaluated and before it is applied to an actual real argument. *) fun fmtFunction {sciFmt, ...} (SCI NONE) = sciFmt 6 | fmtFunction {sciFmt, ...} (SCI (SOME d) ) = if d < 0 orelse d > 200 then raise General.Size else sciFmt d | fmtFunction {fixFmt, ...} (FIX NONE) = fixFmt 6 | fmtFunction {fixFmt, ...} (FIX (SOME d) ) = if d < 0 orelse d > 200 then raise General.Size else fixFmt d | fmtFunction {genFmt, ...}(GEN NONE) = genFmt 12 | fmtFunction {genFmt, ...} (GEN (SOME d) ) = if d < 1 orelse d > 200 then raise General.Size else genFmt d | fmtFunction {exactFmt, ...} EXACT = exactFmt end open RealNumbersAsBits val floatMaxFiniteExp: FixedInt.int = 254 and doubleMaxFiniteExp: FixedInt.int = 2046 in structure Real: REAL = struct open IEEEReal val fromLargeInt: LargeInt.int -> real = Real.rtsCallFastI_R "PolyFloatArbitraryPrecision" val fromInt: int -> real = (* We have to select the appropriate conversion. This will be reduced down to the appropriate function but has to be type-correct whether int is arbitrary precision or fixed precision. Hence the "o Large/FixedInt.fromInt". *) if Bootstrap.intIsArbitraryPrecision then fromLargeInt o LargeInt.fromInt else Real.fromFixedInt o FixedInt.fromInt (* These are needed because we don't yet have conversion from string to real. They are filtered out by the signature. *) val zero = fromInt 0 and one = fromInt 1 and four = fromInt 4 type real = real (* Pick up from globals. *) structure Math: MATH = struct type real = real (* Pick up from globals. *) val sqrt = Real.rtsCallFastR_R "PolyRealSqrt" and sin = Real.rtsCallFastR_R "PolyRealSin" and cos = Real.rtsCallFastR_R "PolyRealCos" and atan = Real.rtsCallFastR_R "PolyRealArctan" and exp = Real.rtsCallFastR_R "PolyRealExp" and ln = Real.rtsCallFastR_R "PolyRealLog" and tan = Real.rtsCallFastR_R "PolyRealTan" and asin = Real.rtsCallFastR_R "PolyRealArcSin" and acos = Real.rtsCallFastR_R "PolyRealArcCos" and log10 = Real.rtsCallFastR_R "PolyRealLog10" and sinh = Real.rtsCallFastR_R "PolyRealSinh" and cosh = Real.rtsCallFastR_R "PolyRealCosh" and tanh = Real.rtsCallFastR_R "PolyRealTanh" val atan2 = Real.rtsCallFastRR_R "PolyRealAtan2" val pow = Real.rtsCallFastRR_R "PolyRealPow" (* Derived values. *) val e = exp one val pi = four * atan one end; infix 4 == != ?=; val op == = Real.== val op != : real * real -> bool = not o op == val radix = 2 val precision = 53 val maxFinite = doubleFromBinary{sign=false, exp=doubleMaxFiniteExp, mantissa = 0xFFFFFFFFFFFFF} val minNormalPos = doubleFromBinary{sign=false, exp=1, mantissa = 0} val minPos = doubleFromBinary{sign=false, exp=0, mantissa = 1} val posInf : real = one/zero val negInf : real = ~one/zero (* Real is LargeReal. *) fun toLarge (x: real) : (*LargeReal.*)real =x fun fromLarge (_ : IEEEReal.rounding_mode) (x: (*LargeReal.*)real): real = x (* isNan can be defined in terms of unordered. *) fun isNan x = Real.unordered(x, x) fun isFinite x = doubleExponent x <= doubleMaxFiniteExp (* This could be implemented using signBit and doubleFromBinary *) val copySign : (real * real) -> real = Real.rtsCallFastRR_R "PolyRealCopySign" val signBit = doubleSignBit fun isNormal x = let val exp = doubleExponent x in exp > 0 andalso exp <= doubleMaxFiniteExp end fun class x = let val exp = doubleExponent x in if exp > doubleMaxFiniteExp then ( if doubleMantissa x <> 0 then NAN else INF ) else if exp = 0 then ( if doubleMantissa x = 0 then ZERO else SUBNORMAL ) else NORMAL end fun sign x = if isNan x then raise General.Domain else if x == zero then 0 else if x < zero then ~1 else 1 fun sameSign (x, y) = signBit x = signBit y (* Returns the minimum. In the case where one is a NaN it returns the other. In that case the comparison will be false. *) fun min (a: real, b: real): real = if a < b orelse isNan b then a else b (* Similarly for max. *) fun max (a: real, b: real): real = if a > b orelse isNan b then a else b fun checkFloat x = if isFinite x then x else if isNan x then raise General.Div else raise General.Overflow local val frExp: real -> int * real = RunCall.rtsCallFull1 "PolyRealFrexp" val fromManAndExp: real*int -> real = Real.rtsCallFastRI_R "PolyRealLdexp" open Real in fun toManExp r = if not (isFinite r) orelse r == zero (* Nan, infinities and +/-0 all return r in the mantissa. We include 0 to preserve its sign. *) then {man=r, exp=0} else let val (exp, man) = frExp r in {man=man, exp=exp} end fun fromManExp {man, exp} = if not (isFinite man) orelse man == zero (* Nan, infinities and +/-0 in the mantissa all return their argument. *) then man else if LibrarySupport.isShortInt exp then fromManAndExp(man, exp) else (* Long arbitrary precision *) copySign(if Int.>(exp, 0) then posInf else zero, man) end (* Convert to integer. *) local (* The RTS function converts to at most a 64-bit value (even on 32-bits). That will convert all the bits of the mantissa but if the exponent is large we may have to multiply by some power of two. *) val realToInt: real -> LargeInt.int = RunCall.rtsCallFull1 "PolyRealBoxedToLongInt" (* These are defined to raise Domain rather than Overflow on Nans. *) fun checkNan x = if isNan x then raise Domain else x in val realFloor = Real.rtsCallFastR_R "PolyRealFloor" and realCeil = Real.rtsCallFastR_R "PolyRealCeil" and realTrunc = Real.rtsCallFastR_R "PolyRealTrunc" and realRound = Real.rtsCallFastR_R "PolyRealRound" fun toArbitrary x = if isNan x then raise General.Domain else if not (isFinite x) then raise General.Overflow else let val { man, exp } = toManExp x in if exp <= precision then realToInt x else IntInf.<< (realToInt(fromManExp{man=man, exp=precision}), Word.fromInt(exp - precision)) end fun toLargeInt IEEEReal.TO_NEGINF = toArbitrary o realFloor | toLargeInt IEEEReal.TO_POSINF = toArbitrary o realCeil | toLargeInt IEEEReal.TO_ZERO = toArbitrary o realTrunc | toLargeInt IEEEReal.TO_NEAREST = toArbitrary o realRound (* Conversions to FixedInt are put in by the compiler. If int is fixed we can use them otherwise we use the long versions. N.B. FixedInt.toInt is a no-op but is needed so this is type-correct when int is arbitrary. *) val floor = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toArbitrary o realFloor else FixedInt.toInt o Real.floorFix o checkNan and ceil = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toArbitrary o realCeil else FixedInt.toInt o Real.ceilFix o checkNan and trunc = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toArbitrary o realTrunc else FixedInt.toInt o Real.truncFix o checkNan and round = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toArbitrary o realRound else FixedInt.toInt o Real.roundFix o checkNan fun toInt IEEEReal.TO_NEGINF = floor | toInt IEEEReal.TO_POSINF = ceil | toInt IEEEReal.TO_ZERO = trunc | toInt IEEEReal.TO_NEAREST = round end; local val realConv: string->real = RunCall.rtsCallFull1 "PolyRealBoxedFromString" val posNan = abs(zero / zero) val negNan = ~posNan in fun fromDecimal { class = INF, sign=true, ...} = SOME negInf | fromDecimal { class = INF, sign=false, ...} = SOME posInf | fromDecimal { class = ZERO, sign=true, ...} = SOME (~ zero) | fromDecimal { class = ZERO, sign=false, ...} = SOME zero (* Generate signed Nans ignoring the digits and mantissa. There was code here to set the mantissa but there's no reference to that in the current version of the Basis library. *) | fromDecimal { class = NAN, sign=true, ... } = SOME negNan | fromDecimal { class = NAN, sign=false, ... } = SOME posNan | fromDecimal { class = _ (* NORMAL or SUBNORMAL *), sign, digits, exp} = (let fun toChar x = if x < 0 orelse x > 9 then raise General.Domain else Char.chr (x + Char.ord #"0") (* Turn the number into a string. *) val str = "0." ^ String.implode(List.map toChar digits) ^"E" ^ Int.toString exp (* Convert it to a real using the RTS conversion function. Change any Conversion exceptions into Domain. *) val result = realConv str handle RunCall.Conversion _ => raise General.Domain in if sign then SOME (~result) else SOME result end handle General.Domain => NONE ) end - local - (* We need to treat "nan" specially because IEEEReal.toString - is defined to return ~nan for negative nans whereas Real.fmt is defined always - to return "nan". This looks like an inconsistency in the definition but we follow - it. *) - fun realToRealConvert r = - let - val {sign, exponent, mantissa, class} = RealToDecimalConversion.doubleToMinimal r - in - case class of - NAN => RCSpecial "nan" - | INF => if sign then RCSpecial "~inf" else RCSpecial "inf" - | _ (* NORMAL, ZERO, SUBNORMAL *) => - RCNormal(sign, exponent + ndigits mantissa, mantissa) - end + fun toDecimal r = + let + val {sign, exponent, mantissa, class} = RealToDecimalConversion.doubleToMinimal r in - fun toDecimal r = - let - val {sign, exponent, mantissa, class} = RealToDecimalConversion.doubleToMinimal r - in - { class = class, sign = sign, exp = exponent + ndigits mantissa, digits = digitList(mantissa, []) } - end + { class = class, sign = sign, exp = exponent + ndigits mantissa, digits = digitList(mantissa, []) } + end + + local + val dble2str = RunCall.rtsCallFull3 "PolyRealDoubleToString" : real * char * int -> string - val fmt = fmtFunction { sciFmt=sciFmt realToRealConvert, fixFmt=fixFmt realToRealConvert, - genFmt=genFmt realToRealConvert, exactFmt=exactFmt realToRealConvert } - val toString = fmt (StringCvt.GEN NONE) + fun fixFmt ndigs r = dble2str(r, #"F", ndigs) + and sciFmt ndigs r = dble2str(r, #"E", ndigs) + and genFmt ndigs r = dble2str(r, #"G", ndigs) + in + val fmt = fmtFunction { sciFmt=sciFmt, fixFmt=fixFmt, + genFmt=genFmt, exactFmt=exactFmt RealToDecimalConversion.doubleToMinimal } end + val toString = fmt (StringCvt.GEN NONE) + (* Define these in terms of IEEEReal.scan since that deals with all the special cases. *) fun scan getc src = case IEEEReal.scan getc src of NONE => NONE | SOME (ieer, rest) => ( case fromDecimal ieer of NONE => NONE | SOME r => SOME(r, rest) ) val fromString = Option.composePartial (fromDecimal, IEEEReal.fromString) (* Converter to real values. This replaces the basic conversion function for reals installed in the bootstrap process. For more information see convInt in Int. *) local fun convReal (s: string) : real = let (* Set the rounding mode to TO_NEAREST whatever the current rounding mode. Otherwise the result of compiling a piece of code with a literal constant could depend on what the rounding mode was set to. We should always support TO_NEAREST. *) val oldRounding = IEEEReal.getRoundingMode() val () = IEEEReal.setRoundingMode IEEEReal.TO_NEAREST val scanResult = StringCvt.scanString scan s val () = IEEEReal.setRoundingMode oldRounding in case scanResult of NONE => raise RunCall.Conversion "Invalid real constant" | SOME res => res end in (* Install this as a conversion function for real literals. *) val (): unit = RunCall.addOverload convReal "convReal" end open Real (* Get the other definitions. *) fun compare (r1, r2) = if r1 == r2 then General.EQUAL else if r1 < r2 then General.LESS else if r1 > r2 then General.GREATER else raise Unordered fun compareReal (r1, r2) = if r1 == r2 then EQUAL else if r1 < r2 then LESS else if r1 > r2 then GREATER else UNORDERED (* This seems to be similar to == except that where == always returns false if either argument is a NaN this returns true. The implementation of == treats the unordered case specially so it would be possible to implement this in the same way. *) fun op ?= (x, y) = unordered(x, y) orelse x == y (* Although these may be built in in some architectures it's probably not worth treating them specially at the moment. *) fun *+ (x: real, y: real, z: real): real = x*y+z and *- (x: real, y: real, z: real): real = x*y-z val rem = Real.rtsCallFastRR_R "PolyRealRem" (* Split a real into whole and fractional parts. The fractional part must have the same sign as the number even if it is zero. *) fun split r = let val whole = realTrunc r val frac = r - whole in { whole = whole, frac = if not (isFinite r) then if isNan r then r else (* Infinity *) if r < zero then ~zero else zero else if frac == zero then if signBit r then ~zero else zero else frac } end (* Get the fractional part of a real. *) fun realMod r = #frac(split r) (* nextAfter: This was previously implemented in ML but, at the very least, needed to work with rounding to something other than TO_NEAREST. *) val nextAfter = Real.rtsCallFastRR_R "PolyRealNextAfter" end (* Real *) structure Real32: REAL where type real = Real32.real = (* Real32 uses some definitions from the Real structure above. *) struct open IEEEReal (* On both the X86 and ARM there is only a single conversion from double to float using the current rounding mode. If we want a specific rounding mode we need to set the rounding. *) fun fromLarge mode value = let val current = getRoundingMode() val () = setRoundingMode mode val result = Real32.fromReal value val () = setRoundingMode current in result end val fromRealRound = fromLarge TO_NEAREST (* Defined to use the current rounding mode. *) val fromLargeInt = Real32.fromReal o Real.fromLargeInt val fromInt: int -> Real32.real = (* We have to select the appropriate conversion. This will be reduced down to the appropriate function but has to be type-correct whether int is arbitrary precision or fixed precision. Hence the "o Large/FixedInt.fromInt". *) if Bootstrap.intIsArbitraryPrecision then fromLargeInt o LargeInt.fromInt else Real32.fromFixedInt o FixedInt.fromInt val zero = fromInt 0 and one = fromInt 1 and four = fromInt 4 val radix = 2 val precision = 24 val maxFinite = floatFromBinary{sign=false, exp=floatMaxFiniteExp, mantissa = 0x7FFFFF} val minNormalPos = floatFromBinary{sign=false, exp=1, mantissa = 0} val minPos = floatFromBinary{sign=false, exp=0, mantissa = 1} local open Real32 in val posInf : real = one/zero val negInf : real = ~one/zero val op != : real * real -> bool = not o op == end infix 4 == != ?=; (* isNan can be defined in terms of unordered. *) fun isNan x = Real32.unordered(x, x) fun isFinite x = floatExponent x <= floatMaxFiniteExp local open Real32 in val copySign : (real * real) -> real = rtsCallFastFF_F "PolyRealFCopySign" end val signBit = floatSignBit fun isNormal x = let val exp = floatExponent x in exp > 0 andalso exp <= floatMaxFiniteExp end fun class x = let val exp = floatExponent x in if exp > floatMaxFiniteExp then ( if floatMantissa x <> 0 then NAN else INF ) else if exp = 0 then ( if floatMantissa x = 0 then ZERO else SUBNORMAL ) else NORMAL end local open Real32 in fun sign x = if isNan x then raise General.Domain else if x == zero then 0 else if x < zero then ~1 else 1 end fun sameSign (x, y) = signBit x = signBit y local open Real32 in (* Returns the minimum. In the case where one is a NaN it returns the other. In that case the comparison will be false. *) fun min (a: real, b: real): real = if a < b orelse isNan b then a else b (* Similarly for max. *) fun max (a: real, b: real): real = if a > b orelse isNan b then a else b fun checkFloat x = if isFinite x then x else if isNan x then raise General.Div else raise General.Overflow (* On certain platforms e.g. mips, toLarge does not preserve the sign on nans. We deal with the non-finite cases here. *) (* Use the Real versions for the moment. *) fun toManExp r = if not (isFinite r) orelse r == zero (* Nan, infinities and +/-0 all return r in the mantissa. We include 0 to preserve its sign. *) then {man=r, exp=0} else let val {man, exp} = Real.toManExp(toLarge r) in {man = fromRealRound man, exp = exp } end and fromManExp {man, exp} = if not (isFinite man) orelse man == zero (* Nan, infinities and +/-0 in the mantissa all return their argument. *) then man else fromRealRound(Real.fromManExp{man=toLarge man, exp=exp}) fun compare (r1, r2) = if r1 == r2 then General.EQUAL else if r1 < r2 then General.LESS else if r1 > r2 then General.GREATER else raise Unordered fun compareReal (r1, r2) = if r1 == r2 then EQUAL else if r1 < r2 then LESS else if r1 > r2 then GREATER else UNORDERED fun op ?= (x, y) = unordered(x, y) orelse x == y (* Although these may be built in in some architectures it's probably not worth treating them specially at the moment. *) fun *+ (x: real, y: real, z: real): real = x*y+z and *- (x: real, y: real, z: real): real = x*y-z val realFloor = rtsCallFastF_F "PolyRealFFloor" and realCeil = rtsCallFastF_F "PolyRealFCeil" and realTrunc = rtsCallFastF_F "PolyRealFTrunc" and realRound = rtsCallFastF_F "PolyRealFRound" val rem = rtsCallFastFF_F "PolyRealFRem" (* Split a real into whole and fractional parts. The fractional part must have the same sign as the number even if it is zero. *) fun split r = let val whole = realTrunc r val frac = r - whole in { whole = whole, frac = if not (isFinite r) then if isNan r then r else (* Infinity *) if r < zero then ~zero else zero else if frac == zero then if signBit r then ~zero else zero else frac } end (* Get the fractional part of a real. *) fun realMod r = #frac(split r) val nextAfter = rtsCallFastFF_F "PolyRealFNextAfter" fun toLargeInt mode r = Real.toLargeInt mode (toLarge r) end local open Real32 (* These are defined to raise Domain rather than Overflow on Nans. *) fun checkNan x = if isNan x then raise Domain else x (* If int is fixed we use the hardware conversions otherwise we convert it to real and use the real to arbitrary conversions. *) in val floor = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toLargeInt IEEEReal.TO_NEGINF else FixedInt.toInt o floorFix o checkNan and ceil = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toLargeInt IEEEReal.TO_POSINF else FixedInt.toInt o ceilFix o checkNan and trunc = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toLargeInt IEEEReal.TO_ZERO else FixedInt.toInt o truncFix o checkNan and round = if Bootstrap.intIsArbitraryPrecision then LargeInt.toInt o toLargeInt IEEEReal.TO_NEAREST else FixedInt.toInt o roundFix o checkNan fun toInt IEEEReal.TO_NEGINF = floor | toInt IEEEReal.TO_POSINF = ceil | toInt IEEEReal.TO_ZERO = trunc | toInt IEEEReal.TO_NEAREST = round end (* Scan input source for a valid number. The format is the same as for double precision. Convert it using the current rounding mode. *) fun scan getc src = case Real.scan getc src of NONE => NONE | SOME (r, a) => SOME(Real32.fromReal r, a) val fromString = StringCvt.scanString scan - local - fun floatToRealConvert r = - let - val {sign, exponent, mantissa, class} = RealToDecimalConversion.floatToMinimal r - in - case class of - NAN => RCSpecial "nan" - | INF => if sign then RCSpecial "~inf" else RCSpecial "inf" - | _ (* NORMAL, ZERO, SUBNORMAL *) => - RCNormal(sign, exponent + ndigits mantissa, mantissa) - end + fun toDecimal r = + let + val {sign, exponent, mantissa, class} = RealToDecimalConversion.floatToMinimal r in - fun toDecimal r = - let - val {sign, exponent, mantissa, class} = RealToDecimalConversion.floatToMinimal r - in - { class = class, sign = sign, exp = exponent + ndigits mantissa, digits = digitList(mantissa, []) } - end + { class = class, sign = sign, exp = exponent + ndigits mantissa, digits = digitList(mantissa, []) } + end - val fmt = fmtFunction { sciFmt=sciFmt floatToRealConvert, fixFmt=fixFmt floatToRealConvert, - genFmt=genFmt floatToRealConvert, exactFmt=exactFmt floatToRealConvert } + local + (* Exact format must be defined specially for Real32.real but we can use the double conversion + for the other functions. *) + val dble2str = RunCall.rtsCallFull3 "PolyRealDoubleToString" : real * char * int -> string + + fun fixFmt ndigs r = dble2str(Real32.toLarge r, #"F", ndigs) + and sciFmt ndigs r = dble2str(Real32.toLarge r, #"E", ndigs) + and genFmt ndigs r = dble2str(Real32.toLarge r, #"G", ndigs) + in + val fmt = fmtFunction { sciFmt=sciFmt, fixFmt=fixFmt, + genFmt=genFmt, exactFmt=exactFmt RealToDecimalConversion.floatToMinimal } end val toString = fmt (StringCvt.GEN NONE) open Real32 (* Inherit the type and the built-in functions. *) (* Convert from decimal. This is defined to use TO_NEAREST. We need to handle NaNs specially because fromRealRound loses the sign on a NaN. *) local val posNan = abs(zero / zero) val negNan = ~posNan in fun fromDecimal { class = INF, sign=true, ...} = SOME negInf | fromDecimal { class = INF, sign=false, ...} = SOME posInf | fromDecimal { class = NAN, sign=true, ... } = SOME negNan | fromDecimal { class = NAN, sign=false, ... } = SOME posNan | fromDecimal arg = Option.map fromRealRound (Real.fromDecimal arg) end structure Math = struct type real = real val sqrt = rtsCallFastF_F "PolyRealFSqrt" and sin = rtsCallFastF_F "PolyRealFSin" and cos = rtsCallFastF_F "PolyRealFCos" and atan = rtsCallFastF_F "PolyRealFArctan" and exp = rtsCallFastF_F "PolyRealFExp" and ln = rtsCallFastF_F "PolyRealFLog" and tan = rtsCallFastF_F "PolyRealFTan" and asin = rtsCallFastF_F "PolyRealFArcSin" and acos = rtsCallFastF_F "PolyRealFArcCos" and log10 = rtsCallFastF_F "PolyRealFLog10" and sinh = rtsCallFastF_F "PolyRealFSinh" and cosh = rtsCallFastF_F "PolyRealFCosh" and tanh = rtsCallFastF_F "PolyRealFTanh" val atan2 = rtsCallFastFF_F "PolyRealFAtan2" val pow = rtsCallFastFF_F "PolyRealFPow" (* Derived values. *) val e = exp one val pi = four * atan one end (* Converter for literal constants. Copied from Real. *) local fun convReal (s: string) : real = let (* Set the rounding mode to TO_NEAREST whatever the current rounding mode. Otherwise the result of compiling a piece of code with a literal constant could depend on what the rounding mode was set to. We should always support TO_NEAREST. *) val oldRounding = IEEEReal.getRoundingMode() val () = IEEEReal.setRoundingMode IEEEReal.TO_NEAREST val scanResult = StringCvt.scanString scan s val () = IEEEReal.setRoundingMode oldRounding in case scanResult of NONE => raise RunCall.Conversion "Invalid real constant" | SOME res => res end in (* Install this as a conversion function for real literals. *) val (): unit = RunCall.addOverload convReal "convReal" end end (* Real32 *) end; structure Math = Real.Math; structure LargeReal: REAL = Real; (* Values available unqualified at the top-level. *) val real : int -> real = Real.fromInt val trunc : real -> int = Real.trunc val floor : real -> int = Real.floor val ceil : real -> int = Real.ceil val round : real -> int =Real.round; (* Overloads for Real32.real. The overloads for real were added in InitialBasis. *) val () = RunCall.addOverload Real32.>= ">=" and () = RunCall.addOverload Real32.<= "<=" and () = RunCall.addOverload Real32.> ">" and () = RunCall.addOverload Real32.< "<" and () = RunCall.addOverload Real32.+ "+" and () = RunCall.addOverload Real32.- "-" and () = RunCall.addOverload Real32.* "*" and () = RunCall.addOverload Real32.~ "~" and () = RunCall.addOverload Real32.abs "abs" and () = RunCall.addOverload Real32./ "/"; (* Install print functions. *) local fun print_real32 _ _ (r: Real32.real) = PolyML.PrettyString(Real32.fmt (StringCvt.GEN(SOME 7)) r) in val () = PolyML.addPrettyPrinter print_real32 end; local fun print_real _ _ (r: real) = PolyML.PrettyString(Real.fmt (StringCvt.GEN(SOME 10)) r) in val () = PolyML.addPrettyPrinter print_real; end; diff --git a/libpolyml/reals.cpp b/libpolyml/reals.cpp index d590b2ca..74d94ce3 100644 --- a/libpolyml/reals.cpp +++ b/libpolyml/reals.cpp @@ -1,940 +1,1005 @@ /* Title: Real number package. Author: Dave Matthews, Cambridge University Computer Laboratory Copyright (c) 2000 Cambridge University Technical Services Limited Further work copyright David C.J. Matthews 2011, 2016-19, 2023 This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License version 2.1 as published by the Free Software Foundation. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */ #ifdef HAVE_CONFIG_H #include "config.h" #elif defined(_WIN32) #include "winconfig.h" #else #error "No configuration file" #endif #ifdef HAVE_IEEEFP_H /* Other operating systems include "finite" in math.h, but Solaris doesn't? */ #include #endif #ifdef HAVE_FPU_CONTROL_H #include #endif #ifdef HAVE_FENV_H #include #endif #ifdef HAVE_FLOAT_H #include #endif #ifdef HAVE_MATH_H #include #endif #ifdef HAVE_STDIO_H #include #endif #ifdef HAVE_STRING_H #include #endif #ifdef HAVE_ERRNO_H #include #endif #ifdef HAVE_STDLIB_H #include #endif #ifdef HAVE_STDINT_H #include #endif #ifdef HAVE_ASSERT_H #include #define ASSERT(x) assert(x) #else #define ASSERT(x) #endif #include // Currently just for isnan. #include "globals.h" #include "run_time.h" #include "reals.h" #include "arb.h" #include "sys.h" #include "polystring.h" #include "save_vec.h" #include "rts_module.h" #include "machine_dep.h" #include "processes.h" #include "rtsentry.h" /* The Standard Basis Library assumes IEEE representation for reals. Among other things it does not permit equality on reals. That simplifies things considerably since we don't have to worry about there being two different representations of zero as 0 and ~0. We also don't need to check that the result is finite since NaN is allowed as a result. This code could do with being checked by someone who really understands IEEE floating point arithmetic. */ extern "C" { POLYEXTERNALSYMBOL POLYUNSIGNED PolyRealBoxedFromString(POLYUNSIGNED threadId, POLYUNSIGNED str); POLYEXTERNALSYMBOL POLYUNSIGNED PolyRealBoxedToLongInt(POLYUNSIGNED threadId, POLYUNSIGNED arg); + POLYEXTERNALSYMBOL POLYUNSIGNED PolyRealDoubleToString(POLYUNSIGNED threadId, POLYUNSIGNED arg, POLYUNSIGNED kind, POLYUNSIGNED prec); POLYEXTERNALSYMBOL double PolyRealSqrt(double arg); POLYEXTERNALSYMBOL double PolyRealSin(double arg); POLYEXTERNALSYMBOL double PolyRealCos(double arg); POLYEXTERNALSYMBOL double PolyRealArctan(double arg); POLYEXTERNALSYMBOL double PolyRealExp(double arg); POLYEXTERNALSYMBOL double PolyRealLog(double arg); POLYEXTERNALSYMBOL double PolyRealTan(double arg); POLYEXTERNALSYMBOL double PolyRealArcSin(double arg); POLYEXTERNALSYMBOL double PolyRealArcCos(double arg); POLYEXTERNALSYMBOL double PolyRealLog10(double arg); POLYEXTERNALSYMBOL double PolyRealSinh(double arg); POLYEXTERNALSYMBOL double PolyRealCosh(double arg); POLYEXTERNALSYMBOL double PolyRealTanh(double arg); POLYEXTERNALSYMBOL double PolyRealFloor(double arg); POLYEXTERNALSYMBOL double PolyRealCeil(double arg); POLYEXTERNALSYMBOL double PolyRealTrunc(double arg); POLYEXTERNALSYMBOL double PolyRealRound(double arg); POLYEXTERNALSYMBOL double PolyRealRem(double arg1, double arg2); POLYEXTERNALSYMBOL double PolyFloatArbitraryPrecision(POLYUNSIGNED arg); POLYEXTERNALSYMBOL POLYSIGNED PolyGetRoundingMode(POLYUNSIGNED); POLYEXTERNALSYMBOL POLYSIGNED PolySetRoundingMode(POLYUNSIGNED); POLYEXTERNALSYMBOL double PolyRealAtan2(double arg1, double arg2); POLYEXTERNALSYMBOL double PolyRealPow(double arg1, double arg2); POLYEXTERNALSYMBOL double PolyRealCopySign(double arg1, double arg2); POLYEXTERNALSYMBOL double PolyRealNextAfter(double arg1, double arg2); POLYEXTERNALSYMBOL double PolyRealLdexp(double arg1, POLYUNSIGNED arg2); POLYEXTERNALSYMBOL POLYUNSIGNED PolyRealFrexp(POLYUNSIGNED threadId, POLYUNSIGNED arg); POLYEXTERNALSYMBOL float PolyRealFSqrt(float arg); POLYEXTERNALSYMBOL float PolyRealFSin(float arg); POLYEXTERNALSYMBOL float PolyRealFCos(float arg); POLYEXTERNALSYMBOL float PolyRealFArctan(float arg); POLYEXTERNALSYMBOL float PolyRealFExp(float arg); POLYEXTERNALSYMBOL float PolyRealFLog(float arg); POLYEXTERNALSYMBOL float PolyRealFTan(float arg); POLYEXTERNALSYMBOL float PolyRealFArcSin(float arg); POLYEXTERNALSYMBOL float PolyRealFArcCos(float arg); POLYEXTERNALSYMBOL float PolyRealFLog10(float arg); POLYEXTERNALSYMBOL float PolyRealFSinh(float arg); POLYEXTERNALSYMBOL float PolyRealFCosh(float arg); POLYEXTERNALSYMBOL float PolyRealFTanh(float arg); POLYEXTERNALSYMBOL float PolyRealFAtan2(float arg1, float arg2); POLYEXTERNALSYMBOL float PolyRealFPow(float arg1, float arg2); POLYEXTERNALSYMBOL float PolyRealFCopySign(float arg1, float arg2); POLYEXTERNALSYMBOL float PolyRealFFloor(float arg); POLYEXTERNALSYMBOL float PolyRealFCeil(float arg); POLYEXTERNALSYMBOL float PolyRealFTrunc(float arg); POLYEXTERNALSYMBOL float PolyRealFRound(float arg); POLYEXTERNALSYMBOL float PolyRealFRem(float arg1, float arg2); POLYEXTERNALSYMBOL float PolyRealFNextAfter(float arg1, float arg2); } static Handle Real_convc(TaskData *mdTaskData, Handle str); // Positive and negative infinities and (positive) NaN. double posInf, negInf, notANumber; float posInfF, negInfF, notANumberF; /* Real numbers are represented by the address of the value. */ #define DBLE sizeof(double)/sizeof(POLYUNSIGNED) union db { double dble; POLYUNSIGNED words[DBLE]; }; double real_arg(Handle x) { union db r_arg_x; for (unsigned i = 0; i < DBLE; i++) { r_arg_x.words[i] = x->WordP()->Get(i).AsUnsigned(); } return r_arg_x.dble; } Handle real_result(TaskData *mdTaskData, double x) { union db argx; argx.dble = x; PolyObject *v = alloc(mdTaskData, DBLE, F_BYTE_OBJ); /* Copy as words in case the alignment is wrong. */ for(unsigned i = 0; i < DBLE; i++) { v->Set(i, PolyWord::FromUnsigned(argx.words[i])); } return mdTaskData->saveVec.push(v); } // We're using float for Real32 so it needs to be 32-bits. // Assume that's true for the moment. #if (SIZEOF_FLOAT != 4) #error "Float is not 32-bits. Please report this" #endif union flt { float fl; int32_t i; }; #if (SIZEOF_FLOAT < SIZEOF_POLYWORD) // Typically for 64-bit mode. Use a tagged representation. // The code-generator on the X86/64 assumes the float is in the // high order word. #define FLT_SHIFT ((SIZEOF_POLYWORD-SIZEOF_FLOAT)*8) float float_arg(Handle x) { union flt argx; argx.i = x->Word().AsSigned() >> FLT_SHIFT; return argx.fl; } Handle float_result(TaskData *mdTaskData, float x) { union flt argx; argx.fl = x; return mdTaskData->saveVec.push(PolyWord::FromSigned(((POLYSIGNED)argx.i << FLT_SHIFT) + 1)); } #else // Typically for 32-bit mode. Use a boxed representation. float float_arg(Handle x) { union flt argx; argx.i = (int32_t)x->WordP()->Get(0).AsSigned(); return argx.fl; } Handle float_result(TaskData *mdTaskData, float x) { union flt argx; argx.fl = x; PolyObject *v = alloc(mdTaskData, 1, F_BYTE_OBJ); v->Set(0, PolyWord::FromSigned(argx.i)); return mdTaskData->saveVec.push(v); } #endif POLYEXTERNALSYMBOL double PolyFloatArbitraryPrecision(POLYUNSIGNED arg) { return get_arbitrary_precision_as_real(PolyWord::FromUnsigned(arg)); } // Convert a boxed real to a long precision int. POLYUNSIGNED PolyRealBoxedToLongInt(POLYUNSIGNED threadId, POLYUNSIGNED arg) { TaskData *taskData = TaskData::FindTaskForId(threadId); ASSERT(taskData != 0); taskData->PreRTSCall(); Handle reset = taskData->saveVec.mark(); Handle pushedArg = taskData->saveVec.push(arg); double dx = real_arg(pushedArg); int64_t i = (int64_t)dx; Handle result = Make_arbitrary_precision(taskData, i); taskData->saveVec.reset(reset); taskData->PostRTSCall(); if (result == 0) return TAGGED(0).AsUnsigned(); else return result->Word().AsUnsigned(); } // RTS call for square-root. double PolyRealSqrt(double arg) { return sqrt(arg); } // RTS call for sine. double PolyRealSin(double arg) { return sin(arg); } // RTS call for cosine. double PolyRealCos(double arg) { return cos(arg); } // RTS call for arctan. double PolyRealArctan(double arg) { return atan(arg); } // RTS call for exp. double PolyRealExp(double arg) { return exp(arg); } // RTS call for ln. double PolyRealLog(double arg) { // Make sure the result conforms to the definition. // If the argument is a Nan each of the first two tests will fail. if (arg > 0.0) return log(arg); else if (arg == 0.0) // x may be +0.0 or -0.0 return negInf; // -infinity. else return notANumber; } // These were handled by the general dispatch function double PolyRealTan(double arg) { return tan(arg); } double PolyRealArcSin(double arg) { if (arg >= -1.0 && arg <= 1.0) return asin(arg); else return notANumber; } double PolyRealArcCos(double arg) { if (arg >= -1.0 && arg <= 1.0) return acos(arg); else return notANumber; } double PolyRealLog10(double arg) { // Make sure the result conforms to the definition. // If the argument is a Nan each of the first two tests will fail. if (arg > 0.0) return log10(arg); else if (arg == 0.0) // x may be +0.0 or -0.0 return negInf; // -infinity. else return notANumber; } double PolyRealSinh(double arg) { return sinh(arg); } double PolyRealCosh(double arg) { return cosh(arg); } double PolyRealTanh(double arg) { return tanh(arg); } double PolyRealFloor(double arg) { return floor(arg); } double PolyRealCeil(double arg) { return ceil(arg); } double PolyRealTrunc(double arg) { // Truncate towards zero if (arg >= 0.0) return floor(arg); else return ceil(arg); } double PolyRealRound(double arg) { // Round to nearest integral value. double drem = fmod(arg, 2.0); if (drem == 0.5 || drem == -1.5) // If the value was exactly positive even + 0.5 or // negative odd -0.5 round it down, otherwise round it up. return ceil(arg-0.5); else return floor(arg+0.5); } double PolyRealRem(double arg1, double arg2) { return fmod(arg1, arg2); } double PolyRealAtan2(double arg1, double arg2) { return atan2(arg1, arg2); } double PolyRealPow(double x, double y) { /* Some of the special cases are defined and don't seem to match the C pow function (at least as implemented in MS C). */ /* Maybe handle all this in ML? */ if (std::isnan(x)) { if (y == 0.0) return 1.0; else return notANumber; } else if (std::isnan(y)) return y; /* i.e. nan. */ else if (x == 0.0 && y < 0.0) { /* This case is not handled correctly in Solaris. It always returns -infinity. */ int iy = (int)floor(y); /* If x is -0.0 and y is an odd integer the result is -infinity. */ if (copysign(1.0, x) < 0.0 && (double)iy == y && (iy & 1)) return negInf; /* -infinity. */ else return posInf; /* +infinity. */ } return pow(x, y); } double PolyRealCopySign(double arg1, double arg2) { return copysign(arg1, arg2); } double PolyRealNextAfter(double arg1, double arg2) { return nextafter(arg1, arg2); } double PolyRealLdexp(double arg1, POLYUNSIGNED arg2) { POLYSIGNED exponent = PolyWord::FromUnsigned(arg2).UnTagged(); #if (SIZEOF_POLYWORD > SIZEOF_INT) // We've already checked for arbitrary precision values where necessary and // for zero and non-finite mantissa. Check the exponent fits in an int. if (exponent > 2 * DBL_MAX_EXP) return copysign(INFINITY, arg1); if (exponent < -2 * DBL_MAX_EXP) return copysign(0.0, arg1); #endif return ldexp(arg1, (int)exponent); } // Return the normalised fraction and the exponent. POLYUNSIGNED PolyRealFrexp(POLYUNSIGNED threadId, POLYUNSIGNED arg) { TaskData *taskData = TaskData::FindTaskForId(threadId); ASSERT(taskData != 0); taskData->PreRTSCall(); Handle reset = taskData->saveVec.mark(); Handle pushedArg = taskData->saveVec.push(arg); Handle result = 0; try { int exp = 0; // The value of exp is not always defined. Handle mantH = real_result(taskData, frexp(real_arg(pushedArg), &exp)); // Allocate a pair for the result result = alloc_and_save(taskData, 2); result->WordP()->Set(0, TAGGED(exp)); result->WordP()->Set(1, mantH->Word()); } catch (...) {} // If an ML exception is raised taskData->saveVec.reset(reset); taskData->PostRTSCall(); if (result == 0) return TAGGED(0).AsUnsigned(); else return result->Word().AsUnsigned(); } // RTS call for square-root. float PolyRealFSqrt(float arg) { return sqrtf(arg); } // RTS call for sine. float PolyRealFSin(float arg) { return sinf(arg); } // RTS call for cosine. float PolyRealFCos(float arg) { return cosf(arg); } // RTS call for arctan. float PolyRealFArctan(float arg) { return atanf(arg); } // RTS call for exp. float PolyRealFExp(float arg) { return expf(arg); } // RTS call for ln. float PolyRealFLog(float arg) { // Make sure the result conforms to the definition. // If the argument is a Nan each of the first two tests will fail. if (arg > 0.0) return logf(arg); else if (arg == 0.0) // x may be +0.0 or -0.0 return negInfF; // -infinity. else return notANumberF; } float PolyRealFTan(float arg) { return tanf(arg); } float PolyRealFArcSin(float arg) { if (arg >= -1.0 && arg <= 1.0) return asinf(arg); else return notANumberF; } float PolyRealFArcCos(float arg) { if (arg >= -1.0 && arg <= 1.0) return acosf(arg); else return notANumberF; } float PolyRealFLog10(float arg) { // Make sure the result conforms to the definition. // If the argument is a Nan each of the first two tests will fail. if (arg > 0.0) return log10f(arg); else if (arg == 0.0) // x may be +0.0 or -0.0 return negInfF; // -infinity. else return notANumberF; } float PolyRealFSinh(float arg) { return sinhf(arg); } float PolyRealFCosh(float arg) { return coshf(arg); } float PolyRealFTanh(float arg) { return tanhf(arg); } float PolyRealFAtan2(float arg1, float arg2) { return atan2f(arg1, arg2); } float PolyRealFPow(float x, float y) { /* Some of the special cases are defined and don't seem to match the C pow function (at least as implemented in MS C). */ /* Maybe handle all this in ML? */ if (std::isnan(x)) { if (y == 0.0) return 1.0; else return notANumberF; } else if (std::isnan(y)) return y; /* i.e. nan. */ else if (x == 0.0 && y < 0.0) { /* This case is not handled correctly in Solaris. It always returns -infinity. */ int iy = (int)floorf(y); /* If x is -0.0 and y is an odd integer the result is -infinity. */ if (copysign(1.0, x) < 0.0 && (float)iy == y && (iy & 1)) return negInfF; /* -infinity. */ else return posInfF; /* +infinity. */ } return powf(x, y); } float PolyRealFFloor(float arg) { return floorf(arg); } float PolyRealFCeil(float arg) { return ceilf(arg); } float PolyRealFTrunc(float arg) { // Truncate towards zero if (arg >= 0.0) return floorf(arg); else return ceilf(arg); } float PolyRealFRound(float arg) { // Round to nearest integral value. float drem = fmodf(arg, 2.0); if (drem == 0.5 || drem == -1.5) // If the value was exactly positive even + 0.5 or // negative odd -0.5 round it down, otherwise round it up. return ceilf(arg - 0.5f); else return floorf(arg + 0.5f); } float PolyRealFRem(float arg1, float arg2) { return fmodf(arg1, arg2); } float PolyRealFCopySign(float arg1, float arg2) { return copysignf(arg1, arg2); } float PolyRealFNextAfter(float arg1, float arg2) { return nextafterf(arg1, arg2); } /* CALL_IO1(Real_conv, REF, NOIND) */ Handle Real_convc(TaskData *mdTaskData, Handle str) /* string to real */ { double result; int i; char *finish; TempCString string_buffer(Poly_string_to_C_alloc(str->Word())); /* Scan the string turning '~' into '-' */ for(i = 0; string_buffer[i] != '\0'; i ++) { if (string_buffer[i] == '~') string_buffer[i] = '-'; } /* Now convert it */ result = strtod(string_buffer, &finish); // We no longer detect overflow and underflow and instead return // (signed) zeros for underflow and (signed) infinities for overflow. if (*finish != '\0') raise_exception_string(mdTaskData, EXC_conversion, ""); return real_result(mdTaskData, result); }/* Real_conv */ // Convert a string to a boxed real. This should really return an unboxed real. POLYUNSIGNED PolyRealBoxedFromString(POLYUNSIGNED threadId, POLYUNSIGNED str) { TaskData *taskData = TaskData::FindTaskForId(threadId); ASSERT(taskData != 0); taskData->PreRTSCall(); Handle reset = taskData->saveVec.mark(); Handle pushedString = taskData->saveVec.push(str); Handle result = 0; try { result = Real_convc(taskData, pushedString); } catch (...) { } // If an ML exception is raised taskData->saveVec.reset(reset); taskData->PostRTSCall(); if (result == 0) return TAGGED(0).AsUnsigned(); else return result->Word().AsUnsigned(); } +// Format a double-precision number using the format specifier. +POLYUNSIGNED PolyRealDoubleToString(POLYUNSIGNED threadId, POLYUNSIGNED arg, POLYUNSIGNED kind, POLYUNSIGNED prec) +{ + TaskData* taskData = TaskData::FindTaskForId(threadId); + ASSERT(taskData != 0); + taskData->PreRTSCall(); + Handle reset = taskData->saveVec.mark(); + Handle pushedArg = taskData->saveVec.push(arg); + Handle pushedPrec = taskData->saveVec.push(prec); + Handle result = 0; + char* format; + switch (UNTAGGED(PolyWord::FromUnsigned(kind))) + { + case 'e': case 'E': format = "%1.*E"; break; + case 'f': case 'F': format = "%1.*F"; break; + default: format = "%1.*G"; break; + } + + try { + int sSize = 100; + do { + TempCString buffer((char*)malloc(sSize)); + if ((char*)buffer == NULL) + raise_fail(taskData, "Insufficient memory for string"); + int precision = get_C_int(taskData, pushedPrec->Word()); + int retSize = snprintf((char* const)buffer, sSize, format, precision, real_arg(pushedArg)); + if (retSize < 0) + raise_fail(taskData, "Conversion error"); + if (retSize <= sSize) + { + char* p = buffer; + char* q = buffer; + // Replace - by ~ and remove any + signs. Suppress leading zeros after E. + bool suppressZero = false; + while (*p) + { + char c = *p++; + switch (c) + { + case '-': *q++ = '~'; break; + case '+': break; + case 'e': case 'E': *q++ = 'E'; suppressZero = true; break; + case '0': if (!suppressZero) *q++ = '0'; break; + default: *q++ = c; suppressZero = false; break; + } + } + if (suppressZero) *q++ = '0'; // We need at least one digit + *q = 0; + result = taskData->saveVec.push(C_string_to_Poly(taskData, buffer)); + break; + } + sSize *= 2; + } while (true); + + } + catch (...) {} // If an ML exception is raised + + taskData->saveVec.reset(reset); + taskData->PostRTSCall(); + if (result == 0) return TAGGED(0).AsUnsigned(); + else return result->Word().AsUnsigned(); +} + #if defined(__SOFTFP__) // soft-float lacks proper rounding mode support // While some systems will support fegetround/fesetround, it will have no // effect on the actual rounding performed, as the software implementation only // ever rounds to nearest. int getrounding() { return POLY_ROUND_TONEAREST; } int setrounding(int rounding) { switch (rounding) { case POLY_ROUND_TONEAREST: return 0; // The only mode supported } return -1; // Error - unsupported } // It would be nice to be able to use autoconf to test for these as functions // but they are frequently inlined #elif defined(HAVE_FENV_H) // C99 version. This is becoming the most common. int getrounding() { switch (fegetround()) { case FE_TONEAREST: return POLY_ROUND_TONEAREST; #ifndef HOSTARCHITECTURE_SH case FE_DOWNWARD: return POLY_ROUND_DOWNWARD; case FE_UPWARD: return POLY_ROUND_UPWARD; #endif case FE_TOWARDZERO: return POLY_ROUND_TOZERO; } return POLY_ROUND_TONEAREST; } int setrounding(int rounding) { switch (rounding) { case POLY_ROUND_TONEAREST: fesetround(FE_TONEAREST); return 0; // Choose nearest #ifndef HOSTARCHITECTURE_SH case POLY_ROUND_DOWNWARD: fesetround(FE_DOWNWARD); return 0; // Towards negative infinity case POLY_ROUND_UPWARD: fesetround(FE_UPWARD); return 0; // Towards positive infinity #endif case POLY_ROUND_TOZERO: fesetround(FE_TOWARDZERO); return 0; // Truncate towards zero default: return -1; } } #elif (defined(HAVE_IEEEFP_H) && ! defined(__CYGWIN__)) // Older FreeBSD. Cygwin has the ieeefp.h header but not the functions! int getrounding() { switch (fpgetround()) { case FP_RN: return POLY_ROUND_TONEAREST; case FP_RM: return POLY_ROUND_DOWNWARD; case FP_RP: return POLY_ROUND_UPWARD; case FP_RZ: return POLY_ROUND_TOZERO; default: return POLY_ROUND_TONEAREST; /* Shouldn't happen. */ } } int setrounding(int rounding) { switch (rounding) { case POLY_ROUND_TONEAREST: fpsetround(FP_RN); break; /* Choose nearest */ case POLY_ROUND_DOWNWARD: fpsetround(FP_RM); break; /* Towards negative infinity */ case POLY_ROUND_UPWARD: fpsetround(FP_RP); break; /* Towards positive infinity */ case POLY_ROUND_TOZERO: fpsetround(FP_RZ); break; /* Truncate towards zero */ } return 0 } #elif defined(_WIN32) // Windows version int getrounding() { switch (_controlfp(0,0) & _MCW_RC) { case _RC_NEAR: return POLY_ROUND_TONEAREST; case _RC_DOWN: return POLY_ROUND_DOWNWARD; case _RC_UP: return POLY_ROUND_UPWARD; case _RC_CHOP: return POLY_ROUND_TOZERO; } return POLY_ROUND_TONEAREST; } int setrounding(int rounding) { switch (rounding) { case POLY_ROUND_TONEAREST: _controlfp(_RC_NEAR, _MCW_RC); return 0; // Choose nearest case POLY_ROUND_DOWNWARD: _controlfp(_RC_DOWN, _MCW_RC); return 0; // Towards negative infinity case POLY_ROUND_UPWARD: _controlfp(_RC_UP, _MCW_RC); return 0; // Towards positive infinity case POLY_ROUND_TOZERO: _controlfp(_RC_CHOP, _MCW_RC); return 0; // Truncate towards zero } return -1; } #elif defined(_FPU_GETCW) && defined(_FPU_SETCW) // Older Linux version int getrounding() { fpu_control_t ctrl; _FPU_GETCW(ctrl); switch (ctrl & _FPU_RC_ZERO) { case _FPU_RC_NEAREST: return POLY_ROUND_TONEAREST; case _FPU_RC_DOWN: return POLY_ROUND_DOWNWARD; case _FPU_RC_UP: return POLY_ROUND_UPWARD; case _FPU_RC_ZERO: return POLY_ROUND_TOZERO; } return POLY_ROUND_TONEAREST; /* Never reached but this avoids warning message. */ } int setrounding(int rounding) { fpu_control_t ctrl; _FPU_GETCW(ctrl); ctrl &= ~_FPU_RC_ZERO; /* Mask off any existing rounding. */ switch (rounding) { case POLY_ROUND_TONEAREST: ctrl |= _FPU_RC_NEAREST; case POLY_ROUND_DOWNWARD: ctrl |= _FPU_RC_DOWN; case POLY_ROUND_UPWARD: ctrl |= _FPU_RC_UP; case POLY_ROUND_TOZERO: ctrl |= _FPU_RC_ZERO; } _FPU_SETCW(ctrl); return 0; } #else // Give up. Assume that we only support TO_NEAREST int getrounding() { return POLY_ROUND_TONEAREST; } int setrounding(int rounding) { if (rounding == POLY_ROUND_TONEAREST) return 0; else return -1; } #endif POLYSIGNED PolyGetRoundingMode(POLYUNSIGNED) { // Get the rounding and turn the result into a tagged integer. return TAGGED(getrounding()).AsSigned(); } POLYSIGNED PolySetRoundingMode(POLYUNSIGNED arg) { return TAGGED(setrounding((int)PolyWord::FromUnsigned(arg).UnTagged())).AsSigned(); } struct _entrypts realsEPT[] = { { "PolyRealBoxedFromString", (polyRTSFunction)&PolyRealBoxedFromString}, { "PolyRealBoxedToLongInt", (polyRTSFunction)&PolyRealBoxedToLongInt}, + { "PolyRealDoubleToString", (polyRTSFunction)&PolyRealDoubleToString}, { "PolyRealSqrt", (polyRTSFunction)&PolyRealSqrt}, { "PolyRealSin", (polyRTSFunction)&PolyRealSin}, { "PolyRealCos", (polyRTSFunction)&PolyRealCos}, { "PolyRealArctan", (polyRTSFunction)&PolyRealArctan}, { "PolyRealExp", (polyRTSFunction)&PolyRealExp}, { "PolyRealLog", (polyRTSFunction)&PolyRealLog}, { "PolyRealTan", (polyRTSFunction)&PolyRealTan}, { "PolyRealArcSin", (polyRTSFunction)&PolyRealArcSin}, { "PolyRealArcCos", (polyRTSFunction)&PolyRealArcCos}, { "PolyRealLog10", (polyRTSFunction)&PolyRealLog10}, { "PolyRealSinh", (polyRTSFunction)&PolyRealSinh}, { "PolyRealCosh", (polyRTSFunction)&PolyRealCosh}, { "PolyRealTanh", (polyRTSFunction)&PolyRealTanh}, { "PolyRealFloor", (polyRTSFunction)&PolyRealFloor}, { "PolyRealCeil", (polyRTSFunction)&PolyRealCeil}, { "PolyRealTrunc", (polyRTSFunction)&PolyRealTrunc}, { "PolyRealRound", (polyRTSFunction)&PolyRealRound}, { "PolyRealRem", (polyRTSFunction)&PolyRealRem }, { "PolyFloatArbitraryPrecision", (polyRTSFunction)&PolyFloatArbitraryPrecision}, { "PolyGetRoundingMode", (polyRTSFunction)&PolyGetRoundingMode}, { "PolySetRoundingMode", (polyRTSFunction)&PolySetRoundingMode}, { "PolyRealAtan2", (polyRTSFunction)&PolyRealAtan2 }, { "PolyRealPow", (polyRTSFunction)&PolyRealPow }, { "PolyRealCopySign", (polyRTSFunction)&PolyRealCopySign }, { "PolyRealNextAfter", (polyRTSFunction)&PolyRealNextAfter }, { "PolyRealLdexp", (polyRTSFunction)&PolyRealLdexp }, { "PolyRealFrexp", (polyRTSFunction)&PolyRealFrexp }, { "PolyRealFSqrt", (polyRTSFunction)&PolyRealFSqrt }, { "PolyRealFSin", (polyRTSFunction)&PolyRealFSin }, { "PolyRealFCos", (polyRTSFunction)&PolyRealFCos }, { "PolyRealFArctan", (polyRTSFunction)&PolyRealFArctan }, { "PolyRealFExp", (polyRTSFunction)&PolyRealFExp }, { "PolyRealFLog", (polyRTSFunction)&PolyRealFLog }, { "PolyRealFTan", (polyRTSFunction)&PolyRealFTan }, { "PolyRealFArcSin", (polyRTSFunction)&PolyRealFArcSin }, { "PolyRealFArcCos", (polyRTSFunction)&PolyRealFArcCos }, { "PolyRealFLog10", (polyRTSFunction)&PolyRealFLog10 }, { "PolyRealFSinh", (polyRTSFunction)&PolyRealFSinh }, { "PolyRealFCosh", (polyRTSFunction)&PolyRealFCosh }, { "PolyRealFTanh", (polyRTSFunction)&PolyRealFTanh }, { "PolyRealFAtan2", (polyRTSFunction)&PolyRealFAtan2 }, { "PolyRealFPow", (polyRTSFunction)&PolyRealFPow }, { "PolyRealFCopySign", (polyRTSFunction)&PolyRealFCopySign }, { "PolyRealFFloor", (polyRTSFunction)&PolyRealFFloor }, { "PolyRealFCeil", (polyRTSFunction)&PolyRealFCeil }, { "PolyRealFTrunc", (polyRTSFunction)&PolyRealFTrunc }, { "PolyRealFRound", (polyRTSFunction)&PolyRealFRound }, { "PolyRealFRem", (polyRTSFunction)&PolyRealFRem }, { "PolyRealFNextAfter", (polyRTSFunction)&PolyRealFNextAfter }, { NULL, NULL} // End of list. }; class RealArithmetic: public RtsModule { public: virtual void Init(void); }; // Declare this. It will be automatically added to the table. static RealArithmetic realModule; void RealArithmetic::Init(void) { /* Some compilers object to overflow in constants so we compute the values here. */ #if (HAVE_DECL_FPSETMASK && ! defined(__CYGWIN__)) /* In FreeBSD 3.4 at least, we sometimes get floating point exceptions if we don't clear the mask. Maybe need to do this on other platforms as well just to be sure. */ // N.B. fpsetmask is defined in the headers on Cygwin but there's no function! fpsetmask(0); #endif // NAN and INFINITY are defined in GCC but not in Visual C++. #if (defined(INFINITY)) posInf = INFINITY; negInf = -(INFINITY); posInfF = INFINITY; negInfF = -(INFINITY); #else { double zero = 0.0; posInf = 1.0 / zero; negInf = -1.0 / zero; float zeroF = 0.0; posInfF = 1.0 / zeroF; negInfF = -1.0 / zeroF; } #endif #if (defined(NAN)) notANumber = NAN; #else { double zero = 0.0; notANumber = zero / zero; float zeroF = 0.0; notANumberF = zeroF / zeroF; } #endif // Make sure this is a positive NaN since we return it from "abs". // "Positive" in this context is copysign(1.0, x) > 0.0 because that's // how we test the sign so we test it first and then try to change the // sign if it's wrong. if (copysign(1.0, notANumber) < 0) notANumber = copysign(notANumber, 1.0); if (copysignf(1.0, notANumberF) < 0) notANumberF = copysignf(notANumberF, 1.0); }