"{ Encoding: utf8 }"

"

 COPYRIGHT (c) 2003 by eXept Software AG

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

LimitedPrecisionReal subclass:#LargeFloat

	instanceVariableNames:'biasedExponent mantissa precision'

	classVariableNames:'Zero One NaN PositiveInfinity NegativeInfinity Pi_1000 E_1000

		Ln10 Pi_130 DefaultPrecision DefaultPrintPrecision'

	poolDictionaries:''

	category:'Magnitude-Numbers'

!

!LargeFloat class methodsFor:'documentation'!

copyright

"

 COPYRIGHT (c) 2003 by eXept Software AG

              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

"

    Attention:

        Experimental & Unfinished Code.

        The implementation is neither complete nor tuned for performance - still being developed.

    This class provides arbitrary precision floats. 

    I store floating point numbers in base 2 with some arbitrary precision (arbitrary number of bits).

    I do inexact arithmetic like Float.

    But I am very slow due to emulated (Large) Integer arithmetic... (compared to IEEE 754 hardwired)

    Unlike Float, mantissa is not normalized to the form 1.mmmmmm

    It is just stored as an integer.

    The sign is stored in the mantissa.

    biasedExponent is the power of two that multiply the mantissa to form the number. 

    There is no limitation of exponent (overflow or underflow), 

    unless you succeed in exhausting the VM memory...

    Like Float, my arithmetic operations are inexact. 

    They will round to nearest precision LargeFloat.

    If two different precisions are used in arithmetic, 

    the result is expressed in the higher precision.

    If the number is exact (i.e. was constructed from an integer),

    I store 0 as precision.

    Default operating mode is rounding, but might be one of the other possibility 

    (truncate floor ceiling).

    Instance Variables:

            mantissa        <Integer>       the bits of mantissa including sign

            exponent        <Integer>       the times two power to multiply the mantissa (floating binary scale)

            precision       <Magnitude>     number of bits to be stored in mantissa when I am normalized

    [author:]

        Claus Gittinger with code merged in from newer squeak code

        (which itself seems to be based on older ST/X code ;-)

    [see also:]

        Number

        Float LongFloat ShortFloat Fraction FixedPoint 

        SmallInteger LargeInteger

"

!

examples

"

  pi should be (500 digits below):

       3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912

       '%500.498e' printf:{ (1.0 asLargeFloatPrecision:500) pi }

  1000 factorial as a LargeFloat:

     (1 to:1000) inject:1 asLargeFloat into:[:p :m | p * m]          

  1000 factorial as an Integer:

     (1 to:1000) inject:1 into:[:p :m | p * m]                 

  compute 20000.0 factorial:

     Time millisecondsToRun:[ 

        (1 to:20000) inject:1 asLargeFloat into:[:p :m | p * m]

     ] -> 210

  compute 20000 factorial:

     Time millisecondsToRun:[ 

        (1 to:20000) inject:1 into:[:p :m | p * m]

     ] -> 160 

"

! !

!LargeFloat class methodsFor:'instance creation'!

fromFraction:aFraction precision:n

    "Note: form below would not be the closest approximation

    ^ (aFraction numerator asLargeFloatPrecision: n)

            inPlaceDivideBy: (aFraction denominator asLargeFloatPrecision: n)"

    | a b mantissa exponent nBits ha hb hm hasTruncatedBits numerator denominator |

    numerator := aFraction numerator.

    denominator := aFraction denominator.

    a := numerator abs.

    b := denominator abs.

    ha := a highBit.

    hb := b 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 := n + 1.

    (ha < nBits and: [hb < nBits]) 

            ifTrue: 

                    [^(numerator asLargeFloatPrecision: n) 

                            inPlaceDivideBy: (denominator asLargeFloatPrecision: n)].

    "Shift the fraction by a power of two exponent so as to obtain a mantissa with n+1 bits.

    First guess is rough, the mantissa might have n+2 bits."

    exponent := ha - hb - nBits.

    exponent > 0 

            ifTrue: [b := b bitShift: exponent]

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

    mantissa := a quo: b.

    hasTruncatedBits := a > (mantissa * b).

    hm := mantissa highBit.

    "Remove excess bits in the mantissa."

    hm > nBits 

            ifTrue: 

                    [exponent := exponent + hm - nBits.

                    hasTruncatedBits := hasTruncatedBits or: [mantissa anyBitOfMagnitudeFrom: 1 to: hm - nBits].

                    mantissa := mantissa bitShift: nBits - hm].

    "Check if mantissa must be rounded upward.

    The case of tie (mantissa odd & hasTruncatedBits not)

    will be handled by Integer>>asLargeFloatPrecision:."

    (hasTruncatedBits and: [mantissa odd])

            ifTrue: [mantissa := mantissa + 1].

    "build the LargeFloat from mantissa and exponent"

    ^(aFraction positive 

            ifTrue: [mantissa asLargeFloatPrecision: n]

            ifFalse: [mantissa negated asLargeFloatPrecision: n]) 

                    inPlaceTimesTwoPower: exponent

    "

     (1/2) asLargeFloat       

    "

    "Created: / 26-05-2019 / 03:54:31 / Claus Gittinger"

    "Modified (format): / 28-05-2019 / 16:13:06 / Claus Gittinger"

