[squeakdev] Exact #sqrt
Juan Vuletich
juan at jvuletich.org
Tue Oct 18 01:06:32 UTC 2011
Nicolas Cellier wrote:
> Hi Juan,
> Finally, why didn't you just tried some sort of Newton ?
> ... great code snipped...
> It's not exactly cheap (repeated raisedTo:) but it should work without Fraction.
> Note that selector is named Truncated not Floor, because it better
> reflects the operation on negative receivers.
>
> Nicolas
>
Good idea! BTW, in #nthRootTruncated you send #nthRootFloor:. I fixed it.
To put it to good use, I wrote a new #nthRoot: that calls it. My
previous implementation was moved to #nthRootAlt:. The new
implementation is better because it is not limited to the Float range.
This is important! OTOH, #nthRootAlt: seems to be faster for n>=9. In a
comment there is the code I used to test this.
The attached also includes a fix to a comment, and a new test.
Cheers,
Juan Vuletich
 next part 
'From Cuis 3.3 of 2 June 2011 [latest update: #1024] on 17 October 2011 at 10:50:43 pm'!
!Fraction methodsFor: 'mathematical functions' stamp: 'jmv 10/17/2011 21:52'!
nthRoot: aPositiveInteger
"Answer the nth root of the receiver."
 d n 
n _numerator nthRoot: aPositiveInteger.
d _ denominator nthRoot: aPositiveInteger.
"The #nthRoot: method in integer will only answer a Float if there's no exact nth root.
So, we need a float anyway."
(n isInfinite or: [ d isInfinite ]) ifTrue: [
^self asFloat nthRoot: aPositiveInteger ].
^n / d! !
!Integer methodsFor: 'mathematical functions' stamp: 'jmv 10/17/2011 22:42'!
nthRoot: aPositiveInteger
"Answer the nth root of the receiver.
See #nthRootAlt: for an alternative implementation."
 selfAsFloat floatResult guess 
selfAsFloat _ self asFloat.
floatResult _ selfAsFloat nthRoot: aPositiveInteger.
floatResult isInfinite ifFalse: [
guess _ floatResult rounded.
"If got an exact answer, answer it."
(guess raisedToInteger: aPositiveInteger) = self
ifTrue: [ ^ guess ]].
"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 abs >= (Float maxExactInteger asFloat raisedToInteger: aPositiveInteger)
ifTrue: [
guess _ self nthRootTruncated: aPositiveInteger.
(guess raisedToInteger: aPositiveInteger) = self
ifTrue: [ ^ guess ]].
"We need an approximate result"
^floatResult! !
!Integer methodsFor: 'mathematical functions' stamp: 'jmv 10/17/2011 22:42'!
nthRootAlt: aPositiveInteger
"Answer the nth root of the receiver.
Alternative implementation. Seems faster than #nthRoot for n >=9
(See test code at the end of the method)
Note, however, that this method is limited to the Float range, while #nthRoot can handle much bigger integers"
 selfAsFloat floatResult guess raised higher lower delta 
selfAsFloat _ self asFloat.
floatResult _ selfAsFloat nthRoot: aPositiveInteger.
"If we can't do Float arithmetic, currently we can't look for an exact answer"
floatResult isInfinite ifTrue: [
^floatResult ].
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
"
 i n 
i _ 1234987687234509123.
n _ i raisedTo: 3.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 3]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 3]] timeToRun print.
n _ i raisedTo: 5.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 5]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 5]] timeToRun print.
n _ i raisedTo: 7.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 7]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 7]] timeToRun print.
n _ i raisedTo: 9.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 9]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 9]] timeToRun print.
n _ i raisedTo: 11.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 11]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 11]] timeToRun print.
n _ i raisedTo: 13.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 13]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 13]] timeToRun print.
n _ i raisedTo: 15.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRoot: 15]] timeToRun print.
Smalltalk garbageCollect.
[ 10000 timesRepeat: [n nthRootAlt: 15]] timeToRun print.
"! !
!Integer methodsFor: 'mathematical functions' stamp: 'jmv 10/17/2011 21:47'!
nthRootTruncated: aPositiveInteger
"Answer the integer part of the nth root of the receiver."
 guess guessToTheNthMinusOne delta 
self = 0 ifTrue: [^0].
self negative
ifTrue:
[aPositiveInteger even ifTrue: [ ArithmeticError signal: 'Negative numbers don''t have even roots.' ].
^(self negated nthRootTruncated: aPositiveInteger) negated].
guess := 1 bitShift: self highBitOfMagnitude + aPositiveInteger  1 // aPositiveInteger.
[
guessToTheNthMinusOne := guess raisedTo: aPositiveInteger  1.
delta := (guess * guessToTheNthMinusOne  self) // (guessToTheNthMinusOne * aPositiveInteger).
delta = 0 ] whileFalse:
[ guess := guess  delta ].
( (guess := guess  1) raisedTo: aPositiveInteger) > self ifTrue:
[ guess := guess  1 ].
^guess! !
!Integer methodsFor: 'mathematical functions' stamp: 'nice 10/17/2011 22:53'!
sqrtFloor
"Return the integer part of the square root of self"
 guess delta 
guess := 1 bitShift: self highBit + 1 // 2.
[
delta := guess squared  self // (guess bitShift: 1).
delta = 0 ] whileFalse: [
guess := guess  delta ].
^guess  1! !
!IntegerTest methodsFor: 'tests  mathematical functions' stamp: 'jmv 10/17/2011 22:47'!
testNthRoot
"
IntegerTest new testNthRoot
"
 i 
i _ 1234987687234509123.
#(3 5 7 9 11 13 15 17 19 21 23 25 27) do: [ :n 
self assert: ((i raisedTo: n) nthRoot: n) = i ].
#(3 5 7 9 11 13 15 17 ) do: [ :n 
self assert: ((i raisedTo: n) nthRootAlt: n) = i ]! !
!LargePositiveInteger methodsFor: 'mathematical functions' stamp: 'nice 10/14/2011 23:50'!
sqrtFloor
"Return the integer part of the square root of self"
 powerOfTwo 
(powerOfTwo := self lowBit  1 // 2) > 1
ifFalse: [^super sqrtFloor].
^(self bitShift: 2 * powerOfTwo) sqrtFloor bitShift: powerOfTwo! !
More information about the Squeakdev
mailing list
