QotD    "To be or not to be" – Shakespeare

benchFib  "Handy send-heavy benchmark"
	"(result // seconds to run) = approx calls per second"
	" | r t |
	  t := Time millisecondsToRun: [r := 26 benchFib].
	  (r * 1000) // t"
	"138000 on a Mac 8100/100"
	^ self < 2
		ifTrue: [1] 
		ifFalse: [(self-1) benchFib + (self-2) benchFib + 1]

fibonacciSequence
	"Yields an infinite sequence of fibonacci numbers."
	
	^ Generator on: [ :generator |
		| a b |
		a := 0. b := 1.
		[ a := b + (b := a).
		  generator yield: a ]
			repeat ]

 squeak my.image bigFib.cs

'From Squeak5.2 of 13 December 2018 [latest update: #18225] on 17 April 2019 at 2:12:53 pm'!

!Integer methodsFor: 'benchmarks' stamp: 'tpr 4/17/2019 13:21'!
fastDoublingFib
"derived from https://www.nayuki.io/page/fast-fibonacci-algorithms"
"(1 to: 20) collect:[:i| i fastDoublingFib]"
"testing a quite large one - "
"8577 fastDoublingFib= 13724780954457889052017147206983806244049002655849289934662148526555403650297300643609932653364630032094175733360509578504423049114004867523161091564824674600978308740378989479162611031273424686573759784740689516901416473328655422390895971263265867635819510949686102571388751758998017349379498918930220900180038324615253426530740473855269056304380498630108126734047701362218655030270360608989081398866489746698916626888016696892691933237089180504631788915650757757515944644732966345269202761471025651071790297611079540881155092137592980230998234868586211881093892491570520577408961869977892273540956424095750855208529072246641778982103984467921868950012668004047986803017482248992968482737462668300684879633714025755790485860328854796518843956263863014632532331145699658530054942590047273403691531821918862996422405159427262092477196755988981309029424760342802374213122162727840557722145891090413688461745240415668189577836068480363407847582529735341950500636735281963089675493707159434777756081146452522323681782226760627277553296721358921412115264845467834979154061137421532609247762981818564578019888974692581079593575783553856910367568474613323528337733872069223030834774749130478360574004172522316484339530942110067893000847800932306298725285623628731149337468217751734165148732164148285915275115006479682658150442259002271790547596033006363411193653337536041106069912826015502035140618407668385378737477702597473151509972754111640855347958033314453349633268543893894677097708945041254623018915871109789412793709229204261914803477697183287924195770678873001065036313926288444791424871512110658175954743584548831946767673488152740675550518235698898217693311515366329280005757014637854214769152690638778904780724293185353992279724740604674926819294787586671833537117545443846365508358918882"
	| a b c |
	a :=  0.
	b := 1.
	self highBit to: 1 by: -1 do:[:i||d e|
		d := (b bitShift: 1) - a karatsubaTimes: a.
		e := a squared + b squared.
		a := d.
		b := e.
		(self bitAt: i) = 1  ifTrue:[
			c := a + b.
			a := b.
			b := c]
		].
	^a! !

!Integer methodsFor: '*FactorialContest' stamp: 'nice 7/23/2008 22:41'!
copyDigitsFrom: start to: stop
	"make a new integer keeping only some digits of self.
	Implementation should be faster than:
	self - ((self bitShift: -8*stop) bitShift: 8*stop) bitShift: -8*start"
	
	| len slice |
	start > 0 ifFalse: [^self error: 'start index should be at least 1'].
	len := self digitLength.
	(start > len or: [start > stop]) ifTrue: [^0].
	stop >= len
		ifTrue: [start = 1 ifTrue: [^self].
				len := len - start + 1]
		ifFalse: [len := stop - start + 1].
	slice := self class new: len neg: self negative.
	slice replaceFrom: 1 to: len with: self startingAt: start.
	^slice normalize! !

!Integer methodsFor: '*FactorialContest' stamp: 'nice 7/23/2008 22:00'!
isLargeEnoughForKaratsuba
	"Answer whether it is interesting to engage a Karatsuba multiplication algorithm with this integer.
	Karatsuba is more efficient than naive algortithm assymptotically for large numbers,
	but it has an overhead justifying to not use it on Small integers"
	
	^false! !

!Integer methodsFor: '*FactorialContest' stamp: 'nice 7/22/2008 21:01'!
karatsubaTimes: anInteger
	"eventually use karatsuba algorithm to perform the multiplication
	default to regular multiplication"
	
	^self * anInteger! !

!Integer methodsFor: '*FactorialContest' stamp: 'nice 7/22/2008 21:33'!
lowestNDigits: N
	"make a new integer keeping only N least significant digits of self.
	Implementation is faster than self - ((self bitShift: -8*N) bitShift: 8*N)"
	
	| low |
	N >= self digitLength ifTrue: [^self].
	low := self class new: N neg: self negative.
	low replaceFrom: 1 to: N with: self startingAt: 1.
	^low normalize! !


