'From Cuis 2.6 of 10 August 2010 [latest update: #540] on 1 September 2010 at 9:53:01 am'!

!Number methodsFor: 'mathematical functions' stamp: 'nice 12/6/2007 21:46'!
raisedTo: aNumber 
	"Answer the receiver raised to aNumber."

	aNumber isInteger ifTrue:
		["Do the special case of integer power"
		^ self raisedToInteger: aNumber].
	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"! !


!Float methodsFor: 'mathematical functions' stamp: 'nice 10/30/2009 22:21'!
arcTan: denominator
	"Answer the angle in radians.
	 Optional. See Object documentation whatIsAPrimitive.
	Implementation note: use sign in order to catch cases of negativeZero"

	^self = 0.0
		ifTrue: [denominator sign >= 0
			ifTrue: [ 0 ]
			ifFalse: [ self sign >= 0
				ifTrue: [ Pi ]
				ifFalse: [ Pi negated ]]]
		ifFalse: [denominator = 0.0
			ifTrue: [self > 0.0
				ifTrue: [ Halfpi ]
				ifFalse: [ Halfpi negated ]]
			ifFalse: [denominator > 0
				ifTrue: [ (self / denominator) arcTan ]
				ifFalse: [self > 0
					ifTrue: [ ((self / denominator) arcTan) + Pi ]
					ifFalse: [ ((self / denominator) arcTan) - Pi ]]]]! !

!Float methodsFor: 'mathematical functions' stamp: 'jmv 8/31/2010 19:04'!
sqrt
	"Answer the square root of the receiver. 
	 Optional. See Object documentation whatIsAPrimitive."
	| exp guess eps delta |
	<primitive: 55>

	"Newton-Raphson"
	self <= 0.0
		ifTrue: [self = 0.0
				ifTrue: [^ 0.0]
				ifFalse: [
					^ FloatingPointException signal: 'undefined if less than zero.']].
	"first guess is half the exponent"
	exp := self exponent // 2.
	guess := self timesTwoPower: 0 - exp.
	"get eps value"
	eps := guess * Epsilon.
	eps := eps * eps.
	delta := self - (guess * guess) / (guess * 2.0).
	[delta * delta > eps]
		whileTrue: 
			[guess := guess + delta.
			delta := self - (guess * guess) / (guess * 2.0)].
	^ guess! !

!Float methodsFor: 'comparing' stamp: 'eem 7/6/2009 11:11'!
>= aNumber 
	"Primitive. Compare the receiver with the argument and return true
	if the receiver is greater than or equal to the argument. Otherwise return
	false. Fail if the argument is not a Float. Optional. See Object documentation 
	whatIsAPrimitive. "

	<primitive: 46>
	^ aNumber adaptToFloat: self andSend: #>=! !

!Float methodsFor: 'comparing' stamp: 'nice 7/19/2009 19:27'!
closeTo: num
 	"are these two numbers close?"
	num isNumber ifFalse: [^[self = num] ifError: [false]].
	self = 0.0 ifTrue: [^num abs < 0.0001].
	num = 0 ifTrue: [^self abs < 0.0001].
	^self = num asFloat
		or: [(self - num) abs / (self abs max: num abs) < 0.0001]! !

!Float methodsFor: 'comparing' stamp: 'nice 1/4/2009 18:02'!
hash
	"Hash is reimplemented because = is implemented. Both words of the float are used. (The bitShift:'s ensure that the intermediate results do not become a large integer.) Care is taken to answer same hash as an equal Integer."

	(self isFinite and: [self fractionPart = 0.0]) ifTrue: [^self truncated hash].
	^ ((self basicAt: 1) bitShift: -4) +
	   ((self basicAt: 2) bitShift: -4)
! !

!Float methodsFor: 'testing' stamp: 'nice 3/14/2008 23:45'!
isFinite
	"simple, byte-order independent test for rejecting Not-a-Number and (Negative)Infinity"

	^(self - self) = 0.0! !

!Float methodsFor: 'testing' stamp: 'nice 3/14/2008 23:49'!
isLiteral
	"There is no literal representation of NaN.
	However, there are literal representations of Infinity, like 1.0e1000.
	But since they are not able to print properly, only case of finite Float is considered."
	
	^self isFinite! !

