 ## [squeak-dev] The Inbox: Kernel-nice.1217.mcz

commits at source.squeak.org commits at source.squeak.org
Fri Apr 26 16:02:06 UTC 2019

```Nicolas Cellier uploaded a new version of Kernel to project The Inbox:
http://source.squeak.org/inbox/Kernel-nice.1217.mcz

==================== Summary ====================

Name: Kernel-nice.1217
Author: nice
Time: 26 April 2019, 6:01:39.036232 pm
UUID: 39478bc1-0f55-4747-8332-45a6d225b281
Ancestors: Kernel-eem.1215

Publish Karatsuba and 3-way Toom-Cooke algorithms for fast large integer multiplication.
These are efficient divide and conquer algorithms that recursively split a problem into simpler problems.

Currently this fast multiplication is used only in LargePositiveInteger>>squared.
Slightly modify raisedToInteger: to use that squared.
This result on a small penalty for small number exponentation, but large speed up for big integers.

Also publish another Karatsuba-like square root algorithm that should accelerate sqrt of huge integers.

Examples on a 64 bits image/VM

[10 raisedToInteger: 30] bench.
'3,770,000 per second. 265 nanoseconds per run.' - before
'3,500,000 per second. 286 nanoseconds per run.' - after
[10.0 raisedToInteger: 30] bench.
'10,900,000 per second. 91.6 nanoseconds per run.' - before
'9,540,000 per second. 105 nanoseconds per run.' - after

| x |
x := 100 factorial.
[x raisedToInteger: 30] bench.
'15,900 per second. 63 microseconds per run.' - before
'19,500 per second. 51.3 microseconds per run.' - after

| x |
x := 1000 factorial.
[x raisedToInteger: 30] bench.
'58.8 per second. 17 milliseconds per run.' - before
'179 per second. 5.6 milliseconds per run.' - after

| x |
x := 500 factorial - 100 factorial.
{
[x sqrtFloor] bench.
[x sqrtFloorNewtonRaphson] bench.
}
#(
'42,300 per second. 23.6 microseconds per run.'
'8,240 per second. 121 microseconds per run.')

=============== Diff against Kernel-eem.1215 ===============

+ ----- Method: Integer>>fastMultiply: (in category 'arithmetic') -----
+ fastMultiply: anInteger
+ 	^self * anInteger!

+ ----- Method: Integer>>karatsubaTimes: (in category 'arithmetic') -----
+ karatsubaTimes: anInteger
+ 	"Eventually use Karatsuba algorithm to perform the multiplication.
+ 	This is efficient only for large integers (see subclass).
+ 	By default, use regular multiplication."
+
+ 	^self * anInteger!

Item was changed:
----- Method: Integer>>sqrtFloor (in category 'mathematical functions') -----
sqrtFloor
+ 	"Return the integer part of the square root of self
+ 	Assume self >= 0
+ 	The following invariants apply:
+ 	1) self sqrtFloor squared <= self
+ 	2) (self sqrtFloor + 1) squared > self"
- 	"Return the integer part of the square root of self"

+ 	^self sqrtFloorNewtonRaphson!
- 	| guess delta |
- 	guess := 1 bitShift: self highBit + 1 // 2.
- 	[
- 		delta := guess squared - self // (guess bitShift: 1).
- 		delta = 0 ] whileFalse: [
- 			guess := guess - delta ].
- 	^guess - 1!

+ ----- Method: Integer>>sqrtFloorNewtonRaphson (in category 'mathematical functions') -----
+ sqrtFloorNewtonRaphson
+ 	"Return the integer part of the square root of self.
+ 	Use a Newton-Raphson iteration since it converges quadratically
+ 	(twice more bits of precision at each iteration)"
+
+ 	| guess delta |
+ 	guess := 1 bitShift: self highBit + 1 // 2.
+ 	[
+ 		delta := guess squared - self // (guess bitShift: 1).
+ 		delta = 0 ] whileFalse: [
+ 			guess := guess - delta ].
+ 	^guess - 1!

+ ----- Method: Integer>>sqrtRem (in category 'mathematical functions') -----
+ sqrtRem
+ 	"Return an array with floor sqrt and sqrt remainder.
+ 	Assume self >= 0.
+ 	The following invariants apply:
+ 	1) self  sqrtRem first squared <= self
+ 	2) (self sqrtRem first + 1) squared > self
+ 	3) self sqrtRem first squared + self sqrtRem last = self"
+
+ 	| s |
+ 	s := self sqrtFloorNewtonRaphson.
+ 	^{ s. self - s squared } !

+ ----- Method: Integer>>toom3Times: (in category 'arithmetic') -----
+ toom3Times: anInteger
+ 	^self * anInteger!

+ ----- Method: LargePositiveInteger>>butLowestNDigits: (in category 'private') -----
+ butLowestNDigits: N
+ 	"make a new integer removing N least significant digits of self."
+
+ 	^self bitShiftMagnitude: -8 * N!

+ ----- Method: LargePositiveInteger>>copyDigitsFrom:to: (in category 'private') -----
+ copyDigitsFrom: start to: stop
+ 	"Make a new integer keeping only some digits of self."
+
+ 	| 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.
+ 	slice replaceFrom: 1 to: len with: self startingAt: start.
+ 	^slice normalize!

+ ----- Method: LargePositiveInteger>>fastMultiply: (in category 'arithmetic') -----
+ fastMultiply: anInteger
+ 	"Eventually use Karatsuba or 3-way Toom-Cook algorithm to perform the multiplication"
+
+ 	| xLen yLen |
+ 	"arrange to have the receiver be the shorter"
+ 	(xLen := self digitLength) > (yLen := anInteger digitLength)
+ 		ifTrue: [^ anInteger fastMultiply: self ].
+
+ 	"If too short to be worth, fallback to naive algorithm"
+ 	(xLen >= 600) ifFalse: [^self * anInteger].
+ 	(xLen >= 6000) ifFalse: [^self karatsubaTimes: anInteger].
+ 	^self toom3Times: anInteger!

+ ----- Method: LargePositiveInteger>>karatsubaTimes: (in category 'arithmetic') -----
+ karatsubaTimes: anInteger
+ 	"Eventually use Karatsuba algorithm to perform the multiplication
+ 	The Karatsuba algorithm divide number in two smaller parts.
+ 	Then it uses tricks to perform only 3 multiplications.
+ 	See https://en.wikipedia.org/wiki/Karatsuba_algorithm"
+
+ 	| half xHigh xLow yHigh yLow low high mid xLen yLen |
+ 	"arrange to have the receiver be the shorter"
+ 	(xLen := self digitLength) > (yLen := anInteger digitLength)
+ 		ifTrue: [^ anInteger karatsubaTimes: self ].
+
+ 	"If too short to be worth, fallback to naive algorithm"
+ 	(xLen >= 600) ifFalse: [^self * anInteger].
+
+ 	"Seek for integers of about the same length, else split the long one"
+ 	yLen >= (7 * xLen bitShift: -2)
+ 		ifTrue:
+ 			[^(0 to: yLen - 1 by: xLen + 1)  digitShiftSum: [:yShift |
+ 				self karatsubaTimes:
+ 					(anInteger copyDigitsFrom: yShift + 1 to: yShift + 1 + xLen)]].
+
+ 	"At this point, lengths are similar, divide each integer in two halves"
+ 	half := (yLen + 1 bitShift: -1)  bitClear: 2r11.
+ 	xLow := self lowestNDigits: half.
+ 	xHigh := self butLowestNDigits: half.
+ 	yLow := anInteger lowestNDigits: half.
+ 	yHigh := anInteger butLowestNDigits: half.
+
+ 	"Karatsuba trick: perform with 3 multiplications instead of 4"
+ 	low := xLow karatsubaTimes: yLow.
+ 	high := xHigh karatsubaTimes: yHigh.
+ 	mid := high + low + (xHigh - xLow karatsubaTimes: yLow - yHigh).
+
+ 	"Sum the parts of decomposition"
+ 	^low + (mid bitShiftMagnitude: 8*half) + (high bitShiftMagnitude: 16*half)
+
+ "
+ | a b |
+ a := 650 factorial-1.
+ b := 700 factorial-1.
+ {a digitLength. b digitLength}.
+ self assert: (a karatsubaTimes: b) - (a * b) = 0.
+ [Smalltalk garbageCollect.
+ [1000 timesRepeat: [a karatsubaTimes: b]] timeToRun] value /
+ [Smalltalk garbageCollect.
+ [1000 timesRepeat: [a * b]] timeToRun] value asFloat
+ "!

+ ----- Method: LargePositiveInteger>>lowestNDigits: (in category 'private') -----
+ lowestNDigits: N
+ 	"Make a new integer keeping only N least significant digits of self"
+
+ 	| low |
+ 	N >= self digitLength ifTrue: [^self].
+ 	low := self class new: N.
+ 	low replaceFrom: 1 to: N with: self startingAt: 1.
+ 	^low normalize!

+ ----- Method: LargePositiveInteger>>sqrtFloor (in category 'mathematical functions') -----
+ sqrtFloor
+ 	"Like super, but use a faster divide and conquer strategy if it's worth"
+
+ 	self digitLength >= 64 ifFalse: [^self sqrtFloorNewtonRaphson].
+ 	^self sqrtRem first!

+ ----- Method: LargePositiveInteger>>sqrtRem (in category 'mathematical functions') -----
+ sqrtRem
+ 	"Like super, but use a divide and conquer method to perform this operation.
+ 	See Paul Zimmermann. Karatsuba Square Root. [Research Report] RR-3805, INRIA. 1999, pp.8. <inria-00072854>
+ 	https://hal.inria.fr/inria-00072854/PDF/RR-3805.pdf"
+
+ 	| n  qr s r sr high mid low |
+ 	n := self digitLength bitShift: -2.
+ 	n >= 16 ifFalse: [^super sqrtRem].
+ 	high := self butLowestNDigits: n * 2.
+ 	mid := self copyDigitsFrom: n + 1 to: n * 2.
+ 	low := self lowestNDigits: n.
+
+ 	sr := high sqrtRem.
+ 	qr := (sr last bitShift: 8 * n) + mid digitDiv: (sr first bitShift: 1) neg: false.
+ 	s := (sr first bitShift: 8 * n) + qr first.
+ 	r := (qr last bitShift: 8 * n) + low - qr first squared.
+ 	r negative
+ 		ifTrue:
+ 			[r := (s bitShift: 1) + r - 1.
+ 			s := s - 1].
+ 	sr at: 1 put: s; at: 2 put: r.
+ 	^sr
+ 	!

+ ----- Method: LargePositiveInteger>>squared (in category 'mathematical functions') -----
+ squared
+ 	"Eventually use a divide and conquer algorithm to perform the multiplication"
+
+ 	(self digitLength >= 400) ifFalse: [^self * self].
+ 	(self digitLength >= 1600) ifFalse: [^self squaredKaratsuba].
+ 	^self squaredToom3!

+ ----- Method: LargePositiveInteger>>squaredKaratsuba (in category 'mathematical functions') -----
+ squaredKaratsuba
+ 	"Use a divide and conquer algorithm to perform the multiplication.
+ 	Split in two parts like Karatsuba, but economize 2 additions by using asymetrical product."
+
+ 	| half xHigh xLow low high mid |
+
+ 	"Divide digits in two halves"
+ 	half := self digitLength + 1 // 2 bitClear: 2r11.
+ 	xLow := self lowestNDigits: half.
+ 	xHigh := self butLowestNDigits: 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 := 440 factorial-1.
+ a digitLength.
+ self assert: a * a - a squaredKaratsuba = 0.
+ [Smalltalk garbageCollect.
+ [2000 timesRepeat: [a squaredKaratsuba]] timeToRun] value /
+ [Smalltalk garbageCollect.
+ [2000 timesRepeat: [a * a]] timeToRun] value asFloat
+ "!

+ ----- Method: LargePositiveInteger>>squaredToom3 (in category 'mathematical functions') -----
+ squaredToom3
+ 	"Use a 3-way Toom-Cook divide and conquer algorithm to perform the multiplication"
+
+ 	| third x0 x1 x2 x20 z0 z1 z2 z3 z4 |
+ 	"divide in 3 parts"
+ 	third := self digitLength + 2 // 3 bitClear: 2r11.
+ 	x2 := self butLowestNDigits: third * 2.
+ 	x1 := self copyDigitsFrom: third + 1 to: third * 2.
+ 	x0 := self lowestNDigits: third.
+
+ 	"Toom-3 trick: 5 multiplications instead of 9"
+ 	z0 := x0 squared.
+ 	z4 := x2 squared.
+ 	x20 := x2 + x0.
+ 	z1 := (x20 + x1) squared.
+ 	x20 := x20 - x1.
+ 	z2 := x20 squared.
+ 	z3 := ((x20 + x2 bitShift: 1) - x0) squared.
+
+ 	"Sum the parts of decomposition"
+ 	z3 := z3 - z1 quo: 3.
+ 	z1 := z1 - z2 bitShift: -1.
+ 	z2 := z2 - z0.
+
+ 	z3 := (z2 - z3 bitShift: -1) + (z4 bitShift: 1).
+ 	z2 := z2 + z1 - z4.
+ 	z1 := z1 - z3.
+ 	^z0 + (z1 bitShift: 8*third) + (z2 bitShift: 16*third) + (z3 + (z4 bitShift: 8*third) bitShift: 24*third)
+
+ "
+ | a |
+ a := 1400 factorial-1.
+ a digitLength.
+ self assert: a * a - a squaredToom3 = 0.
+ [Smalltalk garbageCollect.
+ [1000 timesRepeat: [a squaredToom3]] timeToRun] value /
+ [Smalltalk garbageCollect.
+ [1000 timesRepeat: [a squaredKaratsuba]] timeToRun] value asFloat
+ "!

+ ----- Method: LargePositiveInteger>>toom3Times: (in category 'arithmetic') -----
+ toom3Times: anInteger
+ 	"Eventually use a variant of Toom-Cooke for performing multiplication.
+ 	Toom-Cooke is a generalization of Karatsuba divide and conquer algorithm.
+ 	It divides operands in n parts instead of 2.
+ 	See https://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication
+ 	Here we divide each operand in 3 parts, thus the name Toom-3."
+
+ 	| third x2 x1 x0 y2 y1 y0 xLen yLen y20 z4 z3 z2 z1 z0 x20 |
+ 	"arrange to have the receiver be the shorter"
+ 	(xLen := self digitLength) > (yLen := anInteger digitLength)
+ 		ifTrue: [^anInteger toom3Times: self ].
+
+ 	"If too short to be worth, fallback to Karatsuba algorithm"
+ 	(xLen > 6000) ifFalse: [^self karatsubaTimes: anInteger].
+
+ 	"Seek for well balanced integers"
+ 	yLen > (5 * xLen bitShift: -2)
+ 		ifTrue: [^(0 to: yLen - 1 by: xLen + 1)  digitShiftSum: [:yShift |
+ 				self toom3Times: (anInteger copyDigitsFrom: yShift + 1 to: yShift +1 + xLen)]].
+
+ 	"At this point, lengths are well balanced, divide in 3 parts"
+ 	third := yLen + 2 // 3 bitClear: 2r11.
+ 	x2 := self butLowestNDigits: third * 2.
+ 	x1 := self copyDigitsFrom: third + 1 to: third * 2.
+ 	x0 := self lowestNDigits: third.
+ 	y2 := anInteger butLowestNDigits: third * 2.
+ 	y1 := anInteger copyDigitsFrom: third + 1 to: third * 2.
+ 	y0 := anInteger lowestNDigits: third.
+
+ 	"Toom-3 trick: 5 multiplications instead of 9"
+ 	z0 := x0 toom3Times: y0.
+ 	z4 := x2 toom3Times: y2.
+ 	x20 := x2 + x0.
+ 	y20 := y2 + y0.
+ 	z1 := x20 + x1 toom3Times: y20 + y1.
+ 	x20 := x20 - x1.
+ 	y20 := y20 - y1.
+ 	z2 := x20 toom3Times: y20.
+ 	z3 := (x20 + x2 bitShift: 1) - x0 toom3Times: (y20 + y2 bitShift: 1) - y0.
+
+ 	"Sum the parts of decomposition"
+ 	z3 := z3 - z1 quo: 3.
+ 	z1 := z1 - z2 bitShift: -1.
+ 	z2 := z2 - z0.
+
+ 	z3 := (z2 - z3 bitShift: -1) + (z4 bitShift: 1).
+ 	z2 := z2 + z1 - z4.
+ 	z1 := z1 - z3.
+ 	^z0 + (z1 bitShift: 8*third) + (z2 bitShift: 16*third) + (z3 + (z4 bitShift: 8*third) bitShift: 24*third)
+
+ "
+ | a b |
+ a :=5000 factorial - 1.
+ b := 6000 factorial - 1.
+ a digitLength min: b digitLength.
+ a digitLength max: b digitLength.
+ (a toom3Times: b) = (a * b) ifFalse: [self error].
+ [Smalltalk garbageCollect.
+ [300 timesRepeat: [a toom3Times: b]] timeToRun] value /
+ [Smalltalk garbageCollect.
+ [300 timesRepeat: [a karatsubaTimes: b]] timeToRun] value asFloat
+ "!

Item was changed:
----- Method: Number>>raisedToInteger: (in category 'mathematical functions') -----
raisedToInteger: anInteger
"The 0 raisedToInteger: 0 is an special case. In some contexts must be 1 and in others must
be handled as an indeterminate form.
I take the first context because that's the way that was previously handled.
Maybe further discussion is required on this topic."

| bitProbe result |
anInteger negative ifTrue: [^(self raisedToInteger: anInteger negated) reciprocal].
bitProbe := 1 bitShift: anInteger highBit - 1.
result := self class one.
[
(anInteger bitAnd: bitProbe) > 0 ifTrue: [ result := result * self ].
(bitProbe := bitProbe bitShift: -1) > 0 ]
+ 		whileTrue: [ result := result squared ].
- 		whileTrue: [ result := result * result ].
^result!

```