"{ Encoding: utf8 }"

"

 COPYRIGHT (c) 1994 by Claus Gittinger

	      All Rights Reserved

 This software is furnished under a license and may be used

 only in accordance with the terms of that license and with the

 inclusion of the above copyright notice.   This software may not

 be provided or otherwise made available to, or used by, any

 other person.  No title to or ownership of the software is

 hereby transferred.

"

"{ Package: 'stx:libbasic' }"

Number subclass:#LimitedPrecisionReal

	instanceVariableNames:''

	classVariableNames:''

	poolDictionaries:''

	category:'Magnitude-Numbers'

!

!LimitedPrecisionReal class methodsFor:'documentation'!

copyright

"

 COPYRIGHT (c) 1994 by Claus Gittinger

	      All Rights Reserved

 This software is furnished under a license and may be used

 only in accordance with the terms of that license and with the

 inclusion of the above copyright notice.   This software may not

 be provided or otherwise made available to, or used by, any

 other person.  No title to or ownership of the software is

 hereby transferred.

"

!

documentation

"

    Abstract superclass for any-precision floating point numbers (i.e. Float and Double).

    Due to historic reasons, ST/X's Floats are what Doubles are in ST-80.

    The reason is that in some Smalltalk's, double floats are called Float, and no single float exists (VSE, V'Age).

    This may change in the future.

    Range and Precision of Storage Formats:

      Format |   Class    |   Array Class   | Bits / Significant  | Smallest Pos Number | Largest Pos Number | Significant Digits 

             |            |                 |      (Binary)       |                     |                    |     (Decimal)

      -------+------------+-----------------+---------------------+---------------------+--------------------+--------------------    

      half   |     --     | HalfFloatArray  |    16 / 11          |  6.10... x 10−5     |  6.55...  x 10+5   |      3.3

      -------+------------+-----------------+---------------------+---------------------+--------------------+--------------------    

      single | ShortFloat | FloatArray      |    32 / 24          |  1.175...x 10-38    |  3.402... 10+38    |      6-9

      -------+------------+-----------------+---------------------+---------------------+--------------------+--------------------    

      double | Float      | DoubleArray     |    64 / 53          |  2.225...x 10-308   |  1.797... 10+308   |     15-17

      -------+------------+-----------------+---------------------+---------------------+--------------------+--------------------    

      double | LongFloat  |     --          |   128 / 113         |  3.362...x 10-4932  |  1.189... 10+4932  |     33-36

      ext    |            |                 |                     |                     |                    |

      (SPARC)|            |                 |                     |                     |                    |

      -------+            |                 |---------------------+---------------------+--------------------+--------------------    

      double |            |                 |    96 / 64          |  3.362...x 10-4932  |  1.189... 10+4932  |     18-21

      ext    |            |                 |                     |                     |                    |

      (x86)  |            |                 |                     |                     |                    |

      -------+------------+-----------------+---------------------+---------------------+--------------------+--------------------    

        --   | LargeFloat |     --          |     arbitrary       |  arbitrarily small  |  arbitrarily large |     arbitrary

      -------+------------+-----------------+---------------------+---------------------+--------------------+--------------------    

    HalfFloats are only supported in fixed array containers. This was added for OpenGL and other graphic librariess which allow

    for texture, and vertex data to be passed quickly in that format (see http://www.opengl.org/wiki/Small_Float_Formats).

    Long- and LargeFloat are not supported array containers. These formats are seldom used for buld data.

    [author:]

        Claus Gittinger

    [see also:]

        Fraction FixedPoint

"

! !

!LimitedPrecisionReal class methodsFor:'instance creation'!

coerce:aNumber

    ^ self subclassResponsibility.

!

fromInteger:anInteger

    "return a float with anInteger's value.

     Since floats have a limited precision, you usually loose bits when doing this

     (see Float decimalPrecision, LongFloat decimalPrecision."

    |newFloat sign absVal completeAbsVal delta mask|

    newFloat := self zero.

    sign := anInteger sign.

    (sign ~~ 0) ifTrue:[

        "be careful to do round to nearest even floating point value -

         thanks to Nicolas Cellier for this code"

        absVal := anInteger abs.

        delta := absVal highBit - self precision.

        delta > 0 ifTrue: [

            completeAbsVal := absVal.

            "eliminate insignificant/trailing bits"

            absVal := absVal bitShift: delta negated.

            "inexact := trailingBits ~= 0.

             Round to nearest even"

            (completeAbsVal bitAt:delta) ~= 0 ifTrue:[

                mask := (1 bitShift:delta-1) - 1.

                ((completeAbsVal bitAnd:mask) ~= 0 or:[absVal odd]) ifTrue:[

                    absVal := absVal + 1

                ].

            ].

        ] ifFalse: [

            delta := 0

        ].

        absVal digitLength to:1 by:-1 do:[:i |

            newFloat := (newFloat * 256) + (absVal digitByteAt:i).

        ].

        delta ~~ 0 ifTrue: [

            newFloat := newFloat timesTwoPower:delta.

        ].

        (sign < 0) ifTrue:[

            newFloat := newFloat negated

        ]

    ].

    ^ newFloat

    "

     ShortFloat fromInteger:2

     12345678901234567890 asShortFloat            

     1234567890 asFloat                     

     1234567890 asFloat asInteger                    

     -1234567890 asFloat asInteger                    

     12345678901234567890 asFloat storeString            

     12345678901234567890 asFloat asInteger   

     -12345678901234567890 asFloat asInteger   

     12345678901234567890 asLongFloat           

     12345678901234567890 asLongFloat asInteger 

     -12345678901234567890 asLongFloat asInteger 

     123456789012345678901234567890 asLongFloat           

     123456789012345678901234567890 asLongFloat asInteger  

     -123456789012345678901234567890 asLongFloat asInteger  

     1234567890123456789012345678901234567890 asLongFloat           

     1234567890123456789012345678901234567890 asLongFloat asInteger  

     -1234567890123456789012345678901234567890 asLongFloat asInteger

     'this test is on 65 bits'.

     self assert: 16r1FFFFFFFFFFFF0801 asDouble ~= 16r1FFFFFFFFFFFF0800 asDouble.

     'this test is on 64 bits'.

     self assert: 16r1FFFFFFFFFFFF0802 asDouble ~= 16r1FFFFFFFFFFFF0800 asDouble.

     'nearest even is upper'.

     self assert: 16r1FFFFFFFFFFF1F800 asDouble = 16r1FFFFFFFFFFF20000 asDouble.

     'nearest even is lower'.

     self assert: 16r1FFFFFFFFFFFF0800 asDouble = 16r1FFFFFFFFFFFF0000 asDouble.

    "

!

fromLimitedPrecisionReal:anLPReal

    "return a float with anLPReals value.

     You might loose bits when doing this.

     Slow fallback."

    |fract|

    anLPReal isFinite ifFalse:[

        anLPReal isNaN ifTrue:[^ self NaN].

        anLPReal negative ifTrue:[^ self negativeInfinity].

        ^ self infinity

    ].

    fract := anLPReal asTrueFraction.

    ^ (self fromInteger:fract numerator) / (self fromInteger:fract denominator)

!

fromNumerator:numerator denominator:denominator

    "Create a limited precision real from a Rational.

     This version will answer the nearest flotaing point value,

     according to IEEE 754 round to nearest even default mode"

    |a b q r exponent floatExponent n ha hb hq q1 eMin precision nSubHq|

    a := numerator abs.

    b := denominator abs.

    ha := a highBit.

    hb := b highBit.    

    precision := self precision.

    "If both numerator and denominator are represented exactly in floating point number,

     then fastest thing to do is to use hardwired float division"

    (ha <= precision and:[hb <= precision]) ifTrue: [

        ^ (self coerce:numerator) / (self coerce:denominator)

    ].

    "Try and obtain a mantissa with 1 bit excess precision.

     First guess is rough, we might get one more bit or one less"

    exponent := ha - hb - precision - 1.

    exponent > 0

        ifTrue: [b := b bitShift: exponent]

        ifFalse: [a := a bitShift: exponent negated].

    q := a quo: b.

    r := a - (q * b).

    hq := q highBit.

    "check for gradual underflow, in which case we should use less bits"

    floatExponent := exponent + hq - 1.

    eMin := 2 - (1 bitShift: (self numBitsInExponent - 1)).

    n := floatExponent >= eMin

        ifTrue: [precision + 1]

        ifFalse: [precision + 1 + floatExponent - eMin].

    nSubHq := n - hq.

    nSubHq < 0 ifTrue: [

        exponent := exponent + hq - n.

        r := (q bitAnd: (1 bitShift:nSubHq) - 1) * b + r.

        q := q bitShift:nSubHq

    ] ifFalse:[nSubHq > 0 ifTrue: [

        exponent := exponent + hq - n.

        q1 := (r bitShift:nSubHq) quo: b.

        q := (q bitShift:nSubHq) bitAnd: q1.

        r := r - (q1 * b)

    ]].

    "check if we should round upward.

     The case of exact half (q bitAnd: 1) isZero not & (r isZero)

     will be handled by self fromInteger: conversion"

    ((q bitAnd:1) ~~ 0 and:[r ~~ 0]) ifTrue:[

        q := q + 1

    ].

    (numerator negative xor:denominator negative) ifTrue: [

        q := q negated.

    ].

    ^ (self coerce:q) timesTwoPower:exponent

    "

        Time millisecondsToRun:[

            1000000  timesRepeat:[

                Float fromNumerator:12345678901234567890 denominator:987654321

            ].

        ]

        |fraction|

        fraction := 12345678901234567890//987654321.

        Time millisecondsToRun:[

            1000000  timesRepeat:[

                fraction asFloat

            ].

        ]

    "

!

new:aNumber

    "catch this message - not allowed for floats/doubles"

    self error:'Floats/Doubles cannot be created with new:'

!

readFrom:aStringOrStream onError:exceptionBlock

    "read a float from a string"

    |num|

    num := super readFrom:aStringOrStream onError:exceptionBlock.

    ^ self coerce:num

    "

     Float readFrom:'.1'

     Float readFrom:'0.1'

     Float readFrom:'0'

     ShortFloat readFrom:'.1'

     ShortFloat readFrom:'0.1'

     ShortFloat readFrom:'0'

     LongFloat readFrom:'.1'

     LongFloat readFrom:'0.1'

     LongFloat readFrom:'0'

    "

    "Created: / 07-01-1998 / 16:17:19 / cg"

    "Modified: / 23-11-2010 / 14:35:41 / cg"

! !

!LimitedPrecisionReal class methodsFor:'class initialization'!

initialize

    "initialize ANSI compliant float globals"

    Smalltalk at:#FloatE put:ShortFloat.

    Smalltalk at:#FloatD put:Float.

    Smalltalk at:#FloatQ put:LongFloat.

    "

     self initialize

    "

    "Created: / 7.9.2001 / 14:02:45 / cg"

    "Modified: / 7.9.2001 / 14:02:54 / cg"

! !

!LimitedPrecisionReal class methodsFor:'constants & defaults'!

NaN

    "return the constant NaN (not a Number)."

    ^ self zero uncheckedDivide:self zero

    "

      ShortFloat NaN  

      Float NaN       

      LongFloat NaN   

      LargeFloat NaN   

    "

!

computeEpsilon

    "return the maximum relative spacing"

    ^ self radix asFloat raisedToInteger:(1 - self precision).

    "

      ShortFloat computeEpsilon  

      Float computeEpsilon       

      LongFloat computeEpsilon   

    "

!

e

    "return the closest approximation of the irrational number e"

    self subclassResponsibility

    "Modified: / 7.9.2001 / 14:05:02 / cg"

!

emax

    "return the largest exponent"

    self subclassResponsibility

    "Created: / 7.9.2001 / 14:05:28 / cg"

!

emin

    "return the smallest exponent"

    self subclassResponsibility

    "Created: / 7.9.2001 / 14:05:35 / cg"

!

epsilon

    "return the maximum relative spacing"

    ^ self radix asFloat raisedToInteger:(1 - self precision)

    "

      ShortFloat epsilon  

      Float epsilon       

      LongFloat epsilon   

    "

!

fmax

    "return the largest value allowed"

    self subclassResponsibility

    "Created: / 7.9.2001 / 14:06:56 / cg"

!

fmin

    "return the minimum value allowed"

    self subclassResponsibility

    "Created: / 7.9.2001 / 14:07:06 / cg"

!

infinity

    "return a float representing infinity"

    ^ self unity uncheckedDivide:self zero

    "

      ShortFloat infinity  

      Float infinity       

      LongFloat infinity   

      LargeFloat infinity   

    "

!

negativeInfinity

    "return a float representing negative infinity"

    ^ self unity negated uncheckedDivide:self zero

    "

      ShortFloat negativeInfinity   

      Float negativeInfinity       

      LongFloat negativeInfinity   

      LargeFloat negativeInfinity   

    "

!

pi

    "return the closest approximation of the irrational number pi"

    self subclassResponsibility

    "Created: / 7.9.2001 / 14:07:35 / cg"

! !

!LimitedPrecisionReal class methodsFor:'queries'!

decimalPrecision

    "return the number of valid decimal digits"

    ^ ((self precision + 1) / (10.0 log:self radix)) ceiling   "use ceiling, since the last decimal digit contains some precision"

    "

     ShortFloat decimalPrecision 

     Float decimalPrecision      

     LongFloat decimalPrecision  

    "

!

denormalized

    "Return whether the instances of this class can

     represent values in denormalized format."

    ^ true

!

isAbstract

    "Return if this class is an abstract class.

     True is returned for LimitedPrecisionReal here; false for subclasses."

    ^ self == LimitedPrecisionReal

    "

     1.0 class isAbstract

    "

!

numBitsInExponent

    "return the number of bits in the exponent"

    ^ self subclassResponsibility

!

numBitsInIntegerPart

    ^ 0

!

numBitsInMantissa

    "return the number of bits in the mantissa

     (typically 1 less than the precision due to the hidden bit)"

    ^ self subclassResponsibility

!

precision

    "return the number of valid mantissa bits"

    ^ self numBitsInMantissa + 1 - self numBitsInIntegerPart 

!

radix

    "return the radix (base)"

    self subclassResponsibility

    "Created: / 7.9.2001 / 14:08:20 / cg"

! !

!LimitedPrecisionReal methodsFor:'accessing'!