!Float methodsFor: 'testing' stamp: 'jmv 8/31/2010 00:00'!
sign
	"Answer 1 if the receiver is greater than 0, -1 if less than 0, else 0.
	Handle IEEE-754 negative-zero by reporting a sign of -1
	Warning!! This makes Float negativeZero the only number in the system such that
		x sign negated = x negated sign
	evaluates to false!!
	This precludes the simpler implementation in #signPart
	0.0 sign  ->  0
	0.0 signPart   ->  1
	Float negativeZero sign   ->  -1
	Float negativeZero signPart   ->  -1
	"

	self > 0 ifTrue: [^ 1].
	(self < 0 or: [((self at: 1) bitShift: -31) = 1]) ifTrue: [^ -1].
	^ 0! !

!Float methodsFor: 'truncation and round off' stamp: 'jmv 8/30/2010 23:39'!
exponentPart
	"Alternative implementation for exponent"
	^self partValues: [ :sign :exponent :mantissa | exponent ]! !

!Float methodsFor: 'truncation and round off' stamp: 'jmv 8/30/2010 23:40'!
mantissaPart
	"Equivalent to #significandAsInteger, but nicer implementation."
	^self partValues: [ :sign :exponent :mantissa | mantissa ]! !

!Float methodsFor: 'truncation and round off' stamp: 'nice 7/24/2008 01:32'!
rounded
	"Answer the integer nearest the receiver.
	Implementation note: super would not handle tricky inexact arithmetic"
	
	"self assert: 5000000000000001.0 rounded = 5000000000000001"

	self fractionPart abs < 0.5
		ifTrue: [^self truncated]
		ifFalse: [^self truncated + self sign rounded]! !

!Float methodsFor: 'truncation and round off' stamp: 'jmv 8/30/2010 23:59'!
signPart
	"Alternative implementation for exponent"
	^self partValues: [ :sign :exponent :mantissa | sign ]! !

!Float methodsFor: 'truncation and round off' stamp: 'jmv 8/30/2010 23:49'!
significandAsInteger

	^self partValues: [ :sign :exponent :mantissa | mantissa ]! !

!Float methodsFor: 'truncation and round off' stamp: 'jmv 8/31/2010 14:47'!
truncated
	"Answer with a SmallInteger equal to the value of the receiver without 
	its fractional part. The primitive fails if the truncated value cannot be 
	represented as a SmallInteger. In that case, the code below will compute 
	a LargeInteger truncated value.
	Essential. See Object documentation whatIsAPrimitive. "

	<primitive: 51>

	^ self partValues: [ :sign :exponent :mantissa |
		sign * (mantissa bitShift: exponent - 52) ]! !

!Float methodsFor: 'converting' stamp: 'nice 5/30/2006 02:29'!
asIEEE32BitWord
	"Convert the receiver into a 32 bit Integer value representing the same number in IEEE 32 bit format.
	Used for conversion in FloatArrays only."
	
	| word1 word2 sign mantissa exponent destWord truncatedBits mask roundToUpper |
	
	"skip fast positive and nnegative zero"
	self = 0.0 ifTrue: [^self basicAt: 1].
	
	"retrieve 64 bits of IEEE 754 double"
	word1 := self basicAt: 1.
	word2 := self basicAt: 2.
	
	"prepare sign exponent and mantissa of 32 bits float"
	sign := word1 bitAnd: 16r80000000.
	exponent := ((word1 bitShift: -20) bitAnd: 16r7FF) - 1023 + 127.
	mantissa := (word2 bitShift: -29) + ((word1 bitAnd:  16rFFFFF) bitShift: 3).
	truncatedBits := (word2 bitAnd: 16r1FFFFFFF).

	"We must now honour default IEEE rounding mode (round to nearest even)"
	
	"we are below gradual underflow, even if rounded to upper mantissa"
	exponent < -24 ifTrue: [^sign "this can be negative zero"].
	
	"BEWARE: rounding occurs on less than 23bits when gradual underflow"
	exponent <= 0
		ifTrue:
			[mask := 1 bitShift: exponent negated.
			mantissa := mantissa bitOr: 16r800000.
			roundToUpper := (mantissa bitAnd: mask) isZero not
				and: [truncatedBits isZero not
					or: [(mantissa bitAnd: mask - 1) isZero not
						or: [(mantissa bitAnd: mask*2) isZero not]]].
			mantissa := mantissa bitShift: exponent - 1.
			"exponent := exponent + 1"]
		ifFalse:
			[roundToUpper := (truncatedBits bitAnd: 16r10000000) isZero not
				and: [(mantissa bitAnd: 16r1) isZero not
					or: [(truncatedBits bitAnd: 16r0FFFFFFF) isZero not]]
			].
		
	"adjust mantissa and exponent due to IEEE rounding mode"
	roundToUpper
		ifTrue:
			[mantissa := mantissa + 1.
			mantissa > 16r7FFFFF
				ifTrue:
					[mantissa := 0.
					exponent := exponent+1]].

	exponent > 254 ifTrue: ["Overflow"
		exponent := 255.
		self isNaN
			ifTrue: [mantissa isZero
				ifTrue: ["BEWARE: do not convert a NaN to infinity due to truncatedBits"
					mantissa := 1]]
			ifFalse: [mantissa := 0]].
		
	"Encode the word"
	destWord := (sign bitOr: ((exponent max: 0) bitShift: 23)) bitOr: mantissa.
	^ destWord! !

