LimitedPrecisionReal.st
author Claus Gittinger <cg@exept.de>
Mon, 16 Jun 2003 11:17:33 +0200
changeset 7356 fe8fb0a571f2
parent 7141 033f2c26d8e6
child 7378 e71a28f6b712
permissions -rw-r--r--
double dispatching fixed; many refactorings

"
 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 single and double (and maybe more) 
    precision real numbers (i.e. Float and Double).

    Due to historic reasons, ST/X's Floats are what Doubles are in ST-80.
    This may change soon (implementing LPReal is a first step towards this).

    [author:]
        Claus Gittinger

    [see also:]
        Fraction FixedPoint
"
! !

!LimitedPrecisionReal class methodsFor:'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:'instance creation'!

new:aNumber
    "catch this message - not allowed for floats/doubles"

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

!LimitedPrecisionReal class methodsFor:'constants & defaults'!

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 raisedTo:(1 - self precision)

    "
      ShortFloat epsilon  
      Float epsilon       
      LongFloat epsilon   
    "


    "Created: / 7.9.2001 / 14:06:41 / cg"
!

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

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 / (10.0 log:self radix)) floor

    "
     ShortFloat decimalPrecision 
     Float decimalPrecision      
     LongFloat decimalPrecision  
    "
!

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

at:index
    "redefined to prevent access to individual bytes in a real."

    self error:'not allowed for floats/doubles'
!

at:index put:aValue
    "redefined to prevent access to individual bytes in a real"

    self error:'not allowed for floats/doubles'
! !

!LimitedPrecisionReal methodsFor:'arithmetic'!

* aNumber
    "return the product of the receiver and the argument, aNumber"

"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ aNumber productFromDouble:self asDouble

    ^ aNumber productFromFloat:self asFloat

    "Modified: 17.4.1996 / 12:35:36 / cg"
!

+ aNumber
    "return the sum of the receiver and the argument, aNumber"

"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ aNumber sumFromDouble:self asDouble

    ^ aNumber sumFromFloat:self asFloat

    "Modified: 17.4.1996 / 12:35:55 / cg"
!

- aNumber
    "return the difference of the receiver and the argument, aNumber"

"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ aNumber differenceFromDouble:self asDouble

    ^ aNumber differenceFromFloat:self asFloat

    "Modified: 17.4.1996 / 12:36:07 / cg"
!

/ aNumber
    "return the quotient of the receiver and the argument, aNumber"

    ((aNumber == 0) or:[aNumber = 0.0]) ifTrue:[
        ^ ZeroDivide raiseRequestWith:thisContext.
    ].
"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ aNumber quotientFromDouble:self asDouble

    ^ aNumber quotientFromFloat:self asFloat

    "Modified: / 17.4.1996 / 12:36:21 / cg"
    "Modified: / 26.7.1999 / 10:46:11 / stefan"
!

// aNumber
    "return the integer quotient of dividing the receiver by aNumber with
     truncation towards negative infinity."

    ^ (self / aNumber) floor asInteger

    "Modified: 5.11.1996 / 11:45:37 / cg"
!

ceiling
    self subclassResponsibility
!

floor
    self subclassResponsibility
! !

!LimitedPrecisionReal methodsFor:'coercing & converting'!

asFixedPoint
    "return the receiver as fixedPoint number.
     Q: what should the scale be here ?"

    ^ self asFraction asFixedPoint

    "
     0.3 asFixedPoint
     0.5 asFixedPoint
     (1/5) asFloat asFixedPoint 
     (1/3) asFloat asFixedPoint 
     (2/3) asFloat asFixedPoint 
     (1/8) asFloat asFixedPoint
     3.14159 asFixedPoint
     0.0000001 asFraction
     0.0000001 asFixedPoint
    "

    "Modified: / 25.10.1997 / 15:36:54 / cg"
!

asFixedPoint:scale
    "return the receiver as fixedPoint number with the given
     number of post-decimal-digits."

    ^ self asFraction asFixedPoint:scale

    "
     0.3 asFixedPoint:4     
     0.3 asFixedPoint:3     
     0.3 asFixedPoint:2     
     0.3 asFixedPoint:1     
     0.3 asFixedPoint:0

     0.5 asFixedPoint:3     
     (1/5) asFloat asFixedPoint:1  
     (1/8) asFloat asFixedPoint:1  
     1.0 asFixedPoint:2 
     3.14159 asFixedPoint:2       
     3.14159 asFixedPoint:3       
     (3.14159 asFixedPoint:2) asFixedPoint:5  
    "

    "Modified: / 5.8.1998 / 13:29:51 / cg"
!

