[squeak-dev] Exact #sqrt
Nicolas Cellier
nicolas.cellier.aka.nice at gmail.com
Thu Oct 13 23:37:52 UTC 2011
You might also want another way to quickly discriminate non square
large integers
http://en.wikipedia.org/wiki/Square_number
if (n & 7 == 1) or (n & 31 == 4) or (n & 127 == 16) or (n & 191 == 0) then
return n is probably square
else
return n is definitely not square
LargePositiveInteger>>canBeASquare
| lsb |
lsb := self digitAt: 1.
^((lsb bitAnd: 7) = 1
or: [(lsb bitAnd: 31) = 4
or: [(lsb bitAnd: 127) = 16
or: [(lsb bitAnd: 191) = 0]]])
Nicolas
2011/10/13 Nicolas Cellier <nicolas.cellier.aka.nice at gmail.com>:
> 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
|