!Float methodsFor: 'converting' stamp: 'jmv 8/30/2010 23:32'!
asTrueFraction
	" Answer a fraction that EXACTLY represents self,
	  a double precision IEEE floating point number.
	  By David N. Smith with significant performance
	  improvements by Luciano Esteban Notarfrancesco.
	  (Version of 11April97).
	Refactoring and simplification by jmv"
	
	^self
		partValues: [ :sign :exponent :mantissa | | zeroBitsCount |
			" Prepare result. If exponent is greater than mantissa size, result is an integer"
			(exponent >= 52 or: [
					zeroBitsCount _ mantissa lowBit - 1.
					exponent + zeroBitsCount >= 52 ])
				ifTrue: [
					"result is an integer number"
					 sign * mantissa bitShift: exponent - 52 ]
				ifFalse:	[
					" This is the 'obvious' way. Better do Luciano's trick below:"
					"result := Fraction
						numerator: sign * mantissa
						denominator: (1 bitShift: 52 - exponent)."
					" Form the result. When exp>52, the exponent is adjusted by
					  the number of trailing zero bits in the mantissa to minimize
					  the (huge) time could be spent in #gcd:. "
					Fraction
						numerator: (sign * (mantissa bitShift: 0 - zeroBitsCount))
						denominator: (1 bitShift: 52 - exponent - zeroBitsCount) ]
		]
		ifInfinite: [ self error: 'Cannot represent infinity as a fraction' ]
		ifNaN: [ self error: 'Cannot represent Not-a-Number as a fraction' ].! !

!Float methodsFor: 'converting' stamp: 'jmv 8/30/2010 23:39'!
partValues: aThreeArgumentBlock
	^ self
		partValues: aThreeArgumentBlock
		ifInfinite: [ self error: 'Can not handle infinity' ]
		ifNaN: [ self error: 'Can not handle Not-a-Number' ].! !

