IEEEFloat.st
author Claus Gittinger <cg@exept.de>
Sat, 02 May 2020 21:40:13 +0200
changeset 5476 7355a4b11cb6
parent 5385 a54bf0f05ff5
permissions -rw-r--r--
#FEATURE by cg class: Socket class added: #newTCPclientToHost:port:domain:domainOrder:withTimeout: changed: #newTCPclientToHost:port:domain:withTimeout:

"{ Encoding: utf8 }"

"
 COPYRIGHT (c) 2018 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:libbasic2' }"

"{ NameSpace: Smalltalk }"

LimitedPrecisionReal variableByteSubclass:#IEEEFloat
	instanceVariableNames:'exponentSize'
	classVariableNames:''
	poolDictionaries:''
	category:'Magnitude-Numbers'
!

!IEEEFloat class methodsFor:'documentation'!

copyright
"
 COPYRIGHT (c) 2018 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
"
    Unfinished, ongoing work

    soft float emulation for arbitrary IEEE float formats.
    This is very very slow and should only be used when importing 
    funny sized floating point numbers (such as float24 or float8) from
    external sources, or to simulate computations on otherwise unsupported floating pnt numbers.

    Create one by passing the overall number of bits and the number of exponent-bits:
        IEEEFloat size:64 exponentSize:(Float numBitsInExponent)
    or:
        1.0 asIEEEFloat
    will give you a (slow) variant of a regular float.

    And:
        IEEEFloat size:256 exponentSize:19
    will give you a 256 bit ieee float.

"
! !

!IEEEFloat class methodsFor:'instance creation'!

fromFloat:aFloat
    |nB f mantissaInclInt biasedExponent intBitMask matissaMask|

    nB := aFloat byteSize.

    "/ if there are integer bits in the characteristic (extended x86 floats),
    "/ we need one more byte and there is no hidden integer bit;    
    aFloat numBitsInIntegerPart > 0 ifTrue:[
        mantissaInclInt := aFloat mantissaBits.
        biasedExponent := aFloat exponent + aFloat eBias.
        matissaMask := (1 bitShift:aFloat numBitsInMantissa)-1. 
        intBitMask := (1 bitShift:aFloat numBitsInMantissa-1). 
        (mantissaInclInt bitTest:intBitMask) ifTrue:[
            "/ turn it off to make it hidden
            mantissaInclInt := mantissaInclInt bitClear:intBitMask.
            mantissaInclInt := mantissaInclInt bitShift:1.
            biasedExponent := biasedExponent - 1.
        ] ifFalse:[
            ((mantissaInclInt = 0) and:[biasedExponent == aFloat eBias]) ifTrue:[
                biasedExponent := 0.
            ] ifFalse:[
                self halt.
            ].
        ].
        ^ self size:nB*8 exponentSize:aFloat numBitsInExponent 
            fraction:mantissaInclInt exponent:biasedExponent signBit:(aFloat < 0 ifTrue:1 ifFalse:0)
    ].

    f := self basicNew:nB.
    1 to:nB do:[:i |
        f basicAt:i put:(aFloat byteAt:i)
    ].
    f exponentSize:(aFloat numBitsInExponent).
    ^ f

    "
     self fromFloat:1.0
     self fromFloat:2.0
     self fromFloat:1.0 asShortFloat
     self fromFloat:2.0 asShortFloat
     self fromFloat:1.0 asLongFloat
     self fromFloat:2.0 asLongFloat
     self fromFloat:3.0 asLongFloat
     self fromFloat:101.0 asLongFloat
     self fromFloat:0.0 asLongFloat
     self fromFloat:0.5 asLongFloat
    "
!

fromInteger:anInteger
    "convert to one of 64, 128 or 256 bit IEEEFloat"

    |nM nE h l fraction e bias absInt nBytes|

    absInt := anInteger abs.
    h := absInt highBit.
    l := absInt lowBit.
    nM := Float numBitsInMantissa.
    nE := Float numBitsInExponent.
    bias := Float eBias.
    nBytes := 64 / 8.
    (h - l) > nM ifTrue:[
        nM := QuadFloat numBitsInMantissa.
        nE := QuadFloat numBitsInExponent.
        bias := QuadFloat eBias.
        nBytes := 128 / 8.
        (h - l) > nM ifTrue:[
            nM := OctaFloat numBitsInMantissa.
            nE := OctaFloat numBitsInExponent.
            bias := OctaFloat eBias.
            nBytes := 256 / 8.
        ]
    ].

    "/ cut off precision by shifting right (if h > nM) 
    "/ or add zeros by shifting to the left 
    h == 0 ifTrue:[
        fraction := e := 0.
    ] ifFalse:[
        fraction := absInt bitShift:(nM - h + 1). 
        fraction := fraction bitAnd:((1 bitShift:nM)-1).
        e := bias + h - 1.
    ].
    ^ (self basicNew:nBytes) exponentSize:nE;
        setFraction:fraction exponent:e 
        signBit:(anInteger negative ifTrue:[1] ifFalse:[0]);
        yourself

    "
     self fromInteger:1
     self fromInteger:-1
     self fromInteger:2
     self fromInteger:1024 * 1024 * 1024 * 1024 * 1024 * 1024
     self fromInteger:1e20 asInteger
     self fromInteger:1e100 asInteger
     self fromInteger:2r1010101010101010101010101010101
     self fromInteger:2r1010101010101010101010101010101010101010101010101010101010101010
     self fromInteger:1e-100 asInteger
     1 asIEEEFloat
    "
!

