[PATCH] Fraction and Integer asFloat: answer nearest floating point
Locked
[PATCH] Fraction and Integer asFloat: answer nearest floating point
 
|
Hello,
here are some enhancements that should answer nearest floating point value when converting LargeInteger and Fraction asFloat. For LargePositiveInteger, i have two versions, one un-optimized and one optimized for gst. All this has already been written for Squeak VW ST/X and Dolphin. Some possible tests following (to be completed from above dialects already written tests) Nicolas self assert: ((1 bitShift: 52)+0+(1/4)) asFloatD asExactFraction = ((1 bitShift: 52)+0). self assert: ((1 bitShift: 52)+0+(1/2)) asFloatD asExactFraction = ((1 bitShift: 52)+0). self assert: ((1 bitShift: 52)+0+(3/4)) asFloatD asExactFraction = ((1 bitShift: 52)+1). self assert: ((1 bitShift: 52)+1+(1/4)) asFloatD asExactFraction = ((1 bitShift: 52)+1). self assert: ((1 bitShift: 52)+1+(1/2)) asFloatD asExactFraction = ((1 bitShift: 52)+2). self assert: ((1 bitShift: 52)+1+(3/4)) asFloatD asExactFraction = ((1 bitShift: 52)+2). self assert: ((1 bitShift: 23)+0+(1/4)) asFloatE asExactFraction = ((1 bitShift: 23)+0). self assert: ((1 bitShift: 23)+0+(1/2)) asFloatE asExactFraction = ((1 bitShift: 23)+0). self assert: ((1 bitShift: 23)+0+(3/4)) asFloatE asExactFraction = ((1 bitShift: 23)+1). self assert: ((1 bitShift: 23)+1+(1/4)) asFloatE asExactFraction = ((1 bitShift: 23)+1). self assert: ((1 bitShift: 23)+1+(1/2)) asFloatE asExactFraction = ((1 bitShift: 23)+2). self assert: ((1 bitShift: 23)+1+(3/4)) asFloatE asExactFraction = ((1 bitShift: 23)+2). self assert: ((-1075 to: 1023) allSatisfy: [:i | (1.0 timesTwoPower: i) asExactFraction asFloatD = (1.0 timesTwoPower: i)]). "Note that min denormalized is (1.0 timesTwoPower: -1074)" self assert: ((1 to: 1075) allSatisfy: [:i | (1 bitShift: i) reciprocal negated asFloatD = (1.0 timesTwoPower: i negated) negated]). "Filed out from GNU Smalltalk version 2.3.1 on 27-Jan-2007 23:29:25"! !Fraction methodsFor: 'private'! asFloat: characterization "Answer the receiver converted to a Float" | n d sign hn hd hq nBits q q1 r exponent floatExponent | sign := numerator sign * denominator sign. n := numerator abs. d := denominator abs. hn := n highBit. hd := d highBit. "If both numerator and denominator are represented exactly in floating point number, then fastest thing to do is to use hardwired float division" nBits := characterization precision + 1. (hn < nBits and: [hd < nBits]) ifTrue: [^(characterization coerce: numerator) / (characterization coerce: denominator)]. "Try and obtain a mantissa with characterization precision + 1 bits by integer division. Additional bit is a helper for rounding mode. First guess is rough, we might get one more bit or one less" exponent := hn - hd - nBits. exponent > 0 ifTrue: [d := d bitShift: exponent] ifFalse: [n := n bitShift: exponent negated]. q := n quo: d. r := n - (q * d). hq := q highBit. "check for gradual underflow, in which case we should use less bits" floatExponent := exponent + hq. floatExponent >= (characterization emin - 1) ifFalse: [nBits := nBits + floatExponent - characterization emin+1]. "Use exactly nBits" hq > nBits ifTrue: [exponent := exponent + hq - nBits. r := (q bitAnd: (1 bitShift: hq - nBits) - 1) * d + r. q := q bitShift: nBits - hq]. hq < nBits ifTrue: [exponent := exponent + hq - nBits. q1 := (r bitShift: nBits - hq) quo: d. q := (q bitShift: nBits - hq) bitAnd: q1. r := (r bitShift: nBits - hq) - (q1 * d)]. "check if we should round upward. The case of exact half (q bitAnd: 1) = 1 & (r = 0) will be handled by Integer>>asFloat:" ((q bitAnd: 1) = 0 or: [r = 0]) ifFalse: [q := q + 1]. "build the Float" ^(sign > 0 ifTrue: [characterization coerce: q] ifFalse: [(characterization coerce: q) negated]) timesTwoPower: exponent! ! "Filed out from GNU Smalltalk version 2.3.1 on 28-Jan-2007 0:37:07"! !LargePositiveInteger methodsFor: 'private'! asFloat: characterization "Answer the receiver converted to a Float" | nTruncatedBits result byte | "check for number bigger than maximum mantissa" nTruncatedBits := self highBit - characterization precision. nTruncatedBits > 0 ifTrue: [| mantissa exponent mask trailingBits inexact carry | mantissa := self bitShift: nTruncatedBits negated. exponent := nTruncatedBits. mask := (1 bitShift: nTruncatedBits) - 1. trailingBits := self bitAnd: mask. inexact := trailingBits ~= 0. inexact ifTrue: ["Apply IEEE 754 round to nearest even default rounding mode" carry := self bitAt: nTruncatedBits. (carry = 0 or: [(trailingBits bitAnd: (mask bitShift: -1)) = 0 and: [mantissa even]]) ifFalse: [mantissa := mantissa + 1]]. ^(characterization coerce: mantissa) timesTwoPower: exponent]. "conversion can be exact, construct Float by successive mul add operations" byte := characterization coerce: 256. result := characterization coerce: 0. self size to: 1 by: -1 do: [ :index | result := result * byte + (self at: index). ]. ^result! ! "Filed out from GNU Smalltalk version 2.3.1 on 28-Jan-2007 1:04:10"! !LargePositiveInteger methodsFor: 'private'! asFloat: characterization "Answer the receiver converted to a Float" | nTruncatedBits result byte | "check for number bigger than maximum mantissa" nTruncatedBits := self highBit - characterization precision. nTruncatedBits > 0 ifTrue: [| mantissa exponent carry | mantissa := self bitShift: nTruncatedBits negated. exponent := nTruncatedBits. "Apply IEEE 754 round to nearest even default rounding mode" "inexactFlag := (self checkIfLowBitGreaterThan: nTruncatedBits) not" carry := self bitAt: nTruncatedBits. (carry = 1 and: [mantissa odd or: [(self checkIfLowBitGreaterThan: nTruncatedBits - 1) not]]) ifTrue: [mantissa := mantissa + 1]. ^(characterization coerce: mantissa) timesTwoPower: exponent]. "conversion can be exact, construct Float by successive mul add operations" byte := characterization coerce: 256. result := characterization coerce: 0. self size to: 1 by: -1 do: [ :index | result := result * byte + (self at: index). ]. ^result! ! "Filed out from GNU Smalltalk version 2.3.1 on 28-Jan-2007 1:12:40"! !LargePositiveInteger methodsFor: 'private'! checkIfLowBitGreaterThan: n "Answer true if all bit of self lesser or equal to n are zero" "fast check byte by byte" 1 to: n // 8 do: [:i | (self digitAt: i) = 0 ifFalse: [^false]]. "finish checking bit by bit" (n // 8) * 8 + 1 to: n do: [:i | (self bitAt: i) = 0 ifFalse: [^false]]. "all low bits were zero" ^true ! ! _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
Some more tests:
"check for negative zero" self assert: (1 bitShift: 1075) reciprocal asFloatD positive. self assert: (1 bitShift: 1075) reciprocal negated asFloatD negative. "check for infinity" self assert: (1 bitShift: 1024) asFloatD = FloatD infinity. self assert: (1 bitShift: 1024) negated asFloatD = FloatD negativeInfinity. self assert: ((1 bitShift: 1075)/3) asFloatD = FloatD infinity. self assert: ((1 bitShift: 1075)/3) negated asFloatD = FloatD negativeInfinity. "check for non infinity/nan" self assert: ((1 bitShift: 1024)/3) asFloatD isFinite. self assert: ((1 bitShift: 1024)/3) negated asFloatD isFinite. self assert: (3/(1 bitShift: 1075)) asFloatD isFinite. self assert: (3/(1 bitShift: 1075)) negated asFloatD isFinite. _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
In reply to this post by Nicolas Cellier-3
nicolas cellier wrote:
> Hello, > > here are some enhancements that should answer nearest floating point > value when converting LargeInteger and Fraction asFloat. > > For LargePositiveInteger, i have two versions, one un-optimized and one > optimized for gst. Which is the optimized one (the one whose source does not mention checkIfLowBitGreaterThan: n, or the other)? Paolo _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
In reply to this post by Nicolas Cellier-3
Thanks, I applied a pretty big patch with all your changes and mine.
Paolo 2007-01-27 Paolo Bonzini <[hidden email]> Nicolas Cellier <[hidden email]> * kernel/Float.st: Rewrite #truncated and #asExactFraction. * kernel/Fraction.st: Use #asFloatD in #hash. Rewrite #asFloat:. * kernel/LargeInt.st: Rewrite #asFloat:. Add #trailingZeros. * kernel/SmallInt.st: Add #trailingZeros. * kernel/Integer.st: Add #trailingZeros. * kernel/MappedColl.st: Use map size as collection size. * tests/floatmath.st: Add accuracy tests. * tests/floatmath.ok: Adjust. 2007-01-27 Paolo Bonzini <[hidden email]> * lib-src/truncl.c: New. 2007-01-28 Paolo Bonzini <[hidden email]> * libgst/prims.def: Use truncl and lrint to implement conversion from float to integer. --- orig/NEWS +++ mod/NEWS @@ -27,6 +27,9 @@ o #copyFrom:to: is uniformly 0-based f o Fixed a garbage collection bug that typically occurred when installing GNU Smalltalk, or when launching the installed image. +o Fixed many floating point rounding bugs in LargeIntegers and Fractions, + thanks to Nicolas Cellier. + o gst-package honors the INSTALL command found by configure. o gst-config does not "forget" to prefix the library directories with -L. --- orig/configure.ac +++ mod/configure.ac @@ -204,7 +204,7 @@ AC_CHECK_LIB(m, atan) GST_REPLACE_POLL AC_REPLACE_FUNCS(putenv strdup strerror strsignal mkstemp getpagesize \ getdtablesize strstr ftruncate floorl ceill sqrtl frexpl ldexpl asinl \ - acosl atanl logl expl tanl sinl cosl strsep strpbrk) + acosl atanl logl expl tanl sinl cosl truncl strsep strpbrk) AC_CHECK_FUNCS_ONCE(gethostname memcpy memmove sighold uname sbrk usleep lstat \ grantpt popen getrusage gettimeofday getcwd fork strchr \ sigsetmask alarm select mprotect madvise nl_langinfo waitpid \ --- orig/kernel/Float.st +++ mod/kernel/Float.st @@ -194,27 +194,30 @@ truncated ifFalse: [ float := self negated ]. exponent := float exponent. - bytes := ByteArray new: exponent // 8 + 2. - float := float timesTwoPower: (exponent bitClear: 7) negated. + bytes := LargePositiveInteger new: (self class precision + 7) // 8 + 1. + float := float timesTwoPower: float class precision - exponent - 8. - bytes size - 1 to: 1 by: -1 do: [ :i | - bytes at: i put: float truncated. - float := float fractionPart timesTwoPower: 8 + 1 to: bytes size do: [ :i | + bytes digitAt: i put: (float fractionPart timesTwoPower: 8) truncated. + float := float integerPart timesTwoPower: -8. ]. - ^positive - ifTrue: [ (LargeInteger from: bytes) ] - ifFalse: [ (LargeInteger from: bytes) negated ] + bytes := bytes bitShift: (exponent - float class precision). + positive ifFalse: [ bytes := bytes negated ]. + ^bytes ! asExactFraction "Convert the receiver into a fraction with optimal approximation, but with usually huge terms." + | shift mantissa | self checkCoercion. - - ^(self timesTwoPower: self exponent + self class precision) truncated / - (1 bitShift: self exponent + self class precision) + shift := self exponent negated + self class precision. + mantissa := (self timesTwoPower: shift) truncated. + ^shift negative + ifTrue: [mantissa * (1 bitShift: shift negated)] + ifFalse: [mantissa / (1 bitShift: shift)] ! asFraction --- orig/kernel/Fraction.st +++ mod/kernel/Fraction.st @@ -288,7 +288,7 @@ truncated hash "Answer an hash value for the receiver" denominator = 1 ifTrue: [ ^numerator hash ]. - ^(numerator asFloatD / denominator asFloatD) hash + ^self asFloatD hash ! ! @@ -385,33 +385,59 @@ storeOn: aStream asFloat: characterization "Answer the receiver converted to a Float" - | n d shift sign | - n := numerator. - d := denominator. - sign := n sign * d sign. - - "Avoid answering NaNs and infinite values. - 1e1800 asFloat / (1e1799 + 1) asFloat = NaN, but - (1e1800 / (1e1799 + 1)) asFloat must be 10." - - shift := (characterization emax - numerator highBit). - shift := shift min: (characterization emax - denominator highBit). - shift < 0 ifTrue: [ - "Lose some more precision, but we MUST avoid infinites and NaNs!" - shift := shift - 10. n := n bitShift: shift. d := d bitShift: shift ]. - - d = 0 ifTrue: [ - ^sign > 0 - ifTrue: [ characterization infinity ] - ifFalse: [ characterization negativeInfinity ] ]. - n = 0 ifTrue: [ - ^sign > 0 - ifTrue: [ characterization coerce: 0 ] - ifFalse: [ characterization negativeInfinity reciprocal ] ]. + "Answer the receiver converted to a Float" - ^(characterization coerce: n) - / (characterization coerce: d) -! + | n d sign hn hd hq nBits q q1 r exponent floatExponent | + sign := numerator sign * denominator sign. + n := numerator abs. + d := denominator abs. + hn := n highBit. + hd := d highBit. + + "If both numerator and denominator are represented exactly in floating + point number, then fastest thing to do is to use hardwired float division" + nBits := characterization precision + 1. + (hn < nBits and: [hd < nBits]) + ifTrue: [^(characterization coerce: numerator) / (characterization coerce: denominator)]. + + "Try and obtain a mantissa with characterization precision + 1 bits by integer division. + Additional bit is a helper for rounding mode. + First guess is rough, we might get one more bit or one less" + exponent := hn - hd - nBits. + exponent > 0 + ifTrue: [d := d bitShift: exponent] + ifFalse: [n := n bitShift: exponent negated]. + q := n quo: d. + r := n - (q * d). + hq := q highBit. + + "check for gradual underflow, in which case we should use less bits" + floatExponent := exponent + hq. + floatExponent >= (characterization emin - 1) ifFalse: [nBits := nBits + floatExponent - characterization emin+1]. + + "Use exactly nBits" + hq > nBits + ifTrue: + [exponent := exponent + hq - nBits. + r := (q bitAnd: (1 bitShift: hq - nBits) - 1) * d + r. + q := q bitShift: nBits - hq]. + hq < nBits + ifTrue: + [exponent := exponent + hq - nBits. + q1 := (r bitShift: nBits - hq) quo: d. + q := (q bitShift: nBits - hq) bitAnd: q1. + r := (r bitShift: nBits - hq) - (q1 * d)]. + + "check if we should round upward. + The case of exact half (q bitAnd: 1) = 1 & (r = 0) + will be handled by Integer>>asFloat:" + ((q bitAnd: 1) = 0 or: [r = 0]) ifFalse: [q := q + 1]. + + "build the Float" + ^(sign > 0 + ifTrue: [characterization coerce: q] + ifFalse: [(characterization coerce: q) negated]) + timesTwoPower: exponent! reduce "Reduce the fraction." --- orig/kernel/Integer.st +++ mod/kernel/Integer.st @@ -139,12 +139,18 @@ clearBit: index ! noMask: anInteger - "True if no 1 bits in anInteger are 1 in the receiver" + "Answer true if no 1 bits in anInteger are 1 in the receiver." ^(self bitAnd: anInteger) = 0 ! +trailingZeros + "Return the index of the lowest order 1 bit of the receiver, + minus 1." + self subclassResponsibility +! + highBit - "Return the index of the highest order 1 bit of the receiver" + "Return the index of the highest order 1 bit of the receiver." self subclassResponsibility ! --- orig/kernel/LargeInt.st +++ mod/kernel/LargeInt.st @@ -493,6 +493,17 @@ negated !LargeInteger methodsFor: 'bit operations'! +trailingZeros + "Answer the number of trailing zero bits in the receiver" + + | each | + 1 to: self size do: [ :index | + (each := self digitAt: index) = 0 + ifFalse: [ ^index * 8 - 8 + (TrailingZeros at: each) ]. + ]. + ^self highBit - 1 +! + bitAnd: aNumber "Answer the receiver ANDed with aNumber" | newBytes | @@ -977,17 +988,32 @@ sign asFloat: characterization "Answer the receiver converted to a Float" - | adjust result byte | + | nTruncatedBits mantissa exponent mask trailingBits inexact carry | + + "Check for number bigger than maximum mantissa" + nTruncatedBits := self highBit - characterization precision. + nTruncatedBits <= 0 + ifTrue: [^self fastAsFloat: characterization ]. + + mantissa := self bitShift: nTruncatedBits negated. + exponent := nTruncatedBits. + + "Apply IEEE 754 round to nearest even default rounding mode" + carry := self bitAt: nTruncatedBits. + (carry = 1 and: [mantissa odd or: [self trailingZeros >= nTruncatedBits]]) + ifTrue: [mantissa := mantissa + 1]. + ^(characterization coerce: mantissa) timesTwoPower: exponent! + +fastAsFloat: characterization + "Conversion can be exact, construct Float by successive mul add operations" + | result byte | byte := characterization coerce: 256. - adjust := byte raisedToInteger: self size - 1. - result := 0. + result := characterization coerce: 0. self size to: 1 by: -1 do: [ :index | - result := (self at: index) * adjust + result. - adjust := adjust / byte. + result := result * byte + (self at: index). ]. - ^result -! + ^result! mostSignificantByte "Private - Answer the value of the most significant byte" --- orig/kernel/MappedColl.st +++ mod/kernel/MappedColl.st @@ -89,7 +89,7 @@ at: key put: value size "Answer the receiver's size" - ^domain size + ^map size ! add: anObject --- orig/kernel/SmallInt.st +++ mod/kernel/SmallInt.st @@ -94,6 +94,24 @@ generality !SmallInteger methodsFor: 'bit arithmetic'! +trailingZeros + "Return the index of the trailing zero bits in the receiver" + + | n bit | + self = 0 ifTrue: [ ^0 ]. + + n := self. + bit := 0. + + (n bitAnd: 16r3FFFFFFF) = 0 ifTrue: [ + bit := bit + 30. n := n bitShift: -30 ]. + + (n bitAnd: 16rFFFF) = 0 ifTrue: [ bit := bit + 16. n := n bitShift: -16 ]. + (n bitAnd: 16rFF) = 0 ifTrue: [ bit := bit + 8. n := n bitShift: -8 ]. + (n bitAnd: 16rF) = 0 ifTrue: [ bit := bit + 4. n := n bitShift: -4 ]. + (n bitAnd: 16r3) = 0 ifTrue: [ bit := bit + 2. n := n bitShift: -2 ]. + ^bit + (1 - (n bitAnd: 16r1))! + highBit "Return the index of the highest order 1 bit of the receiver" --- orig/libgst/prims.def +++ mod/libgst/prims.def @@ -1531,14 +1531,9 @@ primitive VMpr_FloatD_truncated [succeed if COMMON (RECEIVER_IS_CLASS (oop1, _gst_floatd_class)) { double oopValue = FLOATD_OOP_VALUE (oop1); - if COMMON (oopValue >= 0.0 && oopValue <= MAX_ST_INT) + if COMMON (oopValue >= MIN_ST_INT && oopValue <= MAX_ST_INT) { - PUSH_INT ((intptr_t) (oopValue + 0.000000000000005)); - PRIM_SUCCEEDED; - } - else if COMMON (oopValue < 0.0 && oopValue >= MIN_ST_INT) - { - PUSH_INT ((intptr_t) (oopValue - 0.000000000000005)); + PUSH_INT (lrint (trunc (oopValue))); PRIM_SUCCEEDED; } } @@ -1732,14 +1727,9 @@ primitive VMpr_FloatE_truncated [succeed if COMMON (RECEIVER_IS_CLASS (oop1, _gst_floate_class)) { double oopValue = FLOATE_OOP_VALUE (oop1); - if COMMON (oopValue >= 0.0 && oopValue <= MAX_ST_INT) + if COMMON (oopValue >= MIN_ST_INT && oopValue <= MAX_ST_INT) { - PUSH_INT ((intptr_t) (oopValue + 0.000000000000005)); - PRIM_SUCCEEDED; - } - else if COMMON (oopValue < 0.0 && oopValue >= MIN_ST_INT) - { - PUSH_INT ((intptr_t) (oopValue - 0.000000000000005)); + PUSH_INT (lrintf (truncf (oopValue))); PRIM_SUCCEEDED; } } @@ -1933,14 +1923,9 @@ primitive VMpr_FloatQ_truncated [succeed if COMMON (RECEIVER_IS_CLASS (oop1, _gst_floatq_class)) { long double oopValue = FLOATQ_OOP_VALUE (oop1); - if COMMON (oopValue >= 0.0 && oopValue <= MAX_ST_INT) - { - PUSH_INT ((intptr_t) (oopValue + 0.000000000000005)); - PRIM_SUCCEEDED; - } - else if COMMON (oopValue < 0.0 && oopValue >= MIN_ST_INT) + if COMMON (oopValue >= MIN_ST_INT && oopValue <= MAX_ST_INT) { - PUSH_INT ((intptr_t) (oopValue - 0.000000000000005)); + PUSH_INT (lrintl (truncl (oopValue))); PRIM_SUCCEEDED; } } --- orig/tests/floatmath.ok +++ mod/tests/floatmath.ok @@ -1,248 +1,74 @@ - -Execution begins... -returned value is 3.10000 - -Execution begins... -returned value is 3.45000 - -Execution begins... -returned value is 30000.0 - -Execution begins... -returned value is 34500.0 - -Execution begins... -returned value is 7.70000 - -Execution begins... -returned value is -8.62000 - -Execution begins... -returned value is false - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is false - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is 0.00180000 - -Execution begins... -returned value is 11250.0 - -Execution begins... -returned value is 3 - -Execution begins... -returned value is 0.141593 - -Execution begins... -returned value is 12 - -Execution begins... -returned value is 720.000 - -Execution begins... -returned value is 2.81250 - -Execution begins... -returned value is inf - -Execution begins... -returned value is 'Inf' - -Execution begins... -returned value is -inf - -Execution begins... -returned value is '-Inf' - -Execution begins... -returned value is nan - -Execution begins... -returned value is 'NaN' - -Execution begins... -returned value is '0.0' - -Execution begins... -returned value is '-0.0' - -Execution begins... -returned value is true - -Execution begins... -returned value is false - -Execution begins... -returned value is nan - -Execution begins... -returned value is nan - -Execution begins... -returned value is nan - -Execution begins... -returned value is 5.00000 - -Execution begins... -returned value is 5.00000 - -Execution begins... -returned value is 7.00000 - -Execution begins... -returned value is 0.00000 - -Execution begins... -returned value is -0.00000 - -Execution begins... -returned value is -0.00000 - -Execution begins... -returned value is -0.00000 - -Execution begins... -returned value is nan - -Execution begins... -returned value is nan - -Execution begins... -returned value is nan - -Execution begins... -returned value is 7.00000 - -Execution begins... -returned value is 7.00000 - -Execution begins... -returned value is 7.00000 - -Execution begins... -returned value is 0.00000 - -Execution begins... -returned value is 0.00000 - -Execution begins... -returned value is 0.00000 - -Execution begins... -returned value is -0.00000 - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -returned value is true - -Execution begins... -(-0.0 -0.0 0.0 0.0 true 0.0 0.0 true ) -(-0.0 0.0 0.0 0.0 true -0.0 -0.0 true ) -(0.0 -0.0 -0.0 -0.0 true 0.0 0.0 true ) -(0.0 0.0 0.0 0.0 true 0.0 0.0 true ) -returned value is Array new: 4 "<0>" - -Execution begins... -(-0.0 -0.0 0.0 0.0 true 0.0 0.0 true ) -(-0.0 0.0 0.0 0.0 true -0.0 -0.0 true ) -(0.0 -0.0 -0.0 -0.0 true 0.0 0.0 true ) -(0.0 0.0 0.0 0.0 true 0.0 0.0 true ) -returned value is Array new: 4 "<0>" - -Execution begins... -(-0.0 -0.0 0.0 0.0 true 0.0 0.0 true ) -(0.0 0.0 0.0 0.0 true 0.0 0.0 true ) -returned value is Array new: 2 "<0>" - -Execution begins... -(-0.0 0.0 true true true ) -(0.0 -0.0 true true true ) -returned value is Array new: 2 "<0>" - -Execution begins... -returned value is Float - -Execution begins... -true->1.0e16 -returned value is 1.00000e+16 - -Execution begins... -true->1.2345e16 -returned value is 1.23450e+16 - -Execution begins... -true->10.0 -returned value is 10.0000 - -Execution begins... -true->17.7674749 -returned value is 17.7675 - -Execution begins... -true->0.12345 -returned value is 0.123450 - -Execution begins... -true->0.0000000012345 -returned value is 1.23450d-09 - -Execution begins... -true->1.2345e-9 -returned value is 1.23450e-09 - -Execution begins... -true->0.83205029433784 -returned value is 0.832050 - -Execution begins... -true->0.832050294337844 -returned value is 0.832050 - -Execution begins... -true->0.55470019622523 -returned value is 0.554700 - -Execution begins... -true->0.554700196225229 -returned value is 0.554700 +*** floatmath.ok Wed Nov 8 11:54:47 2006 +--- /Volumes/disk0s8/devel/gst/+build/tests/floatmath.log Sun Jan 28 22:31:28 2007 +*************** +*** 204,213 **** +--- 204,221 ---- + returned value is Float + + Execution begins... ++ true->10000000000000000.0 ++ returned value is 1.00000d+16 ++ ++ Execution begins... + true->1.0e16 + returned value is 1.00000e+16 + + Execution begins... ++ true->12345000000000000.0 ++ returned value is 1.23450d+16 ++ ++ Execution begins... + true->1.2345e16 + returned value is 1.23450e+16 + +*************** +*** 246,248 **** +--- 254,301 ---- + Execution begins... + true->0.554700196225229 + returned value is 0.554700 ++ ++ Execution begins... ++ returned value is Float class ++ ++ Execution begins... ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ returned value is FloatD ++ ++ Execution begins... ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ true ++ returned value is FloatE --- orig/tests/floatmath.st +++ mod/tests/floatmath.st @@ -8,7 +8,7 @@ "====================================================================== | -| Copyright (C) 1988, 1989, 1999, 2006 Free Software Foundation. +| Copyright (C) 1988, 1989, 1999, 2006, 2007 Free Software Foundation. | Written by Steve Byrne | | This file is part of GNU Smalltalk. @@ -168,7 +168,9 @@ test (((Behavior evaluate: self printString) = self) -> self) printNl. ! ! +1d16 test! 1e16 test! +1.2345d16 test! 1.2345e16 test! 10.0 test! (20 - 2.2325251) test! @@ -179,3 +181,55 @@ test 0.832050294337844 test! 0.55470019622523 test! 0.554700196225229 test! + + +"Fun with rounding" + +!Float class methodsFor: 'testing'! + +assert: aBoolean + aBoolean ifFalse: [ self halt ] ifTrue: [ aBoolean printNl ]! + +test + | p | + p := 1 bitShift: self precision - 1. + self assert: (self coerce: p+0+(1/4)) asExactFraction = (p+0). + self assert: (self coerce: p+0+(1/2)) asExactFraction = (p+0). + self assert: (self coerce: p+0+(3/4)) asExactFraction = (p+1). + self assert: (self coerce: p+1+(1/4)) asExactFraction = (p+1). + self assert: (self coerce: p+1+(1/2)) asExactFraction = (p+2). + self assert: (self coerce: p+1+(3/4)) asExactFraction = (p+2). + + self assert: ((self emin - self precision - 1 to: self emax - 1) allSatisfy: [:i | + p := (self coerce: 1) timesTwoPower: i. + (self coerce: p asExactFraction) = p]). + + self assert: ((1 to: 1 + self precision - self emin) allSatisfy: [:i | + p := (self coerce: 1) timesTwoPower: i negated. + (self coerce: (1 bitShift: i) reciprocal negated) = p negated]). + + "check for negative zero" + p := 1 bitShift: 1 + self precision - self emin. + self assert: (self coerce: p reciprocal) positive. + self assert: (self coerce: p reciprocal negated) negative. + + "check for infinity" + p := 1 bitShift: self emax + 1. + self assert: (self coerce: p) = self infinity. + self assert: (self coerce: p negated) = self negativeInfinity. + + p := 1 bitShift: 1 + self precision - self emin. + self assert: (self coerce: p / 3) = self infinity. + self assert: (self coerce: p / -3) = self negativeInfinity. + + "check for non infinity/nan" + p := 1 bitShift: self emax + 1. + self assert: (self coerce: p / 3) isFinite. + self assert: (self coerce: p / -3) isFinite. + + p := 1 bitShift: 1 + self precision - self emin. + self assert: (self coerce: 3 / p) isFinite. + self assert: (self coerce: -3 / p) isFinite! ! + +FloatD test! +FloatE test! * added files --- /dev/null +++ /Volumes/disk0s8/devel/gst/,,[hidden email]--2004b/new-files-archive/./lib-src/.arch-ids/truncl.c.id @@ -0,0 +1 @@ +Paolo Bonzini <[hidden email]> Sun Jan 28 22:17:40 2007 18819.0 --- /dev/null +++ /Volumes/disk0s8/devel/gst/,,[hidden email]--2004b/new-files-archive/./lib-src/truncl.c @@ -0,0 +1,79 @@ +/******************************** -*- C -*- **************************** + * + * Emulation for truncl + * + * + ***********************************************************************/ + +/*********************************************************************** + * + * Copyright 2007 Free Software Foundation, Inc. + * Written by Paolo Bonzini. + * + * This file is part of GNU Smalltalk. + * + * GNU Smalltalk is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 2, or (at your option) any later + * version. + * + * Linking GNU Smalltalk statically or dynamically with other modules is + * making a combined work based on GNU Smalltalk. Thus, the terms and + * conditions of the GNU General Public License cover the whole + * combination. + * + * In addition, as a special exception, the Free Software Foundation + * give you permission to combine GNU Smalltalk with free software + * programs or libraries that are released under the GNU LGPL and with + * independent programs running under the GNU Smalltalk virtual machine. + * + * You may copy and distribute such a system following the terms of the + * GNU GPL for GNU Smalltalk and the licenses of the other code + * concerned, provided that you include the source code of that other + * code when and as the GNU GPL requires distribution of source code. + * + * Note that people who make modified versions of GNU Smalltalk are not + * obligated to grant this special exception for their modified + * versions; it is their choice whether to do so. The GNU General + * Public License gives permission to release a modified version without + * this exception; this exception also makes it possible to release a + * modified version which carries forward this exception. + * + * GNU Smalltalk 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 General Public License for + * more details. + * + * You should have received a copy of the GNU General Public License along with + * GNU Smalltalk; see the file COPYING. If not, write to the Free Software + * Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + ***********************************************************************/ + +#include <float.h> + +#include "mathl.h" + +/* To compute the integer part of X, sum a big enough + integer so that the precision of the floating point + number is exactly 1. */ + +long double +truncl(long double x) +{ + long double y; + if (x < 0.0L) + { + y = -(1.0L / LDBL_EPSILON - x - 1.0 / LDBL_EPSILON); + if (y < x) + y = y + 1.0L; + } + else + { + y = 1.0L / LDBL_EPSILON + x - 1.0 / LDBL_EPSILON; + if (y > x) + y = y - 1.0L; + } + + return y; +} _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
In reply to this post by Paolo Bonzini
Paolo Bonzini a écrit :
> nicolas cellier wrote: >> Hello, >> >> here are some enhancements that should answer nearest floating point >> value when converting LargeInteger and Fraction asFloat. >> >> For LargePositiveInteger, i have two versions, one un-optimized and >> one optimized for gst. > > Which is the optimized one (the one whose source does not mention > checkIfLowBitGreaterThan: n, or the other)? > > Paolo The optimized one is in file named 'LargePositiveIntegerOptimized.asFloat:.st' and it calls #checkIfLowBitGreaterThan: This last message is most of the optimization. We just need to know if some trailing bits are truncated or not in order to do the rounding. Non optimized version is using bit operations with masks like traditional Smalltalk patterns (bitShift: and bitAnd:). This is efficient for SmallIntegers. Unfortunately, in our case this creates a lot of intermediary LargePositiveInteger objects and is slow. Optimized version do this without creating any new LargePositiveInteger, just iterating on bits with bitAt: test, which is further accelerated by testing bytes with digitAt: tests when number of bits > 8 (and would be further optimized in a primitive with 16, 32 or 64 bits by 64 bits tests...). Also, non zero bit search is stopped as soon as a non zero bit is found. I am quite sure that a more general message named #lowBit would do almost as well and could be used in other algorithms. LargePositiveInteger>>lowBit "answer the index of least significant non zero bit. answer 0 if receiver is zero. LSB is indexed 1." 1 to: self size do: [:i | | lowBit | lowBit := (self digitAt: i) lowBit. lowBit = 0 ifFalse: [(î-1)*8 + lowBit]]. ^0 SmallInteger>>lowBit "same as above. Note: there are faster implementations to be found for sure" | shifted lowBit | self = 0 ifTrue: [^0]. self < 0 ifTrue: [self error: 'i do not know how to handle this one... and i do not much like negative integer highBit implementation neither']. shifted := self. lowBit := 1. [(shifted bitAnd: 1) = 0] whileTrue: [shifted := shifted bitShift: -1. lowBit := lowBit+1}. ^lowBit Nicolas _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
Heh, I had just implemented #trailingZeros while working on your reports
yesterday (but then I didn't use it). So I had actually just recognized that your #checkIfLowBitGreaterThan: did the same as my #trailingZeros, and understood that it was the optimized one. Paolo _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
Paolo Bonzini a écrit :
> Heh, I had just implemented #trailingZeros while working on your reports > yesterday (but then I didn't use it). So I had actually just recognized > that your #checkIfLowBitGreaterThan: did the same as my #trailingZeros, > and understood that it was the optimized one. > > Paolo Yes, our posts did cross... your trailing zero must be "my" lowBit - 1, isn't it? I called it lowBit in reference to highBit. And i should not say "my" because I am sure there is already a lowBit implemented in Squeak, probably in other dialects too. _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
Re: Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point
|
|
> Yes, our posts did cross... your trailing zero must be "my" lowBit - 1, > isn't it? > I called it lowBit in reference to highBit. > And i should not say "my" because I am sure there is already a lowBit > implemented in Squeak, probably in other dialects too. Ok, I changed my trailingZeros to lowBit, found a bug in detecting the number of digits necessary to print a number, fixed the lex.c bug with denormals, and made Float>>#raisedToInteger: more accurate by computing the reciprocal (if necessary) at the end of the computation rather than at the beginning. :-) Paolo 2007-01-29 Paolo Bonzini <[hidden email]> * kernel/LargeInt.st: Rename #trailingZeros to #lowBit. * kernel/SmallInt.st: Rename #trailingZeros to #lowBit. * kernel/Integer.st: Rename #trailingZeros to #lowBit. * kernel/Float.st: Fix bug in printing exact floating-point numbers. Support gradual underflow in #raisedToInteger:. * tests/floatmath.st: Add accuracy tests. * tests/floatmath.ok: Adjust. * libgst/lex.c: Rename ipowl to mul_powl, support gradual underflow. --- orig/kernel/Float.st +++ mod/kernel/Float.st @@ -132,6 +132,41 @@ negated integerPart "Return the receiver's integer part" ^self - self fractionPart +! + +raisedToInteger: anInteger + "Return self raised to the anInteger-th power" + + | exp adjustExp val mant | + + "Some special cases first" + anInteger isInteger ifFalse: [ SystemExceptions.WrongClass signalOn: anInteger mustBe: Integer ]. + anInteger = 0 ifTrue: [ ^self unity ]. + anInteger = 1 ifTrue: [ ^self ]. + + "Avoid overflow when the result is denormal and we would have an + unrepresentable intermediate result for its reciprocal." + adjustExp := self exponent. + exp := anInteger abs. + (anInteger > 0 or: [ (adjustExp + 1) * exp < self class emax ]) + ifTrue: [ + mant := self. + adjustExp := 0 ] + ifFalse: [ + mant := self timesTwoPower: 0 - adjustExp. + adjustExp := adjustExp * anInteger ]. + + "Fire the big loop." + val := mant + raisedToInteger: exp + withCache: ((Array new: (255 min: exp)) + at: 1 put: mant; + yourself). + +(mant->exp->val->adjustExp) printNl. + anInteger < 0 ifTrue: [ val := val reciprocal ]. + adjustExp = 0 ifFalse: [ val := val timesTwoPower: adjustExp ]. + ^val ! ! @@ -371,7 +406,7 @@ printOn: aStream special: whatToPrintArr "Private - Print a decimal representation of the receiver on aStream, printing one of the three elements of whatToPrintArray if it is infinity, negative infinity, or a NaN" - | me exponential small num den gcd + | me exponential small num numLog den denLog gcd intFactor precision int rounding digits digitStream exponent dotPrinted | @@ -415,11 +450,17 @@ printOn: aStream special: whatToPrintArr "Get the first `me class decimalDigits' base-10 digits of num // den, appropriately rounded" - intFactor := 10 raisedToInteger: (den ceilingLog: 10). - rounding := (10 raisedToInteger: (num ceilingLog: 10) - - me class decimalDigits) + 1. + numLog := num ceilingLog: 10. + denLog := den ceilingLog: 10. + numLog > me class decimalDigits + ifTrue: [ + intFactor := 10 raisedToInteger: denLog. + rounding := 10 raisedToInteger: numLog - me class decimalDigits ] + ifFalse: [ + intFactor := 10 raisedToInteger: me class decimalDigits - denLog. + rounding := 0 ]. - int := ((num * intFactor) + ((den // 2) * rounding)) // den. + int := ((num * intFactor) + ((den // 2) * (rounding + 1))) // den. digits := int printString. digits size > me class decimalDigits ifTrue: [ digits := digits copyFrom: 1 to: me class decimalDigits ]. --- orig/kernel/Integer.st +++ mod/kernel/Integer.st @@ -143,9 +143,8 @@ noMask: anInteger ^(self bitAnd: anInteger) = 0 ! -trailingZeros - "Return the index of the lowest order 1 bit of the receiver, - minus 1." +lowBit + "Return the index of the lowest order 1 bit of the receiver." self subclassResponsibility ! --- orig/kernel/LargeInt.st +++ mod/kernel/LargeInt.st @@ -493,15 +493,15 @@ negated !LargeInteger methodsFor: 'bit operations'! -trailingZeros - - "Answer the number of trailing zero bits in the receiver" +lowBit + "Return the index of the lowest order 1 bit of the receiver." | each | 1 to: self size do: [ :index | (each := self digitAt: index) = 0 - ifFalse: [ ^index * 8 - 8 + (TrailingZeros at: each) ]. + ifFalse: [ ^index * 8 - 7 + (TrailingZeros at: each) ]. ]. - ^self highBit - 1 + ^self highBit ! bitAnd: aNumber @@ -1000,7 +1000,7 @@ asFloat: characterization "Apply IEEE 754 round to nearest even default rounding mode" carry := self bitAt: nTruncatedBits. - (carry = 1 and: [mantissa odd or: [self trailingZeros >= nTruncatedBits]]) + (carry = 1 and: [mantissa odd or: [self lowBit > nTruncatedBits]]) ifTrue: [mantissa := mantissa + 1]. ^(characterization coerce: mantissa) timesTwoPower: exponent! --- orig/kernel/SmallInt.st +++ mod/kernel/SmallInt.st @@ -94,14 +94,16 @@ generality !SmallInteger methodsFor: 'bit arithmetic'! -trailingZeros - "Return the index of the trailing zero bits in the receiver" +lowBit + "Return the index of the lowest order 1 bit of the receiver." | n bit | self = 0 ifTrue: [ ^0 ]. n := self. - bit := 0. + "The result is 1-based, but we start from 2 to compensate with the + subtraction in the final line." + bit := 2. (n bitAnd: 16r3FFFFFFF) = 0 ifTrue: [ bit := bit + 30. n := n bitShift: -30 ]. @@ -110,7 +112,7 @@ trailingZeros (n bitAnd: 16rFF) = 0 ifTrue: [ bit := bit + 8. n := n bitShift: -8 ]. (n bitAnd: 16rF) = 0 ifTrue: [ bit := bit + 4. n := n bitShift: -4 ]. (n bitAnd: 16r3) = 0 ifTrue: [ bit := bit + 2. n := n bitShift: -2 ]. - ^bit + (1 - (n bitAnd: 16r1))! + ^bit - (n bitAnd: 16r1)! highBit "Return the index of the highest order 1 bit of the receiver" --- orig/libgst/lex.c +++ mod/libgst/lex.c @@ -167,9 +167,10 @@ static int scan_symbol (int c, static int digit_to_int (int c, int base); -/* Raise BASE to the N-th power. */ -static inline long double ipowl (long double base, - int n); +/* Raise BASE to the N-th power and multiply MANT by the result. */ +static inline long double mul_powl (long double mant, + long double base, + int n); #ifdef LEXDEBUG static void print_token (int token, @@ -649,22 +650,33 @@ scan_ident (int c, long double -ipowl (long double base, - int n) +mul_powl (long double mant, + long double base, + int n) { - int k = 1; long double result = 1.0; - while (n) + int n_abs = n < 0 ? -n : n; + int exp; + int k = 1; + + /* Restrict base to the range from 1.0 to 2.0. */ + base = frexpl (base, &exp); + base *= 2; + exp--; + + while (n_abs) { - if (n & k) + if (n_abs & k) { result *= base; - n ^= k; + n_abs ^= k; } base *= base; k <<= 1; } - return result; + + mant = (n < 0) ? mant / result : mant * result; + return ldexpl (mant, exp * n); } int @@ -836,10 +848,8 @@ scan_number (int c, else _gst_unread_char (ic); - if (exponent > 0) - num *= ipowl ((long double) base, exponent); - else - num /= ipowl ((long double) base, -exponent); + if (exponent != 0) + num = mul_powl (num, (long double) base, exponent); if (isNegative) num = -num; --- orig/tests/floatmath.ok +++ mod/tests/floatmath.ok @@ -220,6 +220,10 @@ true->1.2345e16 returned value is 1.23450e+16 Execution begins... +true->1.25 +returned value is 1.25000 + +Execution begins... true->10.0 returned value is 10.0000 --- orig/tests/floatmath.st +++ mod/tests/floatmath.st @@ -172,6 +172,7 @@ test 1e16 test! 1.2345d16 test! 1.2345e16 test! +1.25 test! 10.0 test! (20 - 2.2325251) test! 0.12345 test! _______________________________________________ help-smalltalk mailing list [hidden email] http://lists.gnu.org/mailman/listinfo/help-smalltalk |
«
Return to Gnu
|
1 view|%1 views
| Free forum by Nabble | Edit this page |