asFraction
    "Answer a rational number (Integer or Fraction) representing the receiver.
     This conversion uses the continued fraction method to approximate 
     a floating point number."

    |num1 denom1 num2 denom2 int frac newD temp limit|

    limit := (self class unity * 10) raisedTo:self defaultNumberOfDigits.
    num1 := self truncated. 
    denom1 := 1.            "The first of two alternating denominators"
    num2 := 1.              "The second numerator"
    denom2 := 0.            "The second denominator--will update"
    int := num1.            "The integer part of self"
    frac := self fractionPart.             
    [frac = 0] whileFalse:[                
        newD := 1.0 / frac.                   
        int := newD truncated.        
        frac := newD fractionPart.      "save the fractional part for next time"
        temp := num2.                   "old numerator and save it"
        num2 := num1.                         
        num1 := num1 * int + temp.      "Update first numerator"
        temp := denom2.                 "old denominator and save it"
        denom2 := denom1.                    
        denom1 := int * denom1 + temp.  "Update first denominator"
        limit < denom1 ifTrue:[
            "Is ratio past float precision?  If so, pick which of the two ratios to use"
            num2 = 0.0 ifTrue:[
                "Is second denominator 0?"
                ^ Fraction numerator:num1 denominator:denom1
            ].
            ^ Fraction numerator:num2 denominator:denom2
        ]
    ].

    "If fractional part is zero, return the first ratio"
    denom1 = 1 ifTrue:[
        "Am i really an Integer?"
        ^ num1 "Yes, return Integer result"
    ].
    "Otherwise return Fraction result"
    ^ Fraction numerator:num1 denominator:denom1

    "
     1.1 asFraction      
     1.2 asFraction      
     0.3 asFraction   
     0.5 asFraction 
     (1/5) asFloat asFraction
     (1/8) asFloat asFraction  
     (1/13) asFloat asFraction 
     3.14159 asFraction        
     3.14159 asFraction asFloat       
     1.3 asFraction  
     1.0 asFraction  
    "

    "Modified: / 25.10.1997 / 16:41:19 / cg"
!

asInteger
    "return an integer with same value - might truncate"

    |max maxF|

    self isNaN ifTrue:[
        ^ self class
            raise:#domainErrorSignal
            receiver:self
            selector:#asInteger
            arguments:#()
            errorString:'receiver is NaN in #asInteger'

    ].

    self abs < 2e16 ifTrue:[
        "/ NOTICE: this must be redefined in float
        "/ subclasses to handle the smallinteger range;
        "/ i.e. this may only be invoked for reals
        "/ which are NOT within the smallInt range.
        "/ otherwise, endless recursion is the consequence.

        max := SmallInteger maxVal // 2 + 1.
        maxF := max asFloat.

        ^ (self quo:maxF) * max + (self rem:maxF) truncated
    ].
    ^ self asTrueFraction

    "
     12345.0 asInteger     
     1e15 asInteger        
     1e33 asInteger asFloat
     1e303 asInteger asFloat
    "

    "Modified: / 16.11.2001 / 14:15:33 / cg"
!

asRational
    "Answer a Rational number--Integer or Fraction--representing the receiver.
     Same as asFraction fro st-80 compatibility."

    ^ self asFraction

    "
     1.1 asRational      
     1.2 asRational      
     0.3 asRational   
     0.5 asRational 
     (1/5) asFloat asRational
     (1/8) asFloat asRational  
     (1/13) asFloat asRational 
     3.14159 asRational        
     3.14159 asRational asFloat       
     1.3 asRational  
     1.0 asRational  
    "
!