!LargePositiveInteger methodsFor: '*FactorialContest' stamp: 'nice 6/13/2011 13:48'!
isLargeEnoughForKaratsuba
	"See super.
	One possibility would be to count non zero, because regular
	multiplication might be more efficient with a sparse integer.
		^self hasMoreNonZeroDigitsThan: 256
	However, this would be called recursively, and would result
	in a cost O(digitLength^2),
	Another idea would be to make the sparsity test random.
	Here, we keep the simplest thing"
	
	^self digitLength >= 160! !

!LargePositiveInteger methodsFor: '*FactorialContest' stamp: 'nice 9/16/2008 21:42'!
karatsubaTimes: anInteger
	"eventually use Karatsuba algorithm to perform the multiplication"
	
	| half xHigh xLow yHigh yLow low high mid xLen yLen |
	(anInteger isLargeEnoughForKaratsuba
		and: [self isLargeEnoughForKaratsuba]) ifFalse: [^self timesInteger: anInteger].
	
	"Check if length ratio is more than 2, and engage a loop
	to operate on integers with well equilibrated lengths.
	Note that we only add overhead at this level,
	but we hope to gain in lower level recursion"
	(xLen := self digitLength) >= (yLen := anInteger digitLength)
		ifTrue: [(half := xLen bitShift: -1) >= yLen
			ifTrue: [^(0 to: xLen by: yLen) detectSum: [:yShift |
				((self copyDigitsFrom: yShift + 1 to: yShift + yLen)
					karatsubaTimes: anInteger)
						bitShift: 8 * yShift]]]
		ifFalse: [(half := yLen bitShift: -1) >= xLen
			ifTrue: [^(0 to: yLen by: xLen) detectSum: [:xShift |
				(self karatsubaTimes:
					(anInteger copyDigitsFrom: xShift + 1 to: xShift + xLen))
						bitShift: 8 * xShift]]].
	
	"At this point, lengths are well equilibrated, Divide each integer in two halves"
	xHigh := self bitShift: -8 * half.
	xLow := self lowestNDigits: half.
	yHigh := anInteger bitShift: -8 * half.
	yLow := anInteger lowestNDigits: half.
	
	"Karatsuba trick: perform with 3 multiplications instead of 4"
	low := xLow karatsubaTimes: yLow.
	high := xHigh karatsubaTimes: yHigh.
	mid := (xHigh + xLow karatsubaTimes: yHigh + yLow) - (low + high).
	
	"Sum the parts of decomposition"
	^low + (mid bitShift: 8*half) + (high bitShift: 16*half)
	
"
| a b |
a := 1000 factorial.
b := 2000 factorial.
{a digitLength. b digitLength}.
self assert: (a karatsubaTimes: b) - (a * b) = 0.
[Smalltalk garbageCollect.
[10 timesRepeat: [a karatsubaTimes: b]] timeToRun] value /
[Smalltalk garbageCollect.
[10 timesRepeat: [a * b]] timeToRun] value asFloat
"! !

!LargePositiveInteger methodsFor: '*FactorialContest' stamp: 'nice 6/13/2011 13:49'!
squared
	"eventually use a divide and conquer algorithm to perform the multiplication"
	
	| half xHigh xLow low high mid |
	(self digitLength >= 160) ifFalse: [^self * self].
	
	"Divide digits in two halves"
	half := self digitLength bitShift: -1.
	xHigh := self bitShift: -8 * half.
	xLow := self lowestNDigits: half.
	
	"eventually use karatsuba"
	low := xLow squared.
	high := xHigh squared.
	mid := xLow karatsubaTimes: xHigh.
	
	"Sum the parts of decomposition"
	^low + (mid bitShift: 8*half+1) + (high bitShift: 16*half)
	
"
| a |
a := 8000 factorial.
a digitLength.
self assert: a * a - a squared = 0.
[Smalltalk garbageCollect.
[2 timesRepeat: [a squared]] timeToRun] value /
[Smalltalk garbageCollect.
[2 timesRepeat: [a * a]] timeToRun] value asFloat
"! !

!LargePositiveInteger methodsFor: '*FactorialContest' stamp: 'nice 9/16/2008 21:41'!
timesInteger: anInteger 
	"Primitive. Multiply the receiver by the argument and answer with an
	Integer result. Fail if either the argument or the result is not a
	SmallInteger or a LargePositiveInteger less than 2-to-the-30th (1073741824). Optional. See
	Object documentation whatIsAPrimitive. "

	<primitive: 29>
	^ self digitMultiply: anInteger 
					neg: self negative ~~ anInteger negative! !

"Postscript:
Load a big fibonnaci calulator and then run 4784969 fastDoublingFib, write the time taken to stdout  and quit"
FileStream  stdout nextPutAll: [4784969 fastDoublingFib] timeToRun printString; flush.
Smalltalk snapshot: false andQuit: true
!