"{ 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' }"

"{ NameSpace: Smalltalk }"

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. IEEE floats and doubles).

    Short summary for beginners (find details in wikipedia):

    ========================================================

    Floating point numbers are represented with a mantissa and an exponent, and the number's value is: 

        mantissa * (2 raisedTo: exponent) 

    with (1 > mantissa >= 0) and exponent adjusted as required for the mantissa to be in that range

    (so called ''normalized'')

    therefore,

        13 asFloat mantissa -> 0.8125

        13 asFloat exponent ->  4  

        0.8125 * (2 raisedTo:4) -> 13

    and:    

        104 asFloat mantissa -> 0.8125

        104 asFloat exponent -> 7  

        0.8125 * (2 raisedTo:7) -> 104

    and:    

        0.1 mantissa -> 0.8

        0.1 exponent -> -3  

        0.8 * (2 raisedTo:-3) -> 0.1

    however:    

        (1 / 3.0) mantissa -> 0.666666666666667

        (1 / 3.0) exponent -> -1  

        0.666666666666667 * (2 raisedTo:-3) -> 0.1

    Danger in using Floats:

    =======================

    Beginners seem to forget (or never learn?) that flt. point numbers are always APPROXIMATIONs of some value.

    You may never ever use them when exact results are neeed (i.e. when computing money!!)

    Take a look at the FixedPoint class for that.

    See also 'Float comparison' below.

    The Float/Double confusion in ST/X:

    ===================================

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

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

    whereas in others, there are both Float and Double classes (VisualWorks).

    In order to allow code from both families to be loaded into ST/X without a missing class error, and without

    loosing precision, we decided to use IEEE doubles as the internal representation of Float 

    and make Double an alias to it.

    This should work for either family (except for the unexpected additional precision in some cases).

    If you really only want single precision floating point numbers, use ShortFloat instances.

    But be aware that there is usually no advantage (neither in memory usage, due to memory alignment restrictions,

    nor in speed), as these days, the CPUs are just as fast doing double precision operations.

    (There might be a noticable difference when doing bulk operations, and you should consider using FloatArray for those).

    Hardware supported precisions

    =============================

    The only really portable sizes are IEEE-single and IEEE-double floats (i.e. ShortFloat and Float instances).

    These are supported on all architectures.

    Some do provide an extended precision floating pnt. number,

    however, the downside is that CPU-architects did not agree on a common format and precision: 

    some use 80 bits, others 96 and others even 128.

    See the comments in the LongFloat class for more details.

    We recommend using Float (i.e. IEEE doubles) unless absolutely required,

    and care for machine dependencies in the code otherwise.

    For higher precision needs, you may also try the new QDouble class, which gives you >200bits (60digits) 

    of precision on all machines (at a noticable performance price, though).

    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... x 10+38   |      6-9

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

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

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

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

      ext    |            |                 |                     |                     |                    |

      (SPARC)|            |                 |                     |                     |                    |

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

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

      ext    |            |                 |                     |                     |                    |

      (x86)  |            |                 |                     |                     |                    |

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

        --   | QDouble    |     --          |   256 / 212         |  2.225... x 10-308  | 1.797... x 10+308  |     >=60

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

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

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

    HalfFloats are only supported in fixed array containers. 

    This was added for OpenGL and other graphic libraries 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 as array containers.

    These formats are seldom used for bulk data.

    QDoubles are special soft floats; slower in performance, but providing 4 times the precision of regular doubles.

    To see the differences in precision:

        '%60.58f' printf:{ 1 asShortFloat exp } -> '2.718281828459045*090795598298427648842334747314453125'          (32 bits)

        '%60.58f' printf:{ 1 asFloat exp }      -> '2.718281828459045*090795598298427648842334747314453125'          (64 bits)

        '%60.58f' printf:{ 1 asLongFloat exp }  -> '2.718281828459045235*4281681079939403389289509505033493041992'   (only 80 valid bits on x86)

        '%60.58f' printf:{ 1 asQDouble exp }    -> '2.71828182845904523536028747135266249775724709369995957496698'   (>200 bits)

        correct value is:                           2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746

    Bulk Containers:

    ================

    If you have a vector or matrix (and especially: large ones) of floating point numbers, the well known

    Array is a very inperformant choice. The reason is that it keeps pointers to each of its elements, and each element

    (if it is a float) is itself stored somewhere in the object memory.

    Thus, there is both a space overhead (every float object has an object header, for class and other information), and

    also a performance overhead (extra indirection, cache misses and alignment inefficiencies).

    For this, the bulk numeric containers are provided, which keep the elements unboxed and properly aligned.

    Use them for matrices and large numeric vectors. They also provide some optimized bulk operation methods,

    such as adding, multiplying etc.

    Take a look at FloatArray, DoubleArray, HalfFloatArray etc.

    Comparing Floats:

    =================

    Due to rounding errors (usually on the last bit(s)), you shalt not compare two floating point numbers

    using the #= operator. For example, the value 0.1 cannot be represented as a sum of powers-of-two fractions,

    and will therefore always be an approximation with a half bit error in the last bit of the mantissa.

    Usually, the print functions take this into consideration and return a (faked) '0.1'.

    However, this half bit error may accumulate, for example, when multiplying that by 0.1 then by 100, 

    the error may get large enough to be no longer pushed under the rug by the print function, 

    and you will get '0.9999999999999' from it.

    Also, comparing against a proper 1.0 (which is representable as an exact power of 2), 

    you will get a false result.

    i.e. (0.1 * 0.1 * 100 ~= 1.0) and (0.1 * 0.1 * 100 - 1.0) ~= 0.0

    This often confuses non-computer scientists (and occasionally even some of those).

    For this, you should always provide an epsilon value, when comparing two non-integer numbers. 

    The epsilon value is the distance you accept two number to be apart to be still considered equal. 

    Effectively the epsilon says are those nearer than this epsilon?.

    Now we could say is the delta between two numbers smaller than 0.00001,

    and get a reasonable answer for big numbers. But what if we compare two tiny numbers?

    Then a reasonable epsilon must also be much smaller!!

    Actually, the epsilon should always be computed dynamically depending on the two values compared.

    That is what the #isAlmostEqualTo:nEpsilon: method does for you. It does not take an absolute epsilon,

    but instead the number of distinct floating point numbers that the two compared floats may be apart.

    That is: the number of actually representable numbers between those two. 

    Effectively, that is the difference between the two mantissas, 

    when the numbers are scaled to the same exponent, taking the number of mantissa bits into account.

    [author:]

        Claus Gittinger

    [see also:]

        Fraction FixedPoint