asTrueFraction
    "Answer a fraction or integer that EXACTLY represents self,
     an any-precision IEEE floating point number, consisting of:
        numMantissaBits bits of normalized mantissa (i.e. with hidden leading 1-bit)
        optional numExtraBits between mantissa and exponent (normalized flag for ext-real)
        numExponentBits bits of 2s complement exponent
        1 sign bit.
     Taken from Floats asTrueFraction"

    |shifty sign expPart exp fraction fractionPart result zeroBitsCount
     numBytes numBits numBitsInMantissa maskMantissa numBitsInExponent maskExponent
     biasExponent numIntegerBits|

    self isFinite ifFalse:[
        ^ self class
            raise:#conversionErrorSignal
            receiver:self
            selector:#asTrueFraction
            arguments:#()
            errorString:'Cannot represent non-finite float as a fraction'.
    ].

    "Extract the bits of an IEEE anySize float "
    numBytes := self basicSize.
    numBitsInMantissa := self class numBitsInMantissa. maskMantissa := (1 bitShift:numBitsInMantissa) - 1.
    numBitsInExponent := self class numBitsInExponent. maskExponent := (1 bitShift:numBitsInExponent) - 1.
    numIntegerBits := self class numBitsInIntegerPart.
    numBits := numBitsInMantissa + numBitsInExponent. "not (numBytes * 8 - 1) -- there could be 80 float bits in 96 overall bits"
    biasExponent := maskExponent bitShift:-1.

    shifty := LargeInteger basicNew numberOfDigits:numBytes.
    UninterpretedBytes isBigEndian ifTrue:[
        1 to:numBytes do:[:i | shifty digitAt:(numBytes+1-i) put:(self basicAt:i)].
    ] ifFalse:[
        1 to:numBytes do:[:i | shifty digitAt:i put:(self basicAt:i)].
    ].

    " Extract the sign and the biased exponent "
    sign := (shifty bitAt:numBits+1) == 0 ifTrue: [1] ifFalse: [-1].
    expPart := (shifty bitShift:numBitsInMantissa negated) bitAnd: maskExponent.

    " Extract fractional part; answer 0 if this is a true 0.0 value "
    fractionPart := shifty bitAnd:maskMantissa.
    ( expPart=0 and: [ fractionPart=0 ] ) ifTrue: [ ^ 0  ].

    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: 16r3FF is bias; numBitsInMantissa (52 for floats) is fraction width"
    exp := biasExponent - expPart + (numBitsInMantissa - numIntegerBits).

    " Form the result. When exp>52, the exponent is adjusted by
      the number of trailing zero bits in the fraction to minimize
      the (huge) time otherwise spent in #gcd:. "
    exp negative ifTrue: [
        result := sign * fraction bitShift: exp negated 
    ] ifFalse:[
        zeroBitsCount := fraction lowBit - 1.
        exp := exp - zeroBitsCount.
        exp <= 0 ifTrue: [
            zeroBitsCount := zeroBitsCount + exp.
            result := sign * fraction bitShift:zeroBitsCount negated 
        ] ifFalse: [
            result := Fraction
                    numerator: (sign * fraction bitShift: zeroBitsCount negated)
                    denominator: (1 bitShift:exp) 
        ] 
    ].

    "Low cost validation omitted after extensive testing"
    "(result asFloat = self) ifFalse: [self error: 'asTrueFraction validation failed']."

    ^ result 

    "
     0.3 asFloat asTrueFraction  
     0.3 asShortFloat asTrueFraction  
     0.3 asLongFloat asTrueFraction  

     1.25 asTrueFraction
     0.25 asTrueFraction  
     -0.25 asTrueFraction  
     3e37 asTrueFraction
     Float NaN asTrueFraction
     Float infinity asTrueFraction
    "
!

fractionalPart
    "This has been renamed to #fractionPart for ST80 compatibility.

     extract the after-decimal fraction part.
     the floats value is 
        float truncated + float fractionalPart"

    <resource:#obsolete>

    self obsoleteMethodWarning:'please use #fractionPart'.
    ^ self fractionPart

    "Modified: / 28.10.1998 / 17:10:12 / cg"
    "Created: / 28.10.1998 / 17:10:32 / cg"
!

generality
    "return the generality value - see ArithmeticValue>>retry:coercing:"

    ^ 80
!

i
    "return a complex number, with the receiver as imaginary part, 0 as real part"

    ^ Complex
        real:0
        imaginary:self

    "
     3.0i
     (1+1i)
    "
! !

!LimitedPrecisionReal methodsFor:'comparing'!

< aNumber
    "return true, if the argument is greater"

"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ aNumber lessFromDouble:self asDouble

    ^ aNumber lessFromFloat:self asFloat

    "Modified: 17.4.1996 / 13:34:50 / cg"
!

<= aNumber
    "return true, if the argument is greater or equal"

    ^ self retry:#<= coercing:aNumber
!

> aNumber
    "return true, if the argument is less"

    ^ self retry:#> coercing:aNumber
!

>= aNumber
    "return true, if the argument is less or equal"

    ^ self retry:#>= coercing:aNumber
! !

!LimitedPrecisionReal methodsFor:'copying'!

deepCopy
    "return a deep copy of myself
     - because storing into floats is not recommended/allowed, its ok to return the receiver"

    ^ self
!

deepCopyUsing:aDictionary
    "return a deep copy of myself
     - because storing into floats is not recommended/allowed, its ok to return the receiver"

    ^ self
!

shallowCopy
    "return a shallow copy of the receiver"

    ^ self
!

simpleDeepCopy
    "return a deep copy of the receiver
     - because storing into floats is not recommended/allowed, its ok to return the receiver"

    ^ self
! !

!LimitedPrecisionReal methodsFor:'double dispatching'!

differenceFromFraction:aFraction
    "sent when a fraction does not know how to subtract the receiver, a float"

    |d|

    d := aFraction denominator.
    ^ (aFraction numerator - (self * d)) / d
