[squeak-dev] Exact #sqrt

Juan Vuletich juan at jvuletich.org
Wed Oct 12 01:47:46 UTC 2011 

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.

Feel free to comment, enhance, fix, etc, especially Nicolas, our 
numerics guru.

Oh, I saw the ongoing debate on automatically converting values to 
Complex, and therefore didn't include any change wrt that.

Cheers,
Juan Vuletich
-------------- next part --------------
'From Cuis 3.3 of 2 June 2011 [latest update: #1024] on 11 October 2011 at 10:27:48 pm'!

!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)
	! !


!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/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: '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/11/2011 21:53'!
nthRoot: aNumber
	"Answer the nth root of the receiver."

	self subclassResponsibility! !

!Number methodsFor: 'mathematical functions' stamp: 'jmv 10/11/2011 21:56'!
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:
		[ self error: self printString, ' raised to a non-integer power' ].
	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/11/2011 22:20'!
nthRoot: aNumber
	"Answer the nth root of the receiver."

	"Answer the nth root of the receiver."
	aNumber = 2 ifTrue: [
		^self sqrt ].
	
	^(self negative and: [ aNumber odd ])
		ifTrue: [ (self negated raisedTo: 1.0 / aNumber) negated ]
		ifFalse: [ self raisedTo: 1.0 / aNumber ]! !


!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! !


!Fraction methodsFor: 'mathematical functions' stamp: 'jmv 10/11/2011 21:54'!
nthRoot: aNumber
	"Answer the nth root of the receiver."

	^ (numerator nthRoot: aNumber) / (denominator nthRoot: aNumber)! !

!Fraction methodsFor: 'mathematical functions' stamp: 'jmv 10/11/2011 10:51'!
sqrt
	^ numerator sqrt / denominator sqrt! !


!Integer methodsFor: 'mathematical functions' stamp: 'jmv 10/11/2011 21:54'!
nthRoot: aNumber
	"Answer the nth root of the receiver."

	| floatResult guess bigIntFloat raised higher lower |
	floatResult _ self asFloat nthRoot: aNumber.
	guess _ floatResult rounded.

	"If got an exact answer, answer it. Otherwise answer float approximate answer."
	raised _ guess raisedToInteger: aNumber.
	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"
	bigIntFloat _ Float maxExactInteger.
	"avoid LargInteger multiplication if possible!!"
	(self abs > bigIntFloat and: [ self abs > (bigIntFloat raisedToInteger: aNumber) ])
		ifTrue: [
			raised > self
				ifTrue: [
					higher _ guess.
					[
						(floatResult - floatResult predecessor) rounded.
						floatResult _ floatResult predecessor.
						lower _ floatResult rounded.
						(lower raisedToInteger: aNumber) > self ] whileTrue ]
				ifFalse: [
					lower _ guess.
					[
						(floatResult successor - floatResult) rounded.
						floatResult _ floatResult successor.
						higher _ floatResult rounded.
						(higher raisedToInteger: aNumber) < self ] whileTrue ].
			[ higher - lower > 1 ] whileTrue: [
				guess _ lower + higher // 2.
				raised _ guess raisedToInteger: aNumber.
				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/11/2011 13:12'!
sqrt
	"Answer the square root of the receiver."

	| floatResult integerResult bigIntFloat |
	floatResult _ self asFloat 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 gets inexact, even if we are a whole square number"
	bigIntFloat _ Float maxExactInteger.
	"avoid LargInteger multiplication if possible!!"
	(self > bigIntFloat and: [ self > bigIntFloat squared ])
		ifTrue: [
			integerResult _ self sqrtFloor.
			integerResult * integerResult = self ifTrue: [
				^integerResult ]].

	"We need an approximate result"
	^floatResult! !


!FractionTest reorganize!
('tests - printing' testFractionPrinting)
('tests - sinuses' testDegreeCos testDegreeSin testReciprocal)
('tests - mathematical functions' testExactRaisedTo testExactSqrt testLn testLog)
('private' assert:classAndValueEquals:)
!

-------------- next part --------------
A non-text attachment was scrubbed...
Name: 1105-sqrtAndnThRoot-jmv.4.cs.zip
Type: application/zip
Size: 2850 bytes
Desc: not available
Url : http://lists.squeakfoundation.org/pipermail/squeak-dev/attachments/20111011/dab59591/1105-sqrtAndnThRoot-jmv.4.cs.zip

More information about the Squeak-dev mailing list