!Float methodsFor: 'converting' stamp: 'jmv 8/30/2010 23:29'!
partValues: aThreeArgumentBlock ifInfinite: aZeroOrOneArgBlock ifNaN: otherZeroOrOneOrTwoArgBlock
	"
	Float pi hex print
	Float pi partValues: [ :sign :exponent :mantissa | { sign hex. exponent hex. mantissa hex} print ]
	0.0 partValues: [ :sign :exponent :mantissa | { sign hex. exponent hex. mantissa hex} print ]
	For 0.0, exponent will be the minimum possible, i.e.  -1023, and mantissa will be 0.
	"
	| shifty sign exponent mantissa expPart fractionPart |

	" Extract the bits of an IEEE double float "
	shifty := ((self basicAt: 1) bitShift: 32) + (self basicAt: 2).

	" Extract the sign and the biased exponent "
	sign := (shifty bitShift: -63) = 0 ifTrue: [1] ifFalse: [-1].
	expPart := (shifty bitShift: -52) bitAnd: 16r7FF.

	" Extract fractional part "
	fractionPart := shifty bitAnd:  16r000FFFFFFFFFFFFF.
	
	" Special cases: infinites and NaN"
	expPart = 16r7FF ifTrue: [
		^fractionPart = 0
			ifTrue: [ aZeroOrOneArgBlock valueWithPossibleArgument: self ]
			ifFalse: [ otherZeroOrOneOrTwoArgBlock valueWithPossibleArgument: self and: fractionPart ]].

	" Replace omitted leading 1 in fraction if appropriate"
	"If expPart = 0, I am +/-zero or a denormal value. In such cases, no implicit leading bit in mantissa"	
	mantissa := expPart = 0
		ifTrue: [ fractionPart ]
		ifFalse: [ fractionPart bitOr: 16r0010000000000000 ].

	" Unbias exponent: 16r3FF is bias"
	exponent := expPart - 16r3FF.

	"Evaluate the block"
	^aThreeArgumentBlock value: sign value: exponent value: mantissa! !

!Float methodsFor: 'private' stamp: 'nice 8/9/2009 21:01'!
absPrintOn: aStream base: base digitCount: digitCount 
	"Print me in the given base, using digitCount significant figures."

	| fuzz x exp q fBase scale logScale xi |
	self isInfinite ifTrue: [^ aStream nextPutAll: 'Inf'].
	fBase := base asFloat.
	"x is myself normalized to [1.0, fBase), exp is my exponent"
	exp := 
		self < 1.0
			ifTrue: [self reciprocalFloorLog: fBase]
			ifFalse: [self floorLog: fBase].
	scale := 1.0.
	logScale := 0.
	[(x := fBase raisedTo: (exp + logScale)) = 0]
		whileTrue:
			[scale := scale * fBase.
			logScale := logScale + 1].
	x := self * scale / x.
	fuzz := fBase raisedTo: 1 - digitCount.
	"round the last digit to be printed"
	x := 0.5 * fuzz + x.
	x >= fBase
		ifTrue: 
			["check if rounding has unnormalized x"
			x := x / fBase.
			exp := exp + 1].
	(exp < 6 and: [exp > -4])
		ifTrue: 
			["decimal notation"
			q := 0.
			exp < 0 ifTrue: [1 to: 1 - exp do: [:i | aStream nextPut: ('0.0000'
at: i)]]]
		ifFalse: 
			["scientific notation"
			q := exp.
			exp := 0].
	[x >= fuzz]
		whileTrue: 
			["use fuzz to track significance"
			xi := x asInteger.
			aStream nextPut: (Character digitValue: xi).
			x := x - xi asFloat * fBase.
			fuzz := fuzz * fBase.
			exp := exp - 1.
			exp = -1 ifTrue: [aStream nextPut: $.]].
	[exp >= -1]
		whileTrue: 
			[aStream nextPut: $0.
			exp := exp - 1.
			exp = -1 ifTrue: [aStream nextPut: $.]].
	q ~= 0
		ifTrue: 
			[aStream nextPut: $e.
			q printOn: aStream]! !


!Fraction methodsFor: 'converting' stamp: 'nice 1/10/2007 02:07'!
asFloat
	"Answer a Float that closely approximates the value of the receiver.
	This implementation will answer the closest floating point number to
	the receiver.
	It uses the IEEE 754 round to nearest even mode
	(can happen in case denominator is a power of two)"
	
	| a b q r exponent floatExponent n ha hb hq q1 |
	a := numerator abs.
	b := denominator abs.
	ha := a highBit.
	hb := b highBit.
	
	"If both numerator and denominator are represented exactly in floating point number,
	then fastest thing to do is to use hardwired float division"
	(ha < 54 and: [hb < 54]) ifTrue: [^numerator asFloat / denominator asFloat].
	
	"Try and obtain a mantissa with 54 bits.
	First guess is rough, we might get one more bit or one less"
	exponent := ha - hb - 54.
	exponent > 0
		ifTrue: [b := b bitShift: exponent]
		ifFalse: [a := a bitShift: exponent negated].
	q := a quo: b.
	r := a - (q * b).
	hq := q highBit.
	
	"check for gradual underflow, in which case we should use less bits"
	floatExponent := exponent + hq - 1.
	n := floatExponent > -1023
		ifTrue: [54]
		ifFalse: [54 + floatExponent + 1022].
	
	hq > n
		ifTrue: [exponent := exponent + hq - n.
			r := (q bitAnd: (1 bitShift: hq - n) - 1) * b + r.
			q := q bitShift: n - hq].
	hq < n
		ifTrue: [exponent := exponent + hq - n.
			q1 := (r bitShift: n - hq) quo: b.
			q := (q bitShift: n - hq) bitAnd: q1.
			r := (r bitShift: n - hq) - (q1 * b)].
		
	"check if we should round upward.
	The case of exact half (q bitAnd: 1) isZero not & (r isZero)
	will be handled by Integer>>asFloat"
	((q bitAnd: 1) isZero or: [r isZero])
		ifFalse: [q := q + 1].
		
	^ (self positive
		ifTrue: [q asFloat]
		ifFalse: [q asFloat negated])
		timesTwoPower: exponent! !