"

! !

!LimitedPrecisionReal class methodsFor:'instance creation'!

fromInteger:anInteger

    "return a float with anInteger's value.

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

     with a large integer 

     (i.e. when numDigits is above the flt. pnt number's precision)

     (see Float decimalPrecision, LongFloat decimalPrecision."

    |newFloat sign absVal completeAbsVal delta mask|

    newFloat := self zero.

    sign := anInteger sign.

    (sign ~~ 0) ifTrue:[

        "be careful to 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 rightShift: delta.

            "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.

    "

    "Modified: / 25-08-2017 / 12:34:44 / cg"

    "Modified (format): / 08-06-2019 / 02:17:43 / Claus Gittinger"

!

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 rightShift: exponent].

    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

            ].

        ]

    "

    "Modified: / 25-08-2017 / 12:34:55 / cg"

!

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 value].

    ^ 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'

     LimitedPrecisionReal readFrom:'bla' onError:nil

     Float readFrom:'bla' onError:nil

     ShortFloat readFrom:'bla' onError:nil

    "

    "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.

    "/ not defined in ANSI, but consistent in its naming

    Smalltalk at:#FloatQD put:QDouble.

    "

     self initialize

    "

    "Created: / 07-09-2001 / 14:02:45 / cg"

    "Modified (comment): / 21-06-2017 / 13:59:08 / cg"

! !

!LimitedPrecisionReal class methodsFor:'constants & defaults'!

NaN

    "return the constant NaN (not a Number) in my representation.

     Here, based on the assumption that division of zero by zero generates a NaN

     (which is defined as sch in the IEEE standard).

     If a subclass does not, it needs to generate a NaN differently"

    |zero|

    zero := self zero.

    ^ zero uncheckedDivide:zero.

    "

      ShortFloat NaN  

      Float NaN       

      LongFloat NaN   

      LargeFloat NaN   

    "

    "Modified (comment): / 16-06-2017 / 11:03:28 / cg"

    "Modified (comment): / 06-06-2019 / 16:52:43 / Claus Gittinger"

!

computeEpsilon

    "compute the maximum relative spacing of instances of mySelf

     (i.e. the value-delta of the least significant bit).

     See https://en.wikipedia.org/wiki/Machine_epsilon"

    ^ (self coerce:self radix) raisedToInteger:(self precision - 1) negated.

    "

      Float radix 

      Float precision

      ShortFloat computeEpsilon -> 1.192093e-07 

      Float computeEpsilon      -> 2.22044604925031E-16