size:numBits exponentSize:exponentSize
    ^ (self basicNew:(numBits // 8)) exponentSize:exponentSize

    "
     self size:256 exponentSize:19
     self size:16 exponentSize:5
     self size:(1.0 basicSize * 8) exponentSize:(1.0 numBitsInExponent)
    "
!

size:numBits exponentSize:exponentSize fraction:normalizedFraction exponent:biasedExponent signBit:signBit
    ^ ((self basicNew:(numBits // 8)) exponentSize:exponentSize) 
            setFraction:normalizedFraction exponent:biasedExponent signBit:signBit

    "
     self size:256 exponentSize:19 fromInteger:1
     self size:256 exponentSize:19 fromInteger:2
     self size:64 exponentSize:11 fraction:0 exponent:(0 + Float eBias) signBit:0
     self size:64 exponentSize:11 fraction:0 exponent:0 signBit:0
    "
!

size:numBits exponentSize:exponentSize fromFloat:aFloat
    ^ ((self basicNew:(numBits // 8)) exponentSize:exponentSize) setValueFromFloat:aFloat

    "
     self size:256 exponentSize:19 fromFloat:1.0
     self size:256 exponentSize:19 fromFloat:2.0
    "
!

size:numBits exponentSize:exponentSize fromInteger:anInteger
    ^ ((self basicNew:(numBits // 8)) exponentSize:exponentSize) setValueFromInteger:anInteger

    "
     self size:256 exponentSize:19 fromInteger:1
     self size:256 exponentSize:19 fromInteger:2
    "
! !

!IEEEFloat class methodsFor:'coercing & converting'!

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

    ^ aNumber asIEEEFloat
! !

!IEEEFloat class methodsFor:'constants'!

pi
    ^ Float pi
!

unity
    ^ 1.0
!

zero
    ^ 0.0
! !

!IEEEFloat methodsFor:'accessing'!

exponentBits
    "return the bits of my exponent.
     These might be biased.
     My bytes are stored LSB first; thus we have to fetch the high exponent bytes."

    |nBytes nExpBytes e|

    nBytes := self basicSize.
    nExpBytes := (exponentSize + 1 + 7) // 8.
    "/ seee eeee eeee mmmm mmmm ....
    "/ BAD: assumes LSB
    e := (self basicAt:nBytes) bitAnd:16r7F.
    2 to:nExpBytes do:[:eI |
         e := (e bitShift:8) bitOr:(self basicAt:nBytes + 1 - eI)
    ].
    "/ 0eee eeee eeee mmmm
    e := e rightShift:(nExpBytes * 8) - 1 "sign" - exponentSize.
    "/ 0000 0eee eeee eeee 
    ^ e

    "
     1.0 asShortFloat asIEEEFloat
     (1 to:4) conform:[:i | (1.0 asShortFloat basicAt:i) = (1.0 asShortFloat asIEEEFloat basicAt:i)]
     1.0 asShortFloat exponentBits 
     1.0 asShortFloat asIEEEFloat exponentBits
     -1.0 asShortFloat exponentBits  
     -1.0 asShortFloat asIEEEFloat exponentBits   

     (self size:32 exponentSize:8) inspect.
     (self size:64 exponentSize:11) inspect.
     (self size:128 exponentSize:16) inspect.
     (self size:256 exponentSize:16) inspect.
    "

    "construct a soft-shortFloat

     |f sf|
     sf := 1.0 asShortFloat.
     f := self size:(sf basicSize * 8) exponentSize:(sf numBitsInExponent).
     1 to:sf basicSize do:[:i | f basicAt:i put:(sf basicAt:i)].
     f inspect
    "

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

exponentSize
    ^ exponentSize
!

exponentSize:something
    exponentSize := something.
!

mantissaBits
    "return the bits of my mantissa (incl. the hidden bit).
     My bytes are stored LSB first."

    |bits hiddenBit mask c|

    bits := self digitBytes asIntegerMSB:false.
    bits == 0 ifTrue:[^ 0].

    hiddenBit := (1 bitShift:self numBitsInMantissa).
    mask := hiddenBit-1.
    c := bits bitAnd:mask.
    "/ TODO: care for subnormals (exponent = 0)!!
    ^ c bitOr:hiddenBit.

    " 
     (self size:32 exponentSize:8) inspect.
    "

   "construct a soft-shortFloat

    |f sf|
    sf := 1.0 asShortFloat.
    f := self size:(sf basicSize * 8) exponentSize:(sf numBitsInExponent).
    1 to:sf basicSize do:[:i | f basicAt:i put:(sf basicAt:i)].
    f inspect
   "

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

setFraction:normalizedFraction exponent:biasedExponent signBit:signBit
    |bits|

    bits := (normalizedFraction 
            bitOr:(biasedExponent bitShift:self numBitsInMantissa))
            bitOr:(signBit bitShift:(self basicSize * 8) - 1). 
    1 to:self basicSize do:[:i | self basicAt:i put:(bits digitAt:i)].

    "
     (self size:64 exponentSize:11) setFraction:0 exponent:1024 signBit:0
    "
!

setValueFromInteger:intValue
    "/ how many bits are there, in this int
    |absValue myNumBits numBitsInNumber shift|

    absValue := intValue abs. 
    numBitsInNumber := absValue highBit.
    myNumBits := (self basicSize * 8) - 1 "sign" - exponentSize.
    shift := myNumBits - numBitsInNumber. 
    numBitsInNumber > myNumBits ifTrue:[
        self halt.
    ] ifFalse:[
        "/ number:
        "/    1xxxxxxx...xxxxx
        "/ myRep:
        "/    seee...eeexxxxxxxxxx
        absValue digitLength to:1 by:-1 do:[:byteIndex |
            
        ].
    ].

    "/ cut off some bits

    "Float numBitsInExponent
     self size:16 exponentSize:4 fromInteger:1
     self size:32 exponentSize:11 fromInteger:1

     self size:256 exponentSize:19 fromInteger:1
     self size:256 exponentSize:19 fromInteger:2
    "
! !

!IEEEFloat methodsFor:'arithmetic'!

* aNumber
    "/ thisContext isRecursive ifTrue:[self halt].
    ^ aNumber productFromIEEEFloat:self
!

+ aNumber
    "/ thisContext isRecursive ifTrue:[self halt].
    ^ aNumber sumFromIEEEFloat:self

    "
     (-2.0 asIEEEFloat)
     (2.0 asIEEEFloat)

     self assert:((4.0 asIEEEFloat) + (2.0 asIEEEFloat)) = (6.0 asIEEEFloat)   
     self assert:((4.0 asIEEEFloat) + (-2.0 asIEEEFloat)) = (2.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) + (-2.0 asIEEEFloat)) = (-6.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) + (2.0 asIEEEFloat)) = (-2.0 asIEEEFloat)  

     self assert:((2.0 asIEEEFloat) + (4.0 asIEEEFloat)) = (6.0 asIEEEFloat)   
     self assert:((2.0 asIEEEFloat) + (-4.0 asIEEEFloat)) = (-2.0 asIEEEFloat)  
     self assert:((-2.0 asIEEEFloat) + (-4.0 asIEEEFloat)) = (-6.0 asIEEEFloat)  
     self assert:((-2.0 asIEEEFloat) + (4.0 asIEEEFloat)) = (2.0 asIEEEFloat)  

     self assert:((IEEEFloat fromFloat:0.1) + (IEEEFloat fromFloat:-0.1)) = (IEEEFloat fromFloat:0.0)  
     self assert:((IEEEFloat fromFloat:-0.1) + (IEEEFloat fromFloat:0.1)) = (IEEEFloat fromFloat:0.0)  
    "
!

- aNumber
    "/ thisContext isRecursive ifTrue:[self halt].
    ^ aNumber differenceFromIEEEFloat:self

    "
     self assert:((4.0 asIEEEFloat) - (2.0 asIEEEFloat)) = (2.0 asIEEEFloat)   
     self assert:((4.0 asIEEEFloat) - (-2.0 asIEEEFloat)) = (6.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) - (-2.0 asIEEEFloat)) = (-2.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) - (2.0 asIEEEFloat)) = (-6.0 asIEEEFloat)  

     self assert:((2.0 asIEEEFloat) - (4.0 asIEEEFloat)) = (-2.0 asIEEEFloat)   
     self assert:((2.0 asIEEEFloat) - (-4.0 asIEEEFloat)) = (6.0 asIEEEFloat)  
     self assert:((-2.0 asIEEEFloat) - (-4.0 asIEEEFloat)) = (2.0 asIEEEFloat)  
     self assert:((-2.0 asIEEEFloat) - (4.0 asIEEEFloat)) = (-6.0 asIEEEFloat)  

     self assert:((IEEEFloat fromFloat:0.1) - (IEEEFloat fromFloat:0.1)) = (IEEEFloat fromFloat:0.0)  
     self assert:((IEEEFloat fromFloat:0.1) - (IEEEFloat fromFloat:-0.1)) = (IEEEFloat fromFloat:0.2)  
     self assert:((IEEEFloat fromFloat:-0.1) - (IEEEFloat fromFloat:-0.1)) = (IEEEFloat fromFloat:0.0)  
    "
!

/ aNumber
    "/ thisContext isRecursive ifTrue:[self halt].
    ^ aNumber quotientFromIEEEFloat:self
!

negated
    |f highByte|

    f := self copy.
    highByte := self basicAt:(self basicSize).
    f basicAt:(self basicSize) put:(highByte bitXor:16r80).
    ^ f.

    "
     self assert:(IEEEFloat fromFloat:2.0) negated = (IEEEFloat fromFloat:-2.0) 
     self assert:(IEEEFloat fromFloat:2.0) negated negated = (IEEEFloat fromFloat:2.0) 
     self assert:(IEEEFloat fromFloat:-2.0) negated = (IEEEFloat fromFloat:2.0) 
    "
! !

!IEEEFloat methodsFor:'coercing & converting'!

asIEEEFloat
    ^ self
!

coerce:aNumber
    "convert the argument aNumber into an instance of the receiver's class and return it.
     Redefined to preserve the size and exponentSize"

    |fraction e h nM|

    aNumber isInteger ifTrue:[
        nM := self numBitsInMantissa.
        "/ cut off precision by shifting right (if h > nM) 
        "/ or add zeros by shifting to the left 
        h := aNumber highBit.
        h == 0 ifTrue:[
            fraction := e := 0.
        ] ifFalse:[
            fraction := aNumber bitShift:(nM - h + 1). 
            fraction := fraction bitAnd:((1 bitShift:nM)-1).
            e := self eBias + h - 1.
        ].
        ^ (self class basicNew:self basicSize) exponentSize:exponentSize;
            setFraction:fraction exponent:e 
            signBit:(aNumber negative ifTrue:[1] ifFalse:[0]);
            yourself
    ].
    ^ self class coerce:aNumber

    "Modified: / 15-06-2017 / 10:27:03 / cg"
!

generality
    ^ 97   "/ between OctaFloat and LargeFloat

    "
     1 asShortFloat generality 70
     1 asFloat generality      80
     1 asLongFloat generality  90
     1 asQuadFloat generality
     1 asQDouble generality    95
     1 asLargeFloat generality 100
    "
! !

!IEEEFloat methodsFor:'comparing'!

< aNumber
    "/ thisContext isRecursive ifTrue:[self halt].
    ^ aNumber lessFromIEEEFloat:self
!

= aNumber
    ^ aNumber equalFromIEEEFloat:self
! !

!IEEEFloat methodsFor:'double dispatching'!

differenceFromIEEEFloat:anIEEEFloat
    "/ anIEEEFloat - self
    self isZero ifTrue:[^ anIEEEFloat]. "/ otherwise, we might get a negative zero
    ^ self negated sumFromIEEEFloat:anIEEEFloat

    "
     0.0 + (0.0 negated)         
     (0.0 negated) + 0.0 
     0.0 asIEEEFloat + (0.0 negated)         
     (0.0 asIEEEFloat negated) + 0.0 
     0.0 asIEEEFloat + (0.0 asIEEEFloat negated)         
     (0.0 asIEEEFloat negated) + 0.0 asIEEEFloat
    "
!

equalFromIEEEFloat:anIEEEFloat
    |nBytes m1 m2 nM1 nM2 e1 e2|

    nBytes := self basicSize.
    anIEEEFloat basicSize == nBytes ifTrue:[
        anIEEEFloat exponentSize == exponentSize ifTrue:[
            1 to:nBytes-1 do:[:i |
                (self basicAt:i) = (anIEEEFloat basicAt:i) ifFalse:[^ false].
            ].
            "/ care for negative zero
            (self basicAt:nBytes) = (anIEEEFloat basicAt:nBytes) ifFalse:[
                ((self basicAt:nBytes) bitAnd:16r7F) == 0 ifTrue:[
                    ((anIEEEFloat basicAt:nBytes) bitAnd:16r7F) == 0 ifTrue:[
                        ^ true.
                    ].
                ].
                ^ false
            ].
            ^ true.
        ].
    ].

    "/ more complicated compare

    self negative = anIEEEFloat negative ifFalse:[
        ^ self isZero and:[anIEEEFloat isZero]
    ].

    "/ bring to same exponent, add m1 + m2, normalize

    m1 := anIEEEFloat mantissaBits. 
    m2 := self mantissaBits.        

    nM1 := anIEEEFloat numBitsInMantissa.
    nM2 := self numBitsInMantissa.
    nM1 ~= nM2 ifTrue:[
        nM1 > nM2 ifTrue:[
            m2 := m2 bitShift:(nM1-nM2).
        ] ifFalse:[
            m1 := m1 bitShift:(nM2-nM1).
        ].
    ].
    m1 = m2 ifFalse:[^ false].

    e1 := anIEEEFloat exponentBits - anIEEEFloat eBias.
    e2 := self exponentBits - self eBias.
    e1 = e2 ifFalse:[^ false].

    ^ true.
!

lessFromIEEEFloat:anIEEEFloat
    |m1 m2 nM1 nM2 e1 e2 meNegative mySize myByte otherByte isLess|

    meNegative := self negative.
    meNegative = anIEEEFloat negative ifFalse:[
        ^ anIEEEFloat negative
    ].

    "/ same sized floats can be compared easily
    "/ (but care for the negative 0)
    mySize := self basicSize.
    mySize == anIEEEFloat basicSize ifTrue:[
        anIEEEFloat exponentSize == exponentSize ifTrue:[
            "/ ignore the sign bit
            myByte := (self basicAt:mySize) bitAnd:16r7F.
            otherByte := (anIEEEFloat basicAt:mySize) bitAnd:16r7F.    
            otherByte < myByte ifTrue:[
                ^ meNegative not
            ].
            otherByte > myByte ifTrue:[
                ^ meNegative
            ].
            mySize-1 to:1 by:-1 do:[:i |
                myByte := self basicAt:i.
                otherByte := anIEEEFloat basicAt:i.
                otherByte < myByte ifTrue:[
                    ^ meNegative not
                ].
                otherByte > myByte ifTrue:[
                    ^ meNegative
                ].
            ].
            "/ same
            ^ false.
        ].
    ].

    "/ more complicated compare

    "/ bring to same exponent, add m1 + m2, normalize

    m1 := anIEEEFloat mantissaBits. 
    m2 := self mantissaBits.        

    nM1 := anIEEEFloat numBitsInMantissa.
    nM2 := self numBitsInMantissa.
    nM1 ~= nM2 ifTrue:[
        nM1 > nM2 ifTrue:[
            m2 := m2 bitShift:(nM1-nM2).
        ] ifFalse:[
            m1 := m1 bitShift:(nM2-nM1).
        ].
    ].

    e1 := anIEEEFloat exponent.
    e2 := self exponent.
    isLess := e1 = e2 ifTrue:[m1 < m2] ifFalse:[e1 < e2].
    meNegative ifTrue:[^ isLess not].
    ^ isLess.

    "
     self assert:(1.0 asIEEEFloat < 2.0 asIEEEFloat).
     self assert:(1.0 asIEEEFloat < 1.0 asIEEEFloat) not.
     self assert:(1.0 asIEEEFloat < 0.0 asIEEEFloat) not.
     self assert:(1.0 asIEEEFloat < -0.0 asIEEEFloat) not.

     self assert:(-1.0 asIEEEFloat < 2.0 asIEEEFloat).
     self assert:(-1.0 asIEEEFloat < 1.0 asIEEEFloat).
     self assert:(-1.0 asIEEEFloat < 0.0 asIEEEFloat).
     self assert:(-1.0 asIEEEFloat < -0.0 asIEEEFloat).

     self assert:(1.0 asIEEEFloat < -2.0 asIEEEFloat) not.
     self assert:(1.0 asIEEEFloat < -1.0 asIEEEFloat) not.
     self assert:(1.0 asIEEEFloat < 0.0 asIEEEFloat) not.
     self assert:(1.0 asIEEEFloat < -0.0 asIEEEFloat) not.

     self assert:(-1.0 asIEEEFloat < -2.0 asIEEEFloat) not.
     self assert:(-1.0 asIEEEFloat < -1.0 asIEEEFloat) not.
     self assert:(-1.0 asIEEEFloat < 0.0 asIEEEFloat).
     self assert:(-1.0 asIEEEFloat < -0.0 asIEEEFloat).
    "
!

productFromIEEEFloat:anIEEEFloat
    "/ anIEEEFloat * self

    |m1 m2 e1 e2 m e signBit hi shift nM1 nM2 nM|

    "/ m1 * m2 / e1 + e2

    m1 := anIEEEFloat mantissaBits. 
    m2 := self mantissaBits.        
    (m1 = 0 or:[m2 = 0]) ifTrue:[
        "/ care for the sign of the zero
        ^ self class zero
    ].

    nM1 := anIEEEFloat numBitsInMantissa.   
    nM2 := nM := self numBitsInMantissa.
    nM1 ~= nM2 ifTrue:[
        nM1 > nM2 ifTrue:[
            m2 := m2 bitShift:(nM1-nM2).
            nM := nM1
        ] ifFalse:[
            m1 := m1 bitShift:(nM2-nM1).
            nM := nM2
        ].
    ].

    e1 := anIEEEFloat exponentBits.
    e2 := self exponentBits.
    e := (e1 - anIEEEFloat eBias) + (e2 - self eBias).

    "/ 1xxxx....xxxx * 1xxxx....xxxx
    m := m1 * m2.
    "/ now have nB bits too many
    e := e - nM.

    signBit := (self negative == anIEEEFloat negative) ifTrue:[0] ifFalse:[1].
    m == 0 ifTrue:[self halt].

    hi := m highBit.    
    shift := (nM + 1) - hi. 
    m := m bitShift:shift.
    e := e - shift.
    m := m bitAnd:(1 bitShift:nM)-1. 

    e := e + self eBias.
    ^ self class 
        size:(self basicSize*8) exponentSize:exponentSize 
        fraction:m exponent:e signBit:signBit.

    "
     self assert:((IEEEFloat fromFloat:1.0) * (IEEEFloat fromFloat:1.0)) = (IEEEFloat fromFloat:1.0) 
     self assert:((IEEEFloat fromFloat:1.0) * (IEEEFloat fromFloat:1.0)) exponent = 1.0 exponent
     self assert:((IEEEFloat fromFloat:1.0) * (IEEEFloat fromFloat:1.0)) mantissaBits = 1.0 mantissaBits

     self assert:((4.0 asIEEEFloat) * (2.0 asIEEEFloat)) = (8.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) * (2.0 asIEEEFloat)) = (-8.0 asIEEEFloat)  
     self assert:((4.0 asIEEEFloat) * (-2.0 asIEEEFloat)) = (-8.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) * (-2.0 asIEEEFloat)) = (8.0 asIEEEFloat)  

     self assert:((IEEEFloat fromFloat:100.0) * (IEEEFloat fromFloat:100.0)) = (IEEEFloat fromFloat:10000.0) 
     self assert:((IEEEFloat fromFloat:1e-20) * (IEEEFloat fromFloat:1e20)) = (IEEEFloat fromFloat:(1e-20 * 1e20)) 
    "
!

quotientFromIEEEFloat:anIEEEFloat
    "/ anIEEEFloat / self

    |m1 m2 e1 e2 m e signBit nM1 nM2 nM hi shift|

    "/ m1 / m2 / e1 - e2

    m1 := anIEEEFloat mantissaBits.
    m1 == 0 ifTrue:[^ anIEEEFloat ].
    m2 := self mantissaBits.        

    nM1 := anIEEEFloat numBitsInMantissa.
    nM2 := nM := self numBitsInMantissa.
    nM1 ~= nM2 ifTrue:[
        nM1 > nM2 ifTrue:[
            m2 := m2 bitShift:(nM1-nM2).
            nM := nM1
        ] ifFalse:[
            m1 := m1 bitShift:(nM2-nM1).
            nM := nM2
        ].
    ].

    e1 := anIEEEFloat exponentBits.
    e2 := self exponentBits.
    e := (e1 - anIEEEFloat eBias) - (e2 - self eBias).

    "/ 1xxxx....xxxx / 1xxxx....xxxx
    m := (m1 bitShift:nM) // m2.
    "/ rest := m1 \\ m2.

    signBit := (self negative == anIEEEFloat negative) ifTrue:[0] ifFalse:[1].
    m == 0 ifTrue:[self halt].

    hi := m highBit.    
    shift := (nM + 1) - hi. 
    m := m bitShift:shift.
    e := e - shift.
    m := m bitAnd:(1 bitShift:nM)-1. 

    e := e + self eBias.
    ^ self class 
        size:(self basicSize*8) exponentSize:exponentSize 
        fraction:m exponent:e signBit:signBit.

    "
     self assert:((IEEEFloat fromFloat:1.0) / (IEEEFloat fromFloat:1.0)) = (IEEEFloat fromFloat:1.0) 
     self assert:((IEEEFloat fromFloat:1.0) / (IEEEFloat fromFloat:1.0)) exponent = 1.0 exponent
     self assert:((IEEEFloat fromFloat:1.0) / (IEEEFloat fromFloat:1.0)) mantissaBits = 1.0 mantissaBits

     self assert:((4.0 asIEEEFloat) / (2.0 asIEEEFloat)) = (2.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) / (2.0 asIEEEFloat)) = (-2.0 asIEEEFloat)  
     self assert:((4.0 asIEEEFloat) / (-2.0 asIEEEFloat)) = (-2.0 asIEEEFloat)  
     self assert:((-4.0 asIEEEFloat) / (-2.0 asIEEEFloat)) = (2.0 asIEEEFloat)  

     self assert:((IEEEFloat fromFloat:100.0) / (IEEEFloat fromFloat:2.0)) = (IEEEFloat fromFloat:50.0) 
     self assert:((IEEEFloat fromFloat:1e-20) / (IEEEFloat fromFloat:1e20)) = (IEEEFloat fromFloat:(1e-20 / 1e20)) 
    "
!

sumFromIEEEFloat:anIEEEFloat
    "/ anIEEEFloat + self

    |m1 m2 e1 e2 m e signBit hi shift nM1 nM2 nM|

    "/ bring to same exponent, add m1 + m2, normalize

    m1 := anIEEEFloat mantissaBits. 
    m1 == 0 ifTrue:[
        "/ always produce a positive zero
        self isZero ifTrue:[
            self negative ifTrue:[ ^ self copy basicAt:(self basicSize) put:0; yourself ].
        ].
        ^ self
    ].
    m2 := self mantissaBits.        
    m2 == 0 ifTrue:[
        anIEEEFloat isZero ifTrue:[
            anIEEEFloat negative ifTrue:[ ^ anIEEEFloat copy basicAt:(anIEEEFloat basicSize) put:0; yourself].
        ].
        ^ anIEEEFloat
    ].

    nM1 := anIEEEFloat numBitsInMantissa.
    nM2 := nM := self numBitsInMantissa.
    nM1 ~= nM2 ifTrue:[
        nM1 > nM2 ifTrue:[
            m2 := m2 bitShift:(nM1-nM2).
            nM := nM1
        ] ifFalse:[
            m1 := m1 bitShift:(nM2-nM1).
            nM := nM2
        ].
    ].

    e1 := anIEEEFloat exponentBits - anIEEEFloat eBias.
    e2 := self exponentBits - self eBias.
    e := e1.
    e1 = e2 ifFalse:[
        e1 > e2 ifTrue:[
            m2 := m2 rightShift:(e1 - e2).
            e := e1.
        ] ifFalse:[
            "/ e2 > e1
            m1 := m1 rightShift:(e2 - e1).
            e := e2.
        ]
    ].

    signBit := 0.
    "/ ok, add the mantissae
    self negative ifTrue:[
        anIEEEFloat negative ifTrue:[
            "/ - (anIEEEFloat + self)
            m := m1 + m2.
            signBit := 1.
        ] ifFalse:[
            "/ anIEEEFloat - self 
            m := m1 - m2
        ].
    ] ifFalse:[
        anIEEEFloat negative ifTrue:[
            "/ self - anIEEEFloat
            m := m2 - m1.
        ] ifFalse:[
            "/ (anIEEEFloat + self)
            m := m1 + m2.
        ].
    ].
    m == 0 ifTrue:[^ self class zero].
    m < 0 ifTrue:[
        m := m negated.
        signBit := 1-signBit
    ].

    hi := m highBit.    
    shift := (nM + 1) - hi. 
    m := m bitShift:shift.
    e := e - shift.
    m := m bitAnd:(1 bitShift:nM)-1. 

    ^ self class 
        size:(self basicSize*8) exponentSize:exponentSize 
        fraction:m exponent:(e + self eBias) signBit:signBit.

    "
     seee eeee  eeee mmmm  mmmm mmmm  mmmm mmmm  mmmm mmmm  mmmm mmmm  mmmm mmmm  mmmm mmmm 
     0100 0000  0000 0000  0000 ......                                                       2.0
     0011 1111  1111 0000  0000 ......                                                       1.0
     0400 0000  0001 0000  0000 ......                                                       4.0

     self assert:(IEEEFloat fromFloat:1.0) exponent = 1.0 exponent.  
     self assert:(IEEEFloat fromFloat:1.0) mantissaBits = 1.0 mantissaBits.
     self assert:(IEEEFloat fromFloat:1.0) eBias = 1.0 eBias.

     self assert:(IEEEFloat fromFloat:2.0) exponent = 2.0 exponent.  
     self assert:(IEEEFloat fromFloat:2.0) mantissaBits = 2.0 mantissaBits.
     self assert:(IEEEFloat fromFloat:2.0) eBias = 2.0 eBias.

     self assert:(IEEEFloat fromFloat:4.0) exponent = 4.0 exponent.  
     self assert:(IEEEFloat fromFloat:4.0) mantissaBits = 4.0 mantissaBits.
     self assert:(IEEEFloat fromFloat:4.0) eBias = 4.0 eBias.

     self assert:(IEEEFloat fromFloat:3.0) exponent = 3.0 exponent.  
     self assert:(IEEEFloat fromFloat:3.0) mantissaBits = 3.0 mantissaBits.
     self assert:(IEEEFloat fromFloat:3.0) eBias = 3.0 eBias.

     self assert:(IEEEFloat fromFloat:10.0) exponent = 10.0 exponent.  
     self assert:(IEEEFloat fromFloat:10.0) mantissaBits = 10.0 mantissaBits.
     self assert:(IEEEFloat fromFloat:10.0) eBias = 10.0 eBias.

     self assert:(IEEEFloat fromFloat:1e-3) exponent = 1e-3 exponent.  
     self assert:(IEEEFloat fromFloat:1e-3) mantissaBits = 1e-3 mantissaBits.
     self assert:(IEEEFloat fromFloat:1e-3) eBias = 1e-3 eBias.

     (IEEEFloat fromFloat:1.0) + (IEEEFloat fromFloat:1.0) = (IEEEFloat fromFloat:2.0) 
     ((IEEEFloat fromFloat:1.0) + (IEEEFloat fromFloat:1.0)) exponent = 2.0 exponent
     ((IEEEFloat fromFloat:1.0) + (IEEEFloat fromFloat:1.0)) mantissaBits = 2.0 mantissaBits

     self assert:((-4.0 asIEEEFloat) + (2.0 asIEEEFloat)) = (-2.0 asIEEEFloat)  

     (IEEEFloat fromFloat:2.0)
     self assert:((IEEEFloat fromFloat:2.0) + (IEEEFloat fromFloat:2.0)) = (IEEEFloat fromFloat:4.0) 
     self assert:((IEEEFloat fromFloat:4.0) + (IEEEFloat fromFloat:4.0)) = (IEEEFloat fromFloat:8.0) 
     self assert:((IEEEFloat fromFloat:0.0) + (IEEEFloat fromFloat:1.0)) = (IEEEFloat fromFloat:1.0) 
     self assert:((IEEEFloat fromFloat:1.0) + (IEEEFloat fromFloat:0.0)) = (IEEEFloat fromFloat:1.0) 
     self assert:((IEEEFloat fromFloat:1.0) + (IEEEFloat fromFloat:2.0)) = (IEEEFloat fromFloat:3.0) 
     self assert:((IEEEFloat fromFloat:100.0) + (IEEEFloat fromFloat:2.0)) = (IEEEFloat fromFloat:102.0) 
     self assert:((IEEEFloat fromFloat:2.0) + (IEEEFloat fromFloat:100.0)) = (IEEEFloat fromFloat:102.0) 

     self assert:((IEEEFloat fromFloat:4.0) + (IEEEFloat fromFloat:-2.0)) = (IEEEFloat fromFloat:2.0)  
     self assert:((IEEEFloat fromFloat:8.0) + (IEEEFloat fromFloat:-2.0)) = (IEEEFloat fromFloat:6.0)  
     self assert:((IEEEFloat fromFloat:100.0) + (IEEEFloat fromFloat:-2.0)) = (IEEEFloat fromFloat:98.0)  
     self assert:((IEEEFloat fromFloat:-2.0) + (IEEEFloat fromFloat:100.0)) = (IEEEFloat fromFloat:98.0)  
     self assert:((IEEEFloat fromFloat:-2.0) + (IEEEFloat fromFloat:-100.0)) = (IEEEFloat fromFloat:-102.0)  

     self assert:((IEEEFloat fromFloat:-0.1) + (IEEEFloat fromFloat:0.1)) = (IEEEFloat fromFloat:0.0)  
    "
! !

!IEEEFloat methodsFor:'mathematical functions'!

ln
    "return the natural logarithm of myself.
     Raises an exception, if the receiver is less or equal to zero.

     Not sure if this is really faster than using a taylor right away:
     the three exp-computations at the end are done in qDouble and are tailors themself..."

    |x|

    self = self class zero ifFalse:[
        self positive ifTrue:[
            "/ initial approx.
            x := self coerce:(self asFloat ln).
            "/ three more iterations of newton...
            x := x + (self * (x negated exp)) - 1.0.
            x := x + (self * (x negated exp)) - 1.0.
            x := x + (self * (x negated exp)) - 1.0.

            ^ x
        ].
    ].

    "/ now done via trapInfinity; was:
    "/ d0 = 0.0 ifTrue:[
    "/     ^ Infinity negative
    "/ ].

    "/ if you need -INF for a zero receiver, try Number trapInfinity:[...]
    ^ self class
        raise:(self = 0 ifTrue:[#infiniteResultSignal] ifFalse:[#domainErrorSignal])
        receiver:self
        selector:#ln
        arguments:#()
        errorString:'bad receiver in ln (not strictly positive)'

    "                                 
     Number trapInfinity:[ 0.0 ln ]   
     Number trapInfinity:[ 0.0 asIEEEFloat ln ]    
     1.0 asIEEEFloat ln         

     2.7 asIEEEFloat ln
     100 asIEEEFloat ln
     3.0 asIEEEFloat ln printfPrintString:'%10.8f'
    "
!

log10
    "return log base-10 of the receiver.
     Raises an exception, if the receiver is less or equal to zero.
     Here, fallback to the general logarithm code."

    ^ self ln / (self coerce:Float ln10)
! !

!IEEEFloat methodsFor:'printing & storing'!

printOn:aStream
    thisContext isRecursive ifTrue:[
        aStream nextPutAll:'IEEEFloat (recursion error while printing)'.
        ^ self.
    ].
    "/ super printOn:aStream.
    PrintfScanf printf:'%g' on:aStream argument:self.
! !

!IEEEFloat methodsFor:'queries'!

eBias
    "Answer the exponent's bias;
     that is the offset of the zero exponent when stored"

    ^ (1 bitShift:exponentSize-1) - 1

    "
     1.0 numBitsInExponent 11
     1.0 eBias             1023
     1.0 emin              -1022
     1.0 emax              1023
     1.0 fmin              2.2250738585072E-308
     1.0 fmax              1.79769313486232E+308
    "
!

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

numBitsInExponent
    ^ exponentSize

    "
     1.0 numBitsInExponent 11
     1.0 numBitsInMantissa 52
     1.0 precision         53
     1.0 decimalPrecision  16
     1.0 eBias             1023
     1.0 emin              -1022
     1.0 emax              1023
     1.0 fmin              2.2250738585072E-308
     1.0 fmax              1.79769313486232E+308
    "
!

numBitsInMantissa
    "answer the number of bits in the mantissa (the significant).
     This is an IEEE float,
        where n*8 - 1 - exponentSize (n=nr of bytes)
     bits are available in the mantissa (the hidden bit is not counted here):
        s ee...ee mmm...mmm
    "

    ^ (self basicSize * 8) - exponentSize - 1 "/ sign
!

precision
    ^ self numBitsInMantissa + 1 - self numBitsInIntegerPart 
!

radix
    "Answer the exponent's radix"

    ^ 2
! !

!IEEEFloat methodsFor:'testing'!

isFinite
    "return true, if the receiver is a finite float (not NaN or +/-INF)"

   ^ self exponentBits ~= ((1 bitShift:self numBitsInExponent)-1)

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

isInfinite
    "return true, if the receiver is an infinite float (+/-INF)"

   ^ (self exponentBits = ((1 bitShift:self numBitsInExponent)-1))
   and:[ self mantissaBits = 0 ]

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


!

isNaN
    "return true, if the receiver is an invalid float (NaN - not a number)"

   ^ (self exponentBits = ((1 bitShift:self numBitsInExponent)-1))
   and:[ self mantissaBits ~= 0 ]

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

isZero
    "return true, if the receiver is zero"

    1 to:self basicSize-1 do:[:i | 
        (self basicAt:i) ~~ 0 ifTrue:[^ false].
    ].
    ((self basicAt:self basicSize) bitAnd:16r7F) ~~ 0 ifTrue:[^ false].     
    ^ true

    "
     0.0 negated = 0.0
     0.0 negated isZero
     0.0 asIEEEFloat isZero
     0.0 asIEEEFloat negated isZero   
     0.0 asIEEEFloat = 0.0        
     0.0 asIEEEFloat negated = 0.0   
    "
!

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

    |sz hi|

    sz := self basicSize.
    ((hi := self basicAt:sz) bitTest:16r80) ifTrue:[
        "/ possibly negative; but could still be a negative zero
        hi == 16r80 ifFalse:[^ true].
        1 to:sz-1 do:[:i | (self basicAt:i) == 0 ifFalse:[^ true]].
        "/ yes - it is -0.0
    ].
    ^ false

    "
     -1.0 asIEEEFloat negative  
     0.0 asIEEEFloat negative   
     -0.0 asIEEEFloat negative  
    "
! !

!IEEEFloat methodsFor:'truncation & rounding'!

ceilingAsFloat
    "return the smallest integer-valued float greater or equal to the receiver.
     This is much like #ceiling, but avoids a (possibly expensive) conversion
     of the result to an integer.
     It may be useful, if the result is to be further used in another float-operation."

    self halt.
    ^ super ceilingAsFloat.

    "
     1.5 asIEEEFloat ceilingAsFloat
     3.0 asIEEEFloat ceilingAsFloat
    "
!

fractionPart
    "extract the after-decimal fraction part.
     such that (self truncated + self fractionPart) = self"

    |e eB n m nM|

    e := self exponentBits.
    eB := self eBias.
    e < eB ifTrue:[
        "/ integer part is zero
        ^ self
    ].
    nM := self numBitsInMantissa. 
    m := self mantissaBits.
    "/ bring in the hidden bit
    m := m bitOr:(1 bitShift:nM).

    n := e - eB + 1.
    e := eB - 1.
    m := (m bitShift:n) bitAnd:(1 bitShift:nM)-1.
    "/ normalize
    [m highBit == nM] whileTrue:[
        m := m bitShift:1.
        e := e - 1.
    ].
    ^ ((self class basicNew:self basicSize) exponentSize:exponentSize)
        setFraction:m exponent:e signBit:(self negative ifTrue:[1] ifFalse:[0]).

    "
     1.5 asIEEEFloat fractionPart
     1.0 asIEEEFloat fractionPart
     3.0 asIEEEFloat fractionPart
    "
!

truncated
    "return the receiver truncated towards zero as an integer"

    ^ super truncated.

    "
     1.5 asIEEEFloat truncated
    "
!

truncatedAsFloat
    "return the receiver truncated towards zero as a float.
     This is much like #truncated, but avoids a (possibly expensive) conversion
     of the result to an integer.
     It may be useful, if the result is to be further used in another
     float-operation."

    ^ self - self fractionPart.

    "
     1.5  truncatedAsFloat
     1.0  truncatedAsFloat
     1.5 asIEEEFloat truncatedAsFloat
     1.0 asIEEEFloat truncatedAsFloat
    "
! !

!IEEEFloat class methodsFor:'documentation'!

version_CVS
    ^ '$Header$'
! !