[FIX] Float readFrom: accumulate round off errors

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

[FIX] Float readFrom: accumulate round off errors

Nicolas Cellier-3
Float>>readFrom: does not always answer same result as sscanf(),
atod(), ...

this is because the construct
  integerPart + fractionalPart asFloat * (10 raisedToInteger: exponent)
is doing successive rounding errors...

When using my correction of Fraction asFloat in previous post, only one
rounding to nearest float is done.

Here after, you will find a test case and a patch.

Implementation does not slow done string -> float conversion.
It even fast it up in case of large number of decimals...

Nicolas

"---------------------------------------------------------------------------------------------------------"
"---------------------------------------------------------------------------------------------------------"
"---------------------------------------------------------------------------------------------------------"

!FloatTest methodsFor!

testFloatPrintString
        "Debug reading/printing a Floating point number without accumulating
round off errors"

        | b r |
        b := ByteArray new: 8.
        r := RandomLinearCongruential
                newModulus: 16r100000000
                multiplier: 16r0F010800F
                increment: 16r00005F2ED.
        r seed: 1234567.
        100
                timesRepeat: [| f str |
                        b basicDwordAtOffset: 4 put: (r next; seed) - 1.
                        b basicDwordAtOffset: 0 put: (r next; seed) - 1.
                        ((b basicDwordAtOffset: 4) bitAnd: 16r7FF00000) = 16r7FF00000
                                ifFalse: ["avoid nan and infinity"
                                        f := b doubleAtOffset: 0.
                                        str := (String new: 64) writeStream.
                                        f printOn: str significantFigures: 17.
                                        self assert: (Float readFrom: str contents readStream) = f]].
        "test big num near infinity"
        10
                timesRepeat: [| f str |
                        b basicDwordAtOffset: 4 put: 16r7FE00000 + ((r next; seed) //
16r1000 - 1).
                        b basicDwordAtOffset: 0 put: (r next; seed) - 1.
                        f := b doubleAtOffset: 0.
                        str := (String new: 64) writeStream.
                        f printOn: str significantFigures: 17.
                        self assert: (Float readFrom: str contents readStream) = f].
        "test infinitesimal (gradual underflow)"
        10
                timesRepeat: [| f str |
                        b basicDwordAtOffset: 4 put: 0 + ((r next; seed) // 16r1000 - 1).
                        b basicDwordAtOffset: 0 put: (r next; seed) - 1.
                        f := b doubleAtOffset: 0.
                        str := (String new: 64) writeStream.
                        f printOn: str significantFigures: 17.
                        self assert: (Float readFrom: str contents readStream) = f].! !

!FloatTest categoriesFor: #testFloatPrintString!public!Testing! !

"---------------------------------------------------------------------------------------------------------"
"---------------------------------------------------------------------------------------------------------"
"---------------------------------------------------------------------------------------------------------"
"---------------------------------------------------------------------------------------------------------"
"---------------------------------------------------------------------------------------------------------"

!Number class methodsFor!

readSmalltalkRealFrom: aStream initialInteger: anInteger
        "Private - Answer a new, positive, <Float> or <ScaledDecimal>, read
from
        the <gettableStream>, aStream. The <integer>, integerPart, has already
been read from
        the stream and we are currently positioned immediately after the
decimal point."

        "Attempt to read positive fractional part"

        | nextChar precision fractionalPart start mantissa |
        start := aStream position.
        (fractionalPart := self readIntegerFrom: aStream radix: 10) isNil
                ifTrue:
                        ["Actually just an Integer with a trailing full stop, which we must
pop back onto the stream"

                        aStream pop.
                        ^anInteger].
        precision := aStream position - start.

        "2006/06/08 nice: use Euclidean division so as to not accumulate round
off errors
        we keep mantissa in integer form"
        fractionalPart = 0
                ifTrue:
                        [mantissa := anInteger.
                        precision := 0]
                ifFalse: [mantissa := anInteger * (10 raisedToInteger: precision) +
fractionalPart].

        "Process any exponent..."
        ((nextChar := aStream peek) == $e or: [nextChar == $d or: [nextChar ==
$q]])
                ifTrue:
                        [| exponent |
                        aStream next. "Skip the exponent character"
                        "Allow plus prefix on the exponent, although not ANSI Smalltalk
syntax"
                        (exponent := self readDecimalIntegerFrom: aStream allowPlus: true)
notNil
                                ifTrue: [^exponent >= precision
                                        ifTrue: [(mantissa * (10 raisedToInteger: exponent - precision))
asFloat]
                                        ifFalse: ["Trick: Fraction do not need to be normalized before
converted to Float
                                                Implied Euclidean division will be cheaper than full gcd
algorithm"
                                                (Fraction
                                                        numerator: mantissa
                                                        denominator: (10 raisedToInteger: exponent negated + precision))
asFloat]].

                        "Found Float with trailing exponent character which is not part of
the number, e.g. 1.5e
                        From ANSI standard p 28: 'An exponentLetter must be followed by an
explicit exponent'"
                        aStream pop]
                ifFalse:
                        [nextChar == $s
                                ifTrue:
                                        [aStream next. "Skip the $s"
                                        ^self
                                                readScaledDecimalFrom: aStream
                                                mantissa: mantissa / (10 raisedToInteger: precision)
                                                precision: precision]].

        "Normal Float, such as 2.5 (perhaps with trailing exponent character)"
        "Trick: Fraction do not need to be normalized before converted to
Float
        Implied Euclidean division will be cheaper than full gcd algorithm"
        ^precision = 0
                ifTrue: [mantissa asFloat]
                ifFalse: [(Fraction
                        numerator: mantissa
                        denominator: (10 raisedToInteger: precision)) asFloat]! !


Reply | Threaded
Open this post in threaded view
|

Re: [FIX] Float readFrom: accumulate round off errors

Blair McGlashan-4
"nicolas cellier" <[hidden email]> wrote in message
news:[hidden email]...
> Float>>readFrom: does not always answer same result as sscanf(),
> atod(), ...
>
> this is because the construct
>  integerPart + fractionalPart asFloat * (10 raisedToInteger: exponent)
> is doing successive rounding errors...
>

Thanks once again for your contribution (and test!). We'll include this in
the next minor release.

Regards

Blair