<div dir="ltr"><div>Err, I messed up with the quo/rem signs...</div><div>Tests pass, but there are not enough tests!<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le sam. 27 avr. 2019 à 10:41, <<a href="mailto:commits@source.squeak.org">commits@source.squeak.org</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Nicolas Cellier uploaded a new version of Kernel to project The Inbox:<br>
<a href="http://source.squeak.org/inbox/Kernel-nice.1218.mcz" rel="noreferrer" target="_blank">http://source.squeak.org/inbox/Kernel-nice.1218.mcz</a><br>
<br>
==================== Summary ====================<br>
<br>
Name: Kernel-nice.1218<br>
Author: nice<br>
Time: 27 April 2019, 10:41:24.794539 am<br>
UUID: 21f74fbe-a0cd-4b6f-86e9-be13465d57fe<br>
Ancestors: Kernel-nice.1217<br>
<br>
Implement the recursive fast division of Burnikel-Ziegler for large integers and connect it to digitDiv:neg: when operands are large enough.<br>
<br>
This is not the fastest known division which is a composition of Barrett and Newton-Raphson inversion - but is easy to implement and should have similar performances for at least a few thousand bytes long integers - see for example <a href="http://bioinfo.ict.ac.cn/~dbu/AlgorithmCourses/Lectures/Lec5-Fast-Division-Hasselstrom2003.pdf" rel="noreferrer" target="_blank">http://bioinfo.ict.ac.cn/~dbu/AlgorithmCourses/Lectures/Lec5-Fast-Division-Hasselstrom2003.pdf</a> <br>
<br>
Use digitDiv:neg: in large integer printString so as to obtain the quotient (head) and remainder (tail) in a single operation. Together with divide and conquer division, this results in a factor of about 3x for 50000 factorial printString.<br>
<br>
Implement the 4-way Toom-Cook squaring variant of Chung-Hasan. This over-performs the symetrical squaredToom3 even for medium size (800 bytes).<br>
<br>
=============== Diff against Kernel-nice.1217 ===============<br>
<br>
Item was added:<br>
+ ----- Method: Integer>>digitDiv21: (in category 'private') -----<br>
+ digitDiv21: anInteger<br>
+ <br>
+ ^(self digitDiv: anInteger neg: false) collect: #normalize!<br>
<br>
Item was added:<br>
+ ----- Method: Integer>>digitDiv32: (in category 'private') -----<br>
+ digitDiv32: anInteger<br>
+ <br>
+ ^(self digitDiv: anInteger neg: false) collect: #normalize!<br>
<br>
Item was changed:<br>
----- Method: Integer>>digitDiv:neg: (in category 'private') -----<br>
+ digitDiv: anInteger neg: aBoolean <br>
+ ^self primDigitDiv: anInteger neg: aBoolean!<br>
- digitDiv: arg neg: ng <br>
- "Answer with an array of (quotient, remainder)."<br>
- | quo rem ql d div dh dnh dl qhi qlo j l hi lo r3 a t divDigitLength remDigitLength |<br>
- <primitive: 'primDigitDivNegative' module:'LargeIntegers'><br>
- arg = 0 ifTrue: [^ (ZeroDivide dividend: self) signal].<br>
- "TFEI added this line"<br>
- l := self digitLength - arg digitLength + 1.<br>
- l <= 0 ifTrue: [^ Array with: 0 with: self].<br>
- "shortcut against #highBit"<br>
- d := 8 - arg lastDigit highBitOfByte.<br>
- div := arg digitLshift: d.<br>
- divDigitLength := div digitLength + 1.<br>
- div := div growto: divDigitLength.<br>
- "shifts so high order word is >=128"<br>
- rem := self digitLshift: d.<br>
- rem digitLength = self digitLength ifTrue: [rem := rem growto: self digitLength + 1].<br>
- remDigitLength := rem digitLength.<br>
- "makes a copy and shifts"<br>
- quo := Integer new: l neg: ng.<br>
- dl := divDigitLength - 1.<br>
- "Last actual byte of data"<br>
- ql := l.<br>
- dh := div digitAt: dl.<br>
- dnh := dl = 1<br>
- ifTrue: [0]<br>
- ifFalse: [div digitAt: dl - 1].<br>
- 1 to: ql do: <br>
- [:k | <br>
- "maintain quo*arg+rem=self"<br>
- "Estimate rem/div by dividing the leading to bytes of rem by dh."<br>
- "The estimate is q = qhi*16+qlo, where qhi and qlo are nibbles."<br>
- j := remDigitLength + 1 - k.<br>
- "r1 := rem digitAt: j."<br>
- (rem digitAt: j)<br>
- = dh<br>
- ifTrue: [qhi := qlo := 15<br>
- "i.e. q=255"]<br>
- ifFalse: <br>
- ["Compute q = (r1,r2)//dh, t = (r1,r2)\\dh. <br>
- Note that r1,r2 are bytes, not nibbles. <br>
- Be careful not to generate intermediate results exceeding 13 <br>
- bits."<br>
- "r2 := (rem digitAt: j - 1)."<br>
- t := ((rem digitAt: j)<br>
- bitShift: 4)<br>
- + ((rem digitAt: j - 1)<br>
- bitShift: -4).<br>
- qhi := t // dh.<br>
- t := (t \\ dh bitShift: 4)<br>
- + ((rem digitAt: j - 1)<br>
- bitAnd: 15).<br>
- qlo := t // dh.<br>
- t := t \\ dh.<br>
- "Next compute (hi,lo) := q*dnh"<br>
- hi := qhi * dnh.<br>
- lo := qlo * dnh + ((hi bitAnd: 15)<br>
- bitShift: 4).<br>
- hi := (hi bitShift: -4)<br>
- + (lo bitShift: -8).<br>
- lo := lo bitAnd: 255.<br>
- "Correct overestimate of q. <br>
- Max of 2 iterations through loop -- see Knuth vol. 2"<br>
- r3 := j < 3<br>
- ifTrue: [0]<br>
- ifFalse: [rem digitAt: j - 2].<br>
- [(t < hi<br>
- or: [t = hi and: [r3 < lo]])<br>
- and: <br>
- ["i.e. (t,r3) < (hi,lo)"<br>
- qlo := qlo - 1.<br>
- lo := lo - dnh.<br>
- lo < 0<br>
- ifTrue: <br>
- [hi := hi - 1.<br>
- lo := lo + 256].<br>
- hi >= dh]]<br>
- whileTrue: [hi := hi - dh].<br>
- qlo < 0<br>
- ifTrue: <br>
- [qhi := qhi - 1.<br>
- qlo := qlo + 16]].<br>
- "Subtract q*div from rem"<br>
- l := j - dl.<br>
- a := 0.<br>
- 1 to: divDigitLength do: <br>
- [:i | <br>
- hi := (div digitAt: i)<br>
- * qhi.<br>
- lo := a + (rem digitAt: l) - ((hi bitAnd: 15)<br>
- bitShift: 4) - ((div digitAt: i)<br>
- * qlo).<br>
- rem digitAt: l put: lo - (lo // 256 * 256).<br>
- "sign-tolerant form of (lo bitAnd: 255)"<br>
- a := lo // 256 - (hi bitShift: -4).<br>
- l := l + 1].<br>
- a < 0<br>
- ifTrue: <br>
- ["Add div back into rem, decrease q by 1"<br>
- qlo := qlo - 1.<br>
- l := j - dl.<br>
- a := 0.<br>
- 1 to: divDigitLength do: <br>
- [:i | <br>
- a := (a bitShift: -8)<br>
- + (rem digitAt: l) + (div digitAt: i).<br>
- rem digitAt: l put: (a bitAnd: 255).<br>
- l := l + 1]].<br>
- quo digitAt: ql + 1 - k put: (qhi bitShift: 4)<br>
- + qlo].<br>
- rem := rem<br>
- digitRshift: d<br>
- bytes: 0<br>
- lookfirst: dl.<br>
- ^ Array with: quo with: rem!<br>
<br>
Item was added:<br>
+ ----- Method: Integer>>primDigitDiv:neg: (in category 'private') -----<br>
+ primDigitDiv: arg neg: ng <br>
+ "Answer with an array of (quotient, remainder)."<br>
+ | quo rem ql d div dh dnh dl qhi qlo j l hi lo r3 a t divDigitLength remDigitLength |<br>
+ <primitive: 'primDigitDivNegative' module:'LargeIntegers'><br>
+ arg = 0 ifTrue: [^ (ZeroDivide dividend: self) signal].<br>
+ "TFEI added this line"<br>
+ l := self digitLength - arg digitLength + 1.<br>
+ l <= 0 ifTrue: [^ Array with: 0 with: self].<br>
+ "shortcut against #highBit"<br>
+ d := 8 - arg lastDigit highBitOfByte.<br>
+ div := arg digitLshift: d.<br>
+ divDigitLength := div digitLength + 1.<br>
+ div := div growto: divDigitLength.<br>
+ "shifts so high order word is >=128"<br>
+ rem := self digitLshift: d.<br>
+ rem digitLength = self digitLength ifTrue: [rem := rem growto: self digitLength + 1].<br>
+ remDigitLength := rem digitLength.<br>
+ "makes a copy and shifts"<br>
+ quo := Integer new: l neg: ng.<br>
+ dl := divDigitLength - 1.<br>
+ "Last actual byte of data"<br>
+ ql := l.<br>
+ dh := div digitAt: dl.<br>
+ dnh := dl = 1<br>
+ ifTrue: [0]<br>
+ ifFalse: [div digitAt: dl - 1].<br>
+ 1 to: ql do: <br>
+ [:k | <br>
+ "maintain quo*arg+rem=self"<br>
+ "Estimate rem/div by dividing the leading to bytes of rem by dh."<br>
+ "The estimate is q = qhi*16+qlo, where qhi and qlo are nibbles."<br>
+ j := remDigitLength + 1 - k.<br>
+ "r1 := rem digitAt: j."<br>
+ (rem digitAt: j)<br>
+ = dh<br>
+ ifTrue: [qhi := qlo := 15<br>
+ "i.e. q=255"]<br>
+ ifFalse: <br>
+ ["Compute q = (r1,r2)//dh, t = (r1,r2)\\dh. <br>
+ Note that r1,r2 are bytes, not nibbles. <br>
+ Be careful not to generate intermediate results exceeding 13 <br>
+ bits."<br>
+ "r2 := (rem digitAt: j - 1)."<br>
+ t := ((rem digitAt: j)<br>
+ bitShift: 4)<br>
+ + ((rem digitAt: j - 1)<br>
+ bitShift: -4).<br>
+ qhi := t // dh.<br>
+ t := (t \\ dh bitShift: 4)<br>
+ + ((rem digitAt: j - 1)<br>
+ bitAnd: 15).<br>
+ qlo := t // dh.<br>
+ t := t \\ dh.<br>
+ "Next compute (hi,lo) := q*dnh"<br>
+ hi := qhi * dnh.<br>
+ lo := qlo * dnh + ((hi bitAnd: 15)<br>
+ bitShift: 4).<br>
+ hi := (hi bitShift: -4)<br>
+ + (lo bitShift: -8).<br>
+ lo := lo bitAnd: 255.<br>
+ "Correct overestimate of q. <br>
+ Max of 2 iterations through loop -- see Knuth vol. 2"<br>
+ r3 := j < 3<br>
+ ifTrue: [0]<br>
+ ifFalse: [rem digitAt: j - 2].<br>
+ [(t < hi<br>
+ or: [t = hi and: [r3 < lo]])<br>
+ and: <br>
+ ["i.e. (t,r3) < (hi,lo)"<br>
+ qlo := qlo - 1.<br>
+ lo := lo - dnh.<br>
+ lo < 0<br>
+ ifTrue: <br>
+ [hi := hi - 1.<br>
+ lo := lo + 256].<br>
+ hi >= dh]]<br>
+ whileTrue: [hi := hi - dh].<br>
+ qlo < 0<br>
+ ifTrue: <br>
+ [qhi := qhi - 1.<br>
+ qlo := qlo + 16]].<br>
+ "Subtract q*div from rem"<br>
+ l := j - dl.<br>
+ a := 0.<br>
+ 1 to: divDigitLength do: <br>
+ [:i | <br>
+ hi := (div digitAt: i)<br>
+ * qhi.<br>
+ lo := a + (rem digitAt: l) - ((hi bitAnd: 15)<br>
+ bitShift: 4) - ((div digitAt: i)<br>
+ * qlo).<br>
+ rem digitAt: l put: lo - (lo // 256 * 256).<br>
+ "sign-tolerant form of (lo bitAnd: 255)"<br>
+ a := lo // 256 - (hi bitShift: -4).<br>
+ l := l + 1].<br>
+ a < 0<br>
+ ifTrue: <br>
+ ["Add div back into rem, decrease q by 1"<br>
+ qlo := qlo - 1.<br>
+ l := j - dl.<br>
+ a := 0.<br>
+ 1 to: divDigitLength do: <br>
+ [:i | <br>
+ a := (a bitShift: -8)<br>
+ + (rem digitAt: l) + (div digitAt: i).<br>
+ rem digitAt: l put: (a bitAnd: 255).<br>
+ l := l + 1]].<br>
+ quo digitAt: ql + 1 - k put: (qhi bitShift: 4)<br>
+ + qlo].<br>
+ rem := rem<br>
+ digitRshift: d<br>
+ bytes: 0<br>
+ lookfirst: dl.<br>
+ ^ Array with: quo with: rem!<br>
<br>
Item was added:<br>
+ ----- Method: LargePositiveInteger>>digitDiv21: (in category 'private') -----<br>
+ digitDiv21: anInteger<br>
+ "This is part of the recursive division algorithm from Burnikel - Ziegler<br>
+ Divide a two limbs receiver by 1 limb dividend<br>
+ Each limb is decomposed in two halves of p bytes (8*p bits)<br>
+ so as to continue the recursion"<br>
+ <br>
+ | p qr1 qr2 |<br>
+ p := anInteger digitLength + 1 bitShift: -1.<br>
+ p <= 256 ifTrue: [^(self primDigitDiv: anInteger neg: false) collect: #normalize].<br>
+ qr1 := (self butLowestNDigits: p) digitDiv32: anInteger.<br>
+ qr2 := (self lowestNDigits: p) + (qr1 last bitShift: 8*p) digitDiv32: anInteger.<br>
+ qr2 at: 1 put: (qr2 at: 1) + ((qr1 at: 1) bitShift: 8*p).<br>
+ ^qr2!<br>
<br>
Item was added:<br>
+ ----- Method: LargePositiveInteger>>digitDiv32: (in category 'private') -----<br>
+ digitDiv32: anInteger<br>
+ "This is part of the recursive division algorithm from Burnikel - Ziegler<br>
+ Divide 3 limb (a2,a1,a0) by 2 limb (b1,b0).<br>
+ Each limb is made of p bytes (8*p bits).<br>
+ This step transforms the division problem into multiplication<br>
+ It must use the fastMultiply: to be worth the overhead costs."<br>
+ <br>
+ | a2 b1 d p q qr r |<br>
+ p := anInteger digitLength + 1 bitShift: -1.<br>
+ (a2 := self butLowestNDigits: 2*p) <br>
+ < (b1 := anInteger butLowestNDigits: p)<br>
+ ifTrue:<br>
+ [qr := (self butLowestNDigits: p) digitDiv21: b1.<br>
+ q := qr first.<br>
+ r := qr last]<br>
+ ifFalse:<br>
+ [q := (1 bitShift: 8*p) - 1.<br>
+ r := (self butLowestNDigits: p) - (b1 bitShift: 8*p) + b1].<br>
+ d := q fastMultiply: (anInteger lowestNDigits: p).<br>
+ r := (self lowestNDigits: p) + (r bitShift: 8*p) - d.<br>
+ [r < 0]<br>
+ whileTrue:<br>
+ [q := q - 1.<br>
+ r := r + anInteger].<br>
+ ^Array with: q with: r<br>
+ !<br>
<br>
Item was added:<br>
+ ----- Method: LargePositiveInteger>>digitDiv:neg: (in category 'private') -----<br>
+ digitDiv: anInteger neg: aBoolean <br>
+ "If length is worth, engage a recursive divide and conquer strategy"<br>
+ | qr |<br>
+ (anInteger digitLength <= 256<br>
+ or: [self digitLength <= anInteger digitLength])<br>
+ ifTrue: [^ self primDigitDiv: anInteger neg: aBoolean].<br>
+ qr := self abs recursiveDigitDiv: anInteger abs.<br>
+ ^ aBoolean<br>
+ ifTrue: [qr collect: #negated]<br>
+ ifFalse: [qr]!<br>
<br>
Item was changed:<br>
----- Method: LargePositiveInteger>>printOn:base: (in category 'printing') -----<br>
printOn: aStream base: b<br>
"Append a representation of this number in base b on aStream.<br>
In order to reduce cost of LargePositiveInteger ops, split the number in approximately two equal parts in number of digits."<br>
<br>
+ | halfDigits halfPower head tail nDigitsUnderestimate qr |<br>
- | halfDigits halfPower head tail nDigitsUnderestimate |<br>
"Don't engage any arithmetic if not normalized"<br>
(self digitLength = 0 or: [(self digitAt: self digitLength) = 0]) ifTrue: [^self normalize printOn: aStream base: b].<br>
<br>
nDigitsUnderestimate := b = 10<br>
ifTrue: [((self highBit - 1) * 1233 >> 12) + 1. "This is because (2 log)/(10 log)*4096 is slightly greater than 1233"]<br>
ifFalse: [self highBit quo: b highBit].<br>
<br>
"splitting digits with a whole power of two is more efficient"<br>
halfDigits := 1 bitShift: nDigitsUnderestimate highBit - 2.<br>
<br>
halfDigits <= 1<br>
ifTrue: ["Hmmm, this could happen only in case of a huge base b... Let lower level fail"<br>
^self printOn: aStream base: b nDigits: (self numberOfDigitsInBase: b)].<br>
<br>
"Separate in two halves, head and tail"<br>
halfPower := b raisedToInteger: halfDigits.<br>
+ qr := self digitDiv: halfPower neg: self negative.<br>
+ head := qr first normalize.<br>
+ tail := qr last normalize.<br>
- head := self quo: halfPower.<br>
- tail := self - (head * halfPower).<br>
<br>
"print head"<br>
head printOn: aStream base: b.<br>
<br>
"print tail without the overhead to count the digits"<br>
tail printOn: aStream base: b nDigits: halfDigits!<br>
<br>
Item was changed:<br>
----- Method: LargePositiveInteger>>printOn:base:nDigits: (in category 'printing') -----<br>
printOn: aStream base: b nDigits: n<br>
"Append a representation of this number in base b on aStream using n digits.<br>
In order to reduce cost of LargePositiveInteger ops, split the number of digts approximatily in two<br>
Should be invoked with: 0 <= self < (b raisedToInteger: n)"<br>
<br>
+ | halfPower half head tail qr |<br>
- | halfPower half head tail |<br>
n <= 1 ifTrue: [<br>
n <= 0 ifTrue: [self error: 'Number of digits n should be > 0'].<br>
<br>
"Note: this is to stop an infinite loop if one ever attempts to print with a huge base<br>
This can happen because choice was to not hardcode any limit for base b<br>
We let Character>>#digitValue: fail"<br>
^aStream nextPut: (Character digitValue: self) ].<br>
halfPower := n bitShift: -1.<br>
half := b raisedToInteger: halfPower.<br>
+ qr := self digitDiv: half neg: self negative.<br>
+ head := qr first normalize.<br>
+ tail := qr last normalize.<br>
- head := self quo: half.<br>
- tail := self - (head * half).<br>
head printOn: aStream base: b nDigits: n - halfPower.<br>
tail printOn: aStream base: b nDigits: halfPower!<br>
<br>
Item was added:<br>
+ ----- Method: LargePositiveInteger>>recursiveDigitDiv: (in category 'private') -----<br>
+ recursiveDigitDiv: anInteger<br>
+ "This is the recursive division algorithm from Burnikel - Ziegler<br>
+ See Fast Recursive Division - Christoph Burnikel, Joachim Ziegler<br>
+ Research Report MPI-I-98-1-022, MPI Saarbrucken, Oct 1998<br>
+ <a href="https://pure.mpg.de/rest/items/item_1819444_4/component/file_2599480/content" rel="noreferrer" target="_blank">https://pure.mpg.de/rest/items/item_1819444_4/component/file_2599480/content</a>"<br>
+ <br>
+ | s m t a b z qr q i |<br>
+ "round digits up to next power of 2"<br>
+ s := anInteger digitLength.<br>
+ m := 1 bitShift: (s - 1) highBit.<br>
+ "shift so that leading bit of leading byte be 1, and digitLength power of two"<br>
+ s := m * 8 - anInteger highBit.<br>
+ a := self bitShift: s.<br>
+ b := anInteger bitShift: s.<br>
+ <br>
+ "Decompose a into t limbs - each limb have m bytes<br>
+ choose t such that leading bit of leading limb of a be 0"<br>
+ t := (a highBit + 1 / (m * 8)) ceiling.<br>
+ z := a butLowestNDigits: t - 2 * m.<br>
+ i := t - 2.<br>
+ q := 0.<br>
+ "and do a division of two limb by 1 limb b for each pair of limb of a"<br>
+ [qr := z digitDiv21: b.<br>
+ q := (q bitShift: 8*m) + qr first. "Note: this naive recomposition of q cost O(t^2) - it is possible in O(t log(t))"<br>
+ (i := i - 1) >= 0] whileTrue:<br>
+ [z := (qr last bitShift: 8*m) + (a copyDigitsFrom: i * m + 1 to: i + 1 * m)].<br>
+ ^Array with: q with: (qr last bitShift: s negated)!<br>
<br>
Item was changed:<br>
----- Method: LargePositiveInteger>>sqrtRem (in category 'mathematical functions') -----<br>
sqrtRem<br>
"Like super, but use a divide and conquer method to perform this operation.<br>
See Paul Zimmermann. Karatsuba Square Root. [Research Report] RR-3805, INRIA. 1999, pp.8. <inria-00072854><br>
<a href="https://hal.inria.fr/inria-00072854/PDF/RR-3805.pdf" rel="noreferrer" target="_blank">https://hal.inria.fr/inria-00072854/PDF/RR-3805.pdf</a>"<br>
<br>
+ | n qr q s r sr high mid low |<br>
- | n qr s r sr high mid low |<br>
n := self digitLength bitShift: -2.<br>
n >= 16 ifFalse: [^super sqrtRem].<br>
high := self butLowestNDigits: n * 2.<br>
mid := self copyDigitsFrom: n + 1 to: n * 2.<br>
low := self lowestNDigits: n.<br>
<br>
sr := high sqrtRem.<br>
qr := (sr last bitShift: 8 * n) + mid digitDiv: (sr first bitShift: 1) neg: false.<br>
+ q := qr first normalize.<br>
+ s := (sr first bitShift: 8 * n) + q.<br>
+ r := (qr last normalize bitShift: 8 * n) + low - q squared.<br>
- s := (sr first bitShift: 8 * n) + qr first.<br>
- r := (qr last bitShift: 8 * n) + low - qr first squared.<br>
r negative<br>
ifTrue:<br>
[r := (s bitShift: 1) + r - 1.<br>
s := s - 1].<br>
sr at: 1 put: s; at: 2 put: r.<br>
^sr<br>
!<br>
<br>
Item was changed:<br>
----- Method: LargePositiveInteger>>squared (in category 'mathematical functions') -----<br>
squared<br>
"Eventually use a divide and conquer algorithm to perform the multiplication"<br>
<br>
(self digitLength >= 400) ifFalse: [^self * self].<br>
+ (self digitLength >= 800) ifFalse: [^self squaredKaratsuba].<br>
+ ^self squaredToom4!<br>
- (self digitLength >= 1600) ifFalse: [^self squaredKaratsuba].<br>
- ^self squaredToom3!<br>
<br>
Item was added:<br>
+ ----- Method: LargePositiveInteger>>squaredToom4 (in category 'mathematical functions') -----<br>
+ squaredToom4<br>
+ "Use a 4-way Toom-Cook divide and conquer algorithm to perform the multiplication.<br>
+ See Asymmetric Squaring Formulae Jaewook Chung and M. Anwar Hasan<br>
+ <a href="https://www.lirmm.fr/arith18/papers/Chung-Squaring.pdf" rel="noreferrer" target="_blank">https://www.lirmm.fr/arith18/papers/Chung-Squaring.pdf</a>"<br>
+ <br>
+ | p a0 a1 a2 a3 a02 a13 s0 s1 s2 s3 s4 s5 s6 t2 t3 |<br>
+ "divide in 4 parts"<br>
+ p := (self digitLength + 3 bitShift: -2) bitClear: 2r11.<br>
+ a3 := self butLowestNDigits: p * 3.<br>
+ a2 := self copyDigitsFrom: p * 2 + 1 to: p * 3.<br>
+ a1 := self copyDigitsFrom: p + 1 to: p * 2.<br>
+ a0 := self lowestNDigits: p.<br>
+ <br>
+ "Toom-4 trick: 7 multiplications instead of 16"<br>
+ a02 := a0 - a2.<br>
+ a13 := a1 - a3.<br>
+ s0 := a0 squared.<br>
+ s1 := (a0 fastMultiply: a1) bitShift: 1.<br>
+ s2 := (a02 + a13) fastMultiply: (a02 - a13).<br>
+ s3 := ((a0 + a1) + (a2 + a3)) squared.<br>
+ s4 := (a02 fastMultiply: a13) bitShift: 1.<br>
+ s5 := (a3 fastMultiply: a2) bitShift: 1.<br>
+ s6 := a3 squared.<br>
+ <br>
+ "Interpolation"<br>
+ t2 := s1 + s5.<br>
+ t3 := (s2 + s3 + s4 bitShift: -1) - t2.<br>
+ s3 := t2 - s4.<br>
+ s4 := t3 - s0.<br>
+ s2 := t3 - s2 - s6.<br>
+ <br>
+ "Sum the parts of decomposition"<br>
+ ^s0 + (s1 bitShift: 8*p) + (s2 + (s3 bitShift: 8*p) bitShift: 16*p)<br>
+ +(s4 + (s5 bitShift: 8*p) + (s6 bitShift: 16*p) bitShift: 32*p)<br>
+ <br>
+ "<br>
+ | a |<br>
+ a := 770 factorial-1.<br>
+ a digitLength.<br>
+ [a * a - a squaredToom4 = 0] assert.<br>
+ [Smalltalk garbageCollect.<br>
+ [1000 timesRepeat: [a squaredToom4]] timeToRun] value /<br>
+ [Smalltalk garbageCollect.<br>
+ [1000 timesRepeat: [a squaredKaratsuba]] timeToRun] value asFloat<br>
+ "!<br>
<br>
<br>
</blockquote></div>