[squeak-dev] 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 Squeak-dev mailing list