!

fromInteger:anInteger

    ^ (self basicNew 

        mantissa:anInteger 

        exponent:0

        precision:nil) normalize

    "

     LargeFloat fromInteger:123456

     1 asLargeFloat       

     2 asLargeFloat       

     16 asLargeFloat       

     1000 factorial asLargeFloat             

    "

    "Modified: / 27-05-2019 / 16:47:15 / Claus Gittinger"

!

fromInteger:anInteger precision:n

    ^ (self basicNew 

        mantissa:anInteger 

        exponent:0

        precision:n) normalize

    "

     LargeFloat fromInteger:123456 precision:6

     1000 factorial asLargeFloat  

     1 asLargeFloat == 1 asLargeFloat

     1.0 asLargeFloat == 1.0 asLargeFloat

    "

    "Created: / 26-05-2019 / 03:54:44 / Claus Gittinger"

    "Modified: / 27-05-2019 / 16:46:59 / Claus Gittinger"

!

fromLimitedPrecisionReal:aLimitedPrecisionReal

    "assumes that aLimitedPrecisionReal is a regular IEEE floating pnt number,

     with a biased base2 exponent and hidden mantissa bit"

    |shifty numBytes numBitsInMantissa maskMantissa numBitsInExponent maskExponent

     numIntegerBits numBits biasExponent sign expPart fractionPart fraction exp|

    aLimitedPrecisionReal isFinite ifFalse:[

        aLimitedPrecisionReal isNaN ifTrue:[^ self NaN].

        aLimitedPrecisionReal > 0 ifTrue:[^ self infinity].

        ^ self negativeInfinity

    ].

    numBytes := aLimitedPrecisionReal basicSize.

    numBitsInMantissa := aLimitedPrecisionReal numBitsInMantissa. maskMantissa := (1 bitShift:numBitsInMantissa) - 1.

    numBitsInExponent := aLimitedPrecisionReal numBitsInExponent. maskExponent := (1 bitShift:numBitsInExponent) - 1.

    numIntegerBits := aLimitedPrecisionReal numBitsInIntegerPart.

    numBits := numBitsInMantissa + numBitsInExponent.

    biasExponent := maskExponent bitShift:-1.

    "/ same as:

    "/ biasExponent := aLimitedPrecisionReal class eBias.

    shifty := LargeInteger basicNew numberOfDigits:numBytes.

    UninterpretedBytes isBigEndian ifTrue:[

        1 to:numBytes do:[:i | shifty digitAt:(numBytes+1-i) put:(aLimitedPrecisionReal basicAt:i)].

    ] ifFalse:[

        1 to:numBytes do:[:i | shifty digitAt:i put:(aLimitedPrecisionReal basicAt:i)].

    ].

    sign := (shifty bitAt:numBits+1) == 0 ifTrue: [1] ifFalse: [-1].

    expPart := (shifty bitShift:numBitsInMantissa negated) bitAnd: maskExponent.

    fractionPart := shifty bitAnd:maskMantissa.

    ( expPart=0 and: [ fractionPart=0 ] ) ifTrue: [ ^ self zero  ].

    numIntegerBits == 0 ifTrue:[

        " Replace omitted leading 1 in fraction (Notice: quadIEEE format does not do this)"

        fraction := fractionPart bitOr: (maskMantissa + 1).

    ] ifFalse:[

        fraction := fractionPart.

    ].

    "Unbias exponent"

    exp := biasExponent - expPart + (numBitsInMantissa - numIntegerBits).

    sign < 0 ifTrue:[ fraction := fraction negated ].

    ^ (self basicNew 

        mantissa:fraction 

        exponent:(exp negated)

        precision:(aLimitedPrecisionReal precision))  normalize    

    "

     1.0 asLargeFloat

     123456 asLargeFloat /7

     take a look at the precision...

     1.0 asShortFloat asLargeFloat precision      

     1.0 asLongFloat asLargeFloat precision      

     1.0 asQDouble asLargeFloat precision      

     1 asLargeFloat       

     (5/3) asLargeFloat       

     (3/5) asLargeFloat       

     (1/2) asLargeFloat       

     2.0 asLargeFloat       

     20000.0 asLargeFloat   

     2e6 asLargeFloat                                  

     1e300 asLargeFloat             

     2e300 asLargeFloat             

     0.5 asLargeFloat      

     0.25 asLargeFloat     

     (1.0/20000.0) asLargeFloat 

     2e-6 asLargeFloat        

     2e-300 asLargeFloat      

     -1.0 asLargeFloat       

     -0.5 asLargeFloat      

     Float NaN asLargeFloat              

     Float infinity asLargeFloat         

     Float negativeInfinity asLargeFloat 

    "

    "Modified (comment): / 17-07-2017 / 15:09:30 / cg"

    "Modified: / 28-05-2019 / 09:04:35 / Claus Gittinger"