!Integer methodsFor: 'converting' stamp: 'nice 7/8/2006 01:05'!
asFloat
	"Answer a Float that represents the value of the receiver.
	Optimized to process only the significant digits of a LargeInteger.
	SqR: 11/30/1998 21:1
	
	This algorithm does honour IEEE 754 round to nearest even mode.
	Numbers are first rounded on nearest integer on 53 bits.
	In case of exact half difference between two consecutive integers (2r0.1),
	there are two possible choices (two integers are as near, 0 and 1)
	In this case, the nearest even integer is chosen.
	examples (with less than 53bits for clarity)
	2r0.00001 is rounded to 2r0
	2r1.00001 is rounded to 2.1
	2r0.1 is rounded to 2r0 (nearest event)
	2r1.1 is rounded to 2.10 (neraest even)
	2r0.10001 is rounded to 2r1
	2r1.10001 is rounded to 2.10"
	
	| abs shift sum delta mask trailingBits carry |
	self isZero
		ifTrue: [^ 0.0].
	abs := self abs.

	"Assume Float is a double precision IEEE 754 number with 53bits mantissa.
	We should better use some Float class message for that (Float precision)..."
	delta := abs highBit - 53.
	delta > 0
		ifTrue: [mask := (1 bitShift: delta) - 1.
			trailingBits := abs bitAnd: mask.
			"inexact := trailingBits isZero not."
			carry := trailingBits bitShift: 1 - delta.
			abs := abs bitShift: delta negated.
			shift := delta.
			(carry isZero
					or: [(trailingBits bitAnd: (mask bitShift: -1)) isZero
							and: [abs even]])
				ifFalse: [abs := abs + 1]]
		ifFalse: [shift := 0].
	
	"now, abs has no more than 53 bits, we can do exact floating point arithmetic"
	sum := 0.0.
	1 to: abs size do:
		[:byteIndex | 
		sum := ((abs digitAt: byteIndex) asFloat timesTwoPower: shift) + sum.
		shift := shift + 8].
	^ self positive
			ifTrue: [sum]
			ifFalse: [sum negated]! !

!Integer methodsFor: 'private' stamp: 'nice 1/26/2008 02:12'!
digitLogic: arg op: op length: len 
	| i result neg1 neg2 rneg z1 z2 rz b1 b2 b |
	neg1 := self negative.
	neg2 := arg negative.
	rneg := ((neg1
				ifTrue: [-1]
				ifFalse: [0])
				perform: op
				with: (neg2
						ifTrue: [-1]
						ifFalse: [0]))
				< 0.
	result := Integer new: len neg: rneg.
	rz := z1 := z2 := true.
	i := 0.
	[(i := i + 1) <= len
		or: ["mind a carry on result that might go past len digits"
			rneg and: [rz
				and: [result := result growby: 1.
					true]]]]
		whileTrue: [b1 := self digitAt: i.
			neg1
				ifTrue: [b1 := z1
								ifTrue: [b1 = 0
										ifTrue: [0]
										ifFalse: [z1 := false.
											256 - b1]]
								ifFalse: [255 - b1]].
			b2 := arg digitAt: i.
			neg2
				ifTrue: [b2 := z2
								ifTrue: [b2 = 0
										ifTrue: [0]
										ifFalse: [z2 := false.
											256 - b2]]
								ifFalse: [255 - b2]].
			b := b1 perform: op with: b2.
			result
				digitAt: i
				put: (rneg
						ifTrue: [rz
								ifTrue: [b = 0
										ifTrue: [0]
										ifFalse: [rz := false.
											256 - b]]
								ifFalse: [255 - b]]
						ifFalse: [b])].
	^ result normalize! !


