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

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

    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

    [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

    ^ (numerator asLargeFloatPrecision: n)

            inPlaceDivideBy: (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)  normalize

    "

     LargeFloat fromInteger:123456

     1 asLargeFloat       

     2 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

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

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

    ^ (self basicNew 

        mantissa:(fraction * sign) 

        exponent:(exp negated)

        precision:(aLimitedPrecisionReal precision))  normalize    

    "

     1.0 asLargeFloat

     take a look at the precision...

     1.0 asShortFloat asLargeFloat       

     1.0 asLongFloat asLargeFloat       

     1.0 asQDouble asLargeFloat       

     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:999.

    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

    "return the argument as a LargeFloat"

    ^ 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_1000 isNil ifTrue:[

        Pi_1000 := FixedPoint pi asLargeFloat

    ].

    ^ Pi_1000.

    "

     LargeFloat pi

    "

!

unity

    ^ One

    "

     LargeFloat unity

    "

!

zero

    ^ Zero

    "

     LargeFloat zero

    "

! !

!LargeFloat class methodsFor:'queries'!

defaultPrecision

    ^ DefaultPrecision

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

!

radix

   "answer the radix of a LargeFloats 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

	^ precision * (2 log: 10)

!

exponent

	"anwser the floating point like exponent e,

	of self normalized as

	1.mmmmmm * (2 raisedTo: e)"

	self isZero ifTrue: [^0].

	^biasedExponent + self precisionInMantissa - 1

!

mantissa

    ^ mantissa

!

nextToward: aNumber 

        "answer the nearest floating point number to self with same precision than self,

        toward the direction of aNumber argument.

        If the nearest one falls on the other side of aNumber, than answer a Number"

        | next |

        "if self is greater, decrease self, but never under aNumber"

        self > aNumber 

                ifTrue: 

                        [next := self predecessor.

                        ^next >= aNumber 

                                ifTrue: [next]

                                ifFalse: [aNumber]].

        "if self is smaller, increase self, but never above aNumber"

        self < aNumber 

                ifTrue: [next := self successor.

                        ^next <= aNumber 

                                ifTrue: [next]

                                ifFalse: [aNumber]].

        "if we are equal, return self"

        ^self

    "Modified (comment): / 28-05-2019 / 16:18:01 / Claus Gittinger"

!

precision

    ^ precision ? 200

!

significand

	^self timesTwoPower: self exponent negated

!

significandAsInteger

	self normalize.

	^mantissa abs

! !

!LargeFloat methodsFor:'arithmetic'!

* aNumber

    ^ aNumber productFromLargeFloat:self

!

+ aNumber

    ^ aNumber sumFromLargeFloat:self