!

productFromFraction:aFraction
    "sent when a fraction does not know how to multiply the receiver, a float"

    ^ self * aFraction numerator / aFraction denominator
!

quotientFromFraction:aFraction
    "sent when a fraction does not know how to divide by the receiver, a float"

    ^ aFraction numerator / (self * aFraction denominator)
!

sumFromFraction:aFraction
    "sent when a fraction does not know how to add the receiver, a float"

    |d|

    d := aFraction denominator.
    ^ (self * d + aFraction numerator) / d
! !

!LimitedPrecisionReal methodsFor:'encoding'!

encodeOn:anEncoder with:aParameter

    anEncoder encodeFloat:self with:aParameter


! !

!LimitedPrecisionReal methodsFor:'printing & storing'!

printOn:aStream
    "append a printed representation of the receiver to
     the argument, aStream.

     LimitedPrecisonReal and its subclasses use #printString instead of
     #printOn: as basic print mechanism."

    aStream nextPutAll:self printString

    "Modified: / 20.1.1998 / 14:10:46 / stefan"
!

printString
    "return a printed representation of the receiver
     LimitedPrecisonReal and its subclasses use #printString instead of
     #printOn: as basic print mechanism."

    ^ self subclassResponsibility

    "Created: / 17.4.1996 / 12:12:20 / cg"
    "Modified: / 20.1.1998 / 14:10:47 / stefan"
! !

!LimitedPrecisionReal methodsFor:'queries'!

defaultNumberOfDigits
    "Answer how many digits of accuracy this class supports"

    ^ self class decimalPrecision - 1 
!

size
   "redefined since reals are kludgy (ByteArry)"

   ^ 0
! !

!LimitedPrecisionReal methodsFor:'testing'!

isFinite
   ^ self subclassResponsibility

    "Created: / 7.1.1998 / 12:02:06 / stefan"
!

isFloat
    "return true, if the receiver is some kind of floating point number;
     false is returned here.
     Same as #isLimitedPrecisionReal, but a better name ;-)"

    ^ true

    "Created: / 14.11.2001 / 14:57:55 / cg"
!

isInfinite
    "return true, if the receiver is an infinite float (Inf).
     These are not created by ST/X float operations (they raise an exception);
     however, inline C-code could produce them ..."

    ^ (self isFinite or:[self isNaN]) not.

    "
        1.0 isInfinite
        (0.0 uncheckedDivide: 0.0) isInfinite
        (1.0 uncheckedDivide: 0.0) isInfinite
    "

    "Modified: / 7.1.1998 / 12:01:30 / stefan"
!

isLimitedPrecisionReal
    "return true, if the receiver is some kind of limited precision real (i.e. floating point) number;
     true is returned here - the method is redefined from Object."

    ^ true
!

isNaN
   ^ self subclassResponsibility

    "Modified: 12.2.1997 / 16:45:27 / cg"
!

isNegativeInfinity
    ^ self negative and:[self isInfinite]
!

isNegativeZero
    "many systems have two float.Pnt zeros"

    ^ (self printString first == $-)

    "
     0.0 isNegativeZero
     -0.0 isNegativeZero       
    "
!

isPositiveInfinity
    ^ self positive and:[self isInfinite]
!

negative
    "return true if the receiver is less than zero"

"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ self asDouble negative

    ^ self asFloat negative

    "Modified: 17.4.1996 / 13:35:00 / cg"
!

numberOfBits
    "return the size (in bits) of the real;
     typically, this is 64 for Floats and 32 for ShortFloats,
     but who knows ..."

    self subclassResponsibility
!

positive
    "return true if the receiver is greater or equal to zero"

"/ as soon as Float are float & Double are doubles,
"/ use:
"/    ^ self asDouble positive

    ^ self asFloat positive

    "Modified: 17.4.1996 / 13:35:10 / cg"
! !

!LimitedPrecisionReal methodsFor:'truncation & rounding'!

ceilingAsFloat
    ^ self coerce:(self ceiling).

    "
     0.4 asLongFloat ceilingAsFloat
    "
!

floorAsFloat
    ^ self coerce:(self floor).

    "
     0.4 asLongFloat floorAsFloat
    "
!

roundedAsFloat
    self negative ifTrue:[ ^ (self - 0.5) ceilingAsFloat ].
    ^ (self + 0.5) floorAsFloat
! !

!LimitedPrecisionReal class methodsFor:'documentation'!

version
    ^ '$Header: /cvs/stx/stx/libbasic/LimitedPrecisionReal.st,v 1.48 2003-06-16 09:16:58 cg Exp $'
! !

LimitedPrecisionReal initialize!