!Number class methodsFor: 'constants' stamp: 'GabrielOmarCotelli 5/23/2009 20:46'!
one

	^1! !


!Float class methodsFor: 'instance creation' stamp: 'nice 5/30/2006 03:13'!
fromIEEE32Bit: word
	"Convert the given 32 bit word (which is supposed to be a positive 32bit value) from a 32bit IEEE floating point representation into an actual Squeak float object (being 64bit wide). Should only be used for conversion in FloatArrays or likewise objects."
	
	| sign mantissa exponent newFloat delta |
	word negative ifTrue: [^ self error:'Cannot deal with negative numbers'].
	word = 0 ifTrue: [^ 0.0].
	sign := word bitAnd: 16r80000000.
	word = sign ifTrue: [^self negativeZero].
	
	exponent := ((word bitShift: -23) bitAnd: 16rFF) - 127.
	mantissa := word bitAnd:  16r7FFFFF.

	exponent = 128 ifTrue:["Either NAN or INF"
		mantissa = 0 ifFalse:[^ Float nan].
		sign = 0 
			ifTrue:[^ Float infinity]
			ifFalse:[^ Float infinity negated]].

	exponent = -127 ifTrue: [
		"gradual underflow (denormalized number)
		Remove first bit of mantissa and adjust exponent"
		delta := mantissa highBit.
		mantissa := (mantissa bitShift: 1) bitAnd: (1 bitShift: delta) - 1.
		exponent := exponent + delta - 23].
	
	"Create new float"
	newFloat := self new: 2.
	newFloat basicAt: 1 put: ((sign bitOr: (1023 + exponent bitShift: 20)) bitOr: (mantissa bitShift: -3)).
	newFloat basicAt: 2 put: ((mantissa bitAnd: 7) bitShift: 29).
	^newFloat! !

!Float class methodsFor: 'constants' stamp: 'GabrielOmarCotelli 5/25/2009 15:42'!
one
	
	^1.0! !


!FloatArray methodsFor: 'arithmetic' stamp: 'nice 11/24/2007 00:10'!
adaptToNumber: rcvr andSend: selector
	"If I am involved in arithmetic with a Number. If possible,
	convert it to a float and perform the (more efficient) primitive operation."
	selector == #+ ifTrue:[^self + rcvr].
	selector == #* ifTrue:[^self * rcvr].
	selector == #- ifTrue:[^self negated += rcvr].
	selector == #/ ifTrue:[
		"DO NOT USE TRIVIAL CODE
			^self reciprocal * rcvr
		BECAUSE OF GRADUAL UNDERFLOW
		self should: (1.0e-39 / (FloatArray with: 1.0e-39)) first < 2."
			^(self class new: self size withAll: rcvr) / self
		].
	^super adaptToNumber: rcvr andSend: selector! !


!SmallInteger methodsFor: 'bit manipulation' stamp: 'nice 7/8/2008 02:39'!
highBitOfMagnitude
	"Answer the index of the high order bit of the receiver, or zero if the  
	receiver is zero. This method is used for negative SmallIntegers as well,  
	since Squeak's LargeIntegers are sign/magnitude."
	
	self < 0 ifTrue: [
		"Beware: do not use highBitOfPositiveReceiver
		because self negated is not necessarily a SmallInteger
		(see SmallInteger minVal)"
		^self negated highBitOfMagnitude].
	
	"Implementation note: this method could be as well inlined here."
	^self highBitOfPositiveReceiver! !

Float removeSelector: #asApproximateFraction!
Float removeSelector: #isInf!
Float removeSelector: #raisedTo:!
Float removeSelector: #shallowCopy!
Number removeSelector: #isInf!