!

mantissa:m exponent:e

    ^ (self basicNew mantissa:m exponent:e) normalize

    "

     LargeFloat mantissa:1 exponent:0 

     LargeFloat mantissa:2 exponent:0 

     LargeFloat mantissa:4 exponent:0   

     LargeFloat mantissa:8 exponent:0 

     LargeFloat mantissa:1 exponent:-1

     LargeFloat mantissa:1 exponent:-2

     LargeFloat mantissa:1 exponent:-3

    "

    "Modified: / 27-05-2019 / 16:48:33 / Claus Gittinger"

!

mantissa:m exponent:e precision:p

    ^ (self basicNew mantissa:m exponent:e precision:p) normalize

    "Modified: / 27-05-2019 / 16:48:39 / Claus Gittinger"

! !

!LargeFloat class methodsFor:'class initialization'!

initialize

    NaN := self mantissa:0 exponent:9999999.

    PositiveInfinity := self mantissa:0 exponent:1.

    NegativeInfinity := self mantissa:0 exponent:-1.

    One := self mantissa:1 exponent:0.

    Zero := self mantissa:0 exponent:0.

    DefaultPrecision := 200.

    "

     LargeFloat initialize

    "

    "Modified: / 10-10-2017 / 15:56:36 / cg"

! !

!LargeFloat class methodsFor:'coercing & converting'!

coerce:aNumber

    "convert the argument aNumber into an instance of the receiver (class) and return it."

    ^ aNumber asLargeFloat

! !

!LargeFloat class methodsFor:'constants'!

NaN

    ^ NaN

    "

     LargeFloat NaN

     (0.0 uncheckedDivide:0.0)

    "

!

infinity

    ^ PositiveInfinity 

    "

     LargeFloat infinity

     (1.0 uncheckedDivide:0.0)

    "

!

ln10

    "return the ln(10) as largeFloat with approx. 500 bits of precision"

    Ln10 isNil ifTrue:[

        Ln10 := (10.0 asLargeFloat ln_withAccuracy:(LongFloat readFrom:'1e-100')) precision:100

    ].

    ^ Ln10

    "

     LargeFloat ln10

    "

    "Created: / 17-07-2017 / 15:15:25 / cg"

    "Modified: / 17-07-2017 / 16:42:27 / cg"

!

negativeInfinity

    ^ NegativeInfinity

    "

     LargeFloat negativeInfinity

     (-1.0 uncheckedDivide:0.0)

    "

!

one

    ^ One

    "

     LargeFloat one

    "

    "Created: / 26-05-2019 / 03:47:18 / Claus Gittinger"

!

pi

    Pi_130 isNil ifTrue:[

       Pi_130 := (self fromString:'3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647')

                    precision:130.

    ].

    ^ Pi_130.

    "

     LargeFloat pi

    "

    "Modified: / 24-07-2019 / 09:53:31 / Claus Gittinger"

!

unity

    ^ One

    "

     LargeFloat unity

    "

!

zero

    ^ Zero

    "

     LargeFloat zero

    "

! !

!LargeFloat class methodsFor:'queries'!

defaultPrecision

    ^ DefaultPrecision

    "Created: / 10-10-2017 / 15:58:03 / cg"

!

defaultPrintPrecision

    ^ DefaultPrintPrecision ? 12

    "Created: / 10-10-2017 / 15:58:03 / cg"

!

eBias

    "Answer the exponent's bias;

     that is the offset of the zero exponent when stored"

    ^ 0

!

emax

    "there is no maximum"

    ^ Infinity positiveInfinity

!

emin

    "there is no minimum"

    ^ Infinity negativeInfinity

!

radix

   "answer the radix of a LargeFloat's exponent"

    ^ 2 

! !

!LargeFloat methodsFor:'accessing'!

biasedExponent

    "anwser the raw exponent;

     this is not what a standard floating point representation expects"

    ^ biasedExponent

    "Created: / 27-05-2019 / 16:38:43 / Claus Gittinger"

!

decimalPrecision

    "Answer how many digits of accuracy this class supports"

    ^ (precision * (2 log: 10)) truncated

!

exponent

    "anwser the floating point like exponent e,

     of self normalized as

         mantissa * (2 raisedTo: e)"

    self isZero ifTrue: [^0].

    ^biasedExponent + self precisionInMantissa - 1

!

mantissa

    ^ mantissa

!

nextToward: aNumber