[squeak-dev] Exact #sqrt

Nicolas Cellier nicolas.cellier.aka.nice at gmail.com
Thu Oct 13 20:48:54 UTC 2011


2011/10/13 Juan Vuletich <juan at jvuletich.org>:
> Juan Vuletich wrote:
>>
>> Hi Folks,
>>
>> Yesterday I realized that there's no good reason for 8 sqrt to answer 2.0
>> and not just 2. To fix this inelegant situation I wrote the attached. This
>> way, for instance, 4 sqrt = (8 raisedTo: 1/3) answers true. I also
>> introduced #nthRoot: , which can be of some use. A rather complete set of
>> tests is included, to illustrate the desired behavior.
>>
>> ...
>
> This new version includes performance improvements and better checking of
> valid domains. Enjoy!
>
> Cheers,
> Juan Vuletich
>
> 'From Cuis 3.3 of 2 June 2011 [latest update: #1024] on 13 October 2011 at
> 9:18:36 am'!
>
> !FloatTest methodsFor: 'testing' stamp: 'jmv 10/11/2011 08:55'!
> testMaxExactInteger
>        "
>        FloatTest new testMaxExactInteger
>        "
>
>        self assert: Float maxExactInteger asFloat truncated = Float
> maxExactInteger.
>        0 to: 10000 do: [ :j |
>                self assert: (Float maxExactInteger-j) asFloat truncated =
> (Float maxExactInteger-j) ].
>        self deny: (Float maxExactInteger+1) asFloat truncated = (Float
> maxExactInteger+1)
>        ! !
>
> !FloatTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/13/2011 09:06'!
> testNthRoot
>        "
>        FloatTest new testNthRoot
>        "
>        self should: [ -1.23 nthRoot: 4 ] raise: ArithmeticError! !
>
> !FloatTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/13/2011 09:06'!
> testRaisedTo
>        "
>        FloatTest new testRaisedTo
>        "
>        self should: [ -1.23 raisedTo: 1/4 ] raise: ArithmeticError! !
>
>
> !FractionTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/11/2011 22:27'!
> testExactRaisedTo
>        "
>        FractionTest new testExactRaisedTo
>        "
>        | f |
>        self assert: (4/9 raisedTo: 1/2) classAndValueEquals: 2/3.
>        self assert: (9/4 raisedTo: 1/2) classAndValueEquals: 3/2.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) pairsDo: [ :a :b |
>                f _ a / b.
>                self assert: (f squared raisedTo: 1/2) classAndValueEquals:
> f.
>                self assert: (f negated squared raisedTo: 1/2)
> classAndValueEquals: f.
>                f _ b / a.
>                self assert: (f squared raisedTo: 1/2) classAndValueEquals:
> f.
>                self assert: (f negated squared raisedTo: 1/2)
> classAndValueEquals: f ].
>
>        self assert: (8/27 raisedTo: 1/3) classAndValueEquals: 2/3.
>        self assert: (27/8 raisedTo: 1/3) classAndValueEquals: 3/2.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) pairsDo: [ :a :b |
>                f _ a / b.
>                self assert: ((f raisedTo: 3) raisedTo: 1/3)
> classAndValueEquals: f.
>                self assert: ((f negated raisedTo: 3) raisedTo: 1/3)
> classAndValueEquals: f negated.
>                f _ b / a.
>                self assert: ((f raisedTo: 3) raisedTo: 1/3)
> classAndValueEquals: f.
>                self assert: ((f negated raisedTo: 3) raisedTo: 1/3)
> classAndValueEquals: f negated ].
>
>        self assert: (4/9 raisedTo: 3/2) classAndValueEquals: 8/27.
>        self assert: (8/27 raisedTo: 2/3) classAndValueEquals: 4/9.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) pairsDo: [ :a :b |
>                f _ a / b.
>                self assert: ((f raisedTo: 3) raisedTo: 2/3)
> classAndValueEquals: f*f.
>                self assert: ((f raisedTo: 2) raisedTo: 3/2)
> classAndValueEquals: f*f*f.
>                self assert: ((f negated raisedTo: 3) raisedTo: 2/3)
> classAndValueEquals: f*f.
>                self assert: ((f negated raisedTo: 2) raisedTo: 3/2)
> classAndValueEquals: f*f*f.
>                f _ b / a.
>                self assert: ((f raisedTo: 3) raisedTo: 2/3)
> classAndValueEquals: f*f.
>                self assert: ((f raisedTo: 2) raisedTo: 3/2)
> classAndValueEquals: f*f*f.
>                self assert: ((f negated raisedTo: 3) raisedTo: 2/3)
> classAndValueEquals: f*f.
>                self assert: ((f negated raisedTo: 2) raisedTo: 3/2)
> classAndValueEquals: f*f*f ].
>
>        self assert: (32/243 raisedTo: 3/5) classAndValueEquals: 8/27.
>        self assert: (8/27 raisedTo: 5/3) classAndValueEquals: 32/243.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) pairsDo: [ :a :b |
>                f _ a / b.
>                self assert: ((f raisedTo: 5) raisedTo: 3/5)
> classAndValueEquals: f*f*f.
>                self assert: ((f raisedTo: 3) raisedTo: 5/3)
> classAndValueEquals: f*f*f*f*f.
>                self assert: ((f negated raisedTo: 5) raisedTo: 3/5)
> classAndValueEquals: (f*f*f) negated.
>                self assert: ((f negated raisedTo: 3) raisedTo: 5/3)
> classAndValueEquals: (f*f*f*f*f) negated.
>
>                self assert: ((f raisedTo: -5) raisedTo: 3/5)
> classAndValueEquals: 1/(f*f*f).
>                self assert: ((f raisedTo: -3) raisedTo: 5/3)
> classAndValueEquals: 1/(f*f*f*f*f).
>                self assert: ((f negated raisedTo: -5) raisedTo: 3/5)
> classAndValueEquals: -1/(f*f*f).
>                self assert: ((f negated raisedTo: -3) raisedTo: 5/3)
> classAndValueEquals: -1/(f*f*f*f*f).
>                self assert: ((f raisedTo: 5) raisedTo: -3/5)
> classAndValueEquals: 1/(f*f*f).
>                self assert: ((f raisedTo: 3) raisedTo: -5/3)
> classAndValueEquals: 1/(f*f*f*f*f).
>                self assert: ((f negated raisedTo: 5) raisedTo: -3/5)
> classAndValueEquals: -1/(f*f*f).
>                self assert: ((f negated raisedTo: 3) raisedTo: -5/3)
> classAndValueEquals: -1/(f*f*f*f*f).
>
>                "No exact result => Float result"
>                self assert: ((f raisedTo: 3) +1 raisedTo: 5/3) isFloat.
>                self assert: ((f negated raisedTo: 3) -1 raisedTo: 5/3)
> isFloat.
>
>                f _ b / a.
>                self assert: ((f raisedTo: 5) raisedTo: 3/5)
> classAndValueEquals: f*f*f.
>                self assert: ((f raisedTo: 3) raisedTo: 5/3)
> classAndValueEquals: f*f*f*f*f.
>                self assert: ((f negated raisedTo: 5) raisedTo: 3/5)
> classAndValueEquals: (f*f*f) negated.
>                self assert: ((f negated raisedTo: 3) raisedTo: 5/3)
> classAndValueEquals: (f*f*f*f*f) negated.
>
>                "No exact result => Float result"
>                self assert: ((f raisedTo: 3) +1 raisedTo: 5/3) isFloat.
>                self assert: ((f negated raisedTo: 3) -1 raisedTo: 5/3)
> isFloat ].! !
>
> !FractionTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/11/2011 22:12'!
> testExactSqrt
>        "
>        FractionTest new testExactSqrt
>        "
>        | f |
>        self assert: (4/9) sqrt classAndValueEquals: 2/3.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) pairsDo: [ :i :j |
>                f _ i / j.
>                self assert: f squared sqrt classAndValueEquals: f.
>                f _ j / i.
>                self assert: f squared sqrt classAndValueEquals: f ]! !
>
> !FractionTest methodsFor: 'private' stamp: 'jmv 10/11/2011 22:12'!
> assert: a classAndValueEquals: b
>        self assert: a class = b class.
>        self assert: a = b! !
>
>
> !IntegerTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/11/2011 22:27'!
> testExactRaisedTo
>        "
>        IntegerTest new testExactRaisedTo
>        "
>        self assert: (4 raisedTo: 1/2) classAndValueEquals: 2.
>        self assert: (9 raisedTo: 1/2) classAndValueEquals: 3.
>        self assert: (9 raisedTo: -1/2) classAndValueEquals: 1/3.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) do: [ :i |
>                self assert: (i squared raisedTo: 1/2) classAndValueEquals:
> i.
>                self assert: (i negated squared raisedTo: 1/2)
> classAndValueEquals: i ].
>
>        self assert: (8 raisedTo: 1/3) classAndValueEquals: 2.
>        self assert: (27 raisedTo: 1/3) classAndValueEquals: 3.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) do: [ :i |
>                self assert: ((i raisedTo: 3) raisedTo: 1/3)
> classAndValueEquals: i.
>                self assert: ((i negated raisedTo: 3) raisedTo: 1/3)
> classAndValueEquals: i negated ].
>
>        self assert: (4 raisedTo: 3/2) classAndValueEquals: 8.
>        self assert: (8 raisedTo: 2/3) classAndValueEquals: 4.
>        self assert: (8 raisedTo: -2/3) classAndValueEquals: 1/4.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) do: [ :i |
>                self assert: ((i raisedTo: 3) raisedTo: 2/3)
> classAndValueEquals: i*i.
>                self assert: ((i raisedTo: 2) raisedTo: 3/2)
> classAndValueEquals: i*i*i.
>                self assert: ((i negated raisedTo: 3) raisedTo: 2/3)
> classAndValueEquals: i*i.
>                self assert: ((i negated raisedTo: 2) raisedTo: 3/2)
> classAndValueEquals: i*i*i ].
>
>        self assert: (32 raisedTo: 3/5) classAndValueEquals: 8.
>        self assert: (8 raisedTo: 5/3) classAndValueEquals: 32.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) do: [ :i |
>                self assert: ((i raisedTo: 5) raisedTo: 3/5)
> classAndValueEquals: i*i*i.
>                self assert: ((i raisedTo: 3) raisedTo: 5/3)
> classAndValueEquals: i*i*i*i*i.
>                self assert: ((i negated raisedTo: 5) raisedTo: 3/5)
> classAndValueEquals: (i*i*i) negated.
>                self assert: ((i negated raisedTo: 3) raisedTo: 5/3)
> classAndValueEquals: (i*i*i*i*i) negated.
>
>                self assert: ((i raisedTo: -5) raisedTo: 3/5)
> classAndValueEquals: 1/(i*i*i).
>                self assert: ((i raisedTo: -3) raisedTo: 5/3)
> classAndValueEquals: 1/(i*i*i*i*i).
>                self assert: ((i negated raisedTo: -5) raisedTo: 3/5)
> classAndValueEquals: -1/(i*i*i).
>                self assert: ((i negated raisedTo: -3) raisedTo: 5/3)
> classAndValueEquals: -1/(i*i*i*i*i).
>
>                self assert: ((i raisedTo: 5) raisedTo: -3/5)
> classAndValueEquals: 1/(i*i*i).
>                self assert: ((i raisedTo: 3) raisedTo: -5/3)
> classAndValueEquals: 1/(i*i*i*i*i).
>                self assert: ((i negated raisedTo: 5) raisedTo: -3/5)
> classAndValueEquals: -1/(i*i*i).
>                self assert: ((i negated raisedTo: 3) raisedTo: -5/3)
> classAndValueEquals: -1/(i*i*i*i*i).
>
>                "No exact result => Float result"
>                self assert: ((i raisedTo: 3) +1 raisedTo: 5/3) isFloat.
>                self assert: ((i negated raisedTo: 3) -1 raisedTo: 5/3)
> isFloat ].! !
>
> !IntegerTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/13/2011 09:07'!
> testExactRaisedToErrorConditions
>        "
>        IntegerTest new testExactRaisedToErrorConditions
>        "
>
>        self should: [ -2 raisedTo: 1/4 ] raise: ArithmeticError.
>        self should: [ -2 raisedTo: 1.24 ] raise: ArithmeticError.! !
>
> !IntegerTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/11/2011 22:09'!
> testExactSqrt
>        "
>        IntegerTest new testExactSqrt
>        "
>        self assert: 4 sqrt classAndValueEquals: 2.
>        self assert: 9 sqrt classAndValueEquals: 3.
>        self assert: Float maxExactInteger squared sqrt classAndValueEquals:
> Float maxExactInteger.
>        self assert: (Float maxExactInteger+1) squared sqrt
> classAndValueEquals: Float maxExactInteger+1.
>        #( 1 5 29 135 1234 567890 123123123 456456456456
> 98765432109876543210987654321
> 987123987123987123987123987123987123987123987123) do: [ :i |
>                self assert: i squared sqrt classAndValueEquals: i ]! !
>
> !IntegerTest methodsFor: 'tests - mathematical functions' stamp: 'jmv
> 10/13/2011 09:09'!
> testNthRootErrorConditions
>        "
>        IntegerTest new testExactRaisedToErrorConditions
>        "
>
>        self should: [ -2 nthRoot: 1/4 ] raise: ArithmeticError.
>        self should: [ -2 nthRoot: 1.24 ] raise: ArithmeticError.! !
>
> !IntegerTest methodsFor: 'private' stamp: 'jmv 10/11/2011 08:14'!
> assert: a classAndValueEquals: b
>        self assert: a class = b class.
>        self assert: a = b! !
>
>
> !Number methodsFor: 'mathematical functions' stamp: 'jmv 10/13/2011 08:36'!
> nthRoot: aPositiveInteger
>        "Answer the nth root of the receiver."
>
>        self subclassResponsibility! !
>
> !Number methodsFor: 'mathematical functions' stamp: 'jmv 10/13/2011 08:50'!
> raisedTo: aNumber
>        "Answer the receiver raised to aNumber."
>
>        aNumber isInteger ifTrue: [
>                "Do the special case of integer power"
>                ^ self raisedToInteger: aNumber].
>        aNumber isFraction ifTrue: [
>                "Special case for fraction power"
>                ^ (self nthRoot: aNumber denominator) raisedToInteger:
> aNumber numerator ].
>        self < 0 ifTrue: [
>                ^ ArithmeticError signal: 'Negative numbers can''t be raised
> to float powers.' ].
>        0 = aNumber ifTrue: [^ self class one]. "Special case of exponent=0"
>        1 = aNumber ifTrue: [^ self].   "Special case of exponent=1"
>        0 = self ifTrue: [                              "Special case of self
> = 0"
>                aNumber < 0
>                        ifTrue: [^ (ZeroDivide dividend: self) signal]
>                        ifFalse: [^ self]].
>        ^ (aNumber * self ln) exp               "Otherwise use logarithms"! !
>
> !Number methodsFor: 'mathematical functions' stamp: 'jmv 10/11/2011 08:34'!
> sqrt
>        "Answer the square root of the receiver."
>
>        self subclassResponsibility! !
>
>
> !Float methodsFor: 'mathematical functions' stamp: 'jmv 10/13/2011 09:03'!
> nthRoot: aPositiveInteger
>        "Answer the nth root of the receiver."
>        aPositiveInteger = 2 ifTrue: [
>                ^self sqrt ].
>
>        (aPositiveInteger isInteger not or: [ aPositiveInteger negative ])
>                ifTrue: [^ ArithmeticError signal: 'nth root only defined for
> positive Integer n.'].
>
>        ^self negative
>                ifTrue: [
>                        aPositiveInteger odd
>                                ifTrue: [ (self negated raisedTo: 1.0 /
> aPositiveInteger) negated ]
>                                ifFalse: [ ArithmeticError signal: 'Negative
> numbers don''t have even roots.' ]]
>                ifFalse: [ self raisedTo: 1.0 / aPositiveInteger ]! !
>
>
> !Float class methodsFor: 'constants' stamp: 'jmv 10/11/2011 08:54'!
> maxExactInteger
>        "Answer the biggest integer such that it is exactly represented in a
> float, and all smaller integers also are"
>        ^1 bitShift: 53! !

or         ^1 bitShift: self precision


>
>
> !Fraction methodsFor: 'mathematical functions' stamp: 'jmv 10/13/2011
> 08:36'!
> nthRoot: aPositiveInteger
>        "Answer the nth root of the receiver."
>
>        ^ (numerator nthRoot: aPositiveInteger) / (denominator nthRoot:
> aPositiveInteger)! !
>
> !Fraction methodsFor: 'mathematical functions' stamp: 'jmv 10/11/2011
> 10:51'!
> sqrt
>        ^ numerator sqrt / denominator sqrt! !
>

This can cause an overflow...

((1 << 1024 + 1) / (1 << 1024 + 3)) sqrt

I'd expect 1.0 (inexact, but not that much).

>
> !Integer methodsFor: 'mathematical functions' stamp: 'jmv 10/13/2011 09:18'!
> nthRoot: aPositiveInteger
>        "Answer the nth root of the receiver."
>
>        | selfAsFloat floatResult guess raised higher lower delta |
>        selfAsFloat _ self asFloat.
>        floatResult _ selfAsFloat nthRoot: aPositiveInteger.
>        guess _ floatResult rounded.
>
>        "If got an exact answer, answer it."
>        raised _ guess raisedToInteger: aPositiveInteger.
>        raised = self
>                ifTrue: [ ^ guess ].
>
>        "In this case, maybe it failed because we are such a big integer that
> the Float
>        method gets inexact, even if we are a whole square number.
>        Note (jmv): The algorithms I found for computing the nthRoot would
> havily use
>        very large fractions. I wrote this one, that doesn't create
> fractions."
>        selfAsFloat abs >= (Float maxExactInteger asFloat raisedToInteger:
> aPositiveInteger)
>                ifTrue: [
>                        raised > self
>                                ifTrue: [
>                                        higher _ guess.
>                                        delta _  floatResult predecessor -
> floatResult.
>                                        [
>                                                floatResult _ floatResult +
> delta.
>                                                lower _ floatResult rounded.
>                                                (lower raisedToInteger:
> aPositiveInteger) > self ] whileTrue: [
>                                                        delta _ delta * 2.
>                                                        higher _ lower ] ]
>                                ifFalse: [
>                                        lower _ guess.
>                                        delta _  floatResult successor -
> floatResult.
>                                        [
>                                                floatResult _ floatResult +
> delta.
>                                                higher _ floatResult rounded.
>                                                (higher raisedToInteger:
> aPositiveInteger) < self ] whileTrue: [
>                                                        delta _ delta * 2.
>                                                        lower _ higher ]].
>                        [ higher - lower > 1 ] whileTrue: [
>                                guess _ lower + higher // 2.
>                                raised _ guess raisedToInteger:
> aPositiveInteger.
>                                raised = self
>                                        ifTrue: [
>                                                ^ guess ].
>                                raised > self
>                                        ifTrue: [ higher _ guess ]
>                                        ifFalse: [ lower _ guess ]]].
>
>        "We need an approximate result"
>        ^floatResult! !
>
> !Integer methodsFor: 'mathematical functions' stamp: 'jmv 10/12/2011 15:17'!
> sqrt
>        "Answer the square root of the receiver."
>
>        | selfAsFloat floatResult integerResult |
>        selfAsFloat _ self asFloat.
>        floatResult _ selfAsFloat sqrt.
>        integerResult _ floatResult truncated.
>
>        "If got an exact answer, answer it. Otherwise answer float
> approximate answer."
>        integerResult squared = self
>                ifTrue: [ ^ integerResult ].
>
>        "In this case, maybe it failed because we are such a big integer that
> the Float method becomes
>        inexact, even if we are a whole square number. So, try the slower but
> more general method"
>        selfAsFloat >= Float maxExactInteger asFloat squared
>                ifTrue: [
>                        integerResult _ self sqrtFloor.
>                        integerResult * integerResult = self ifTrue: [
>                                ^integerResult ]].
>
>        "We need an approximate result"
>        ^floatResult! !
>
>
> !FloatTest reorganize!
> ('testing - arithmetic' testDivide)
> ('infinity behavior' testHugeIntegerCloseTo testInfinity1 testInfinity2
> testInfinityCloseTo)
> ('IEEE 754' test32bitGradualUnderflow test32bitRoundingMode testInfinity3
> testNaN5 testZero2)
> ('zero behavior' testIsZero testNegativeZeroAbs testNegativeZeroSign
> testZero1 testZeroSignificandAsInteger)
> ('NaN behavior' testNaN1 testNaN2 testNaN3 testNaN4 testNaNCompare
> testNaNisLiteral testReciprocal)
> ('testing - conversion' testAsTrueFraction testFloatRounded
> testFloatTruncated testFractionAsFloat testFractionAsFloat2
> testIntegerAsFloat testStringAsNumber)
> ('testing compare' testCloseTo testComparison
> testComparisonWhenPrimitiveFails)
> ('printing' testStoreBase16)
> ('testing' testArcTan testCopy testMaxExactInteger testSetOfFloat
> testStoreOn)
> ('characterization' testCharacterization)
> ('tests - mathematical functions' testNthRoot testRaisedTo)
> !
>
>
>
>
>



More information about the Squeak-dev mailing list