[Pkg] The Trunk: Kernel-nice.741.mcz

[commits at source.squeak.org]

Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernel-nice.741.mcz

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

Name: Kernel-nice.741
Author: nice
Time: 24 February 2013, 5:32:56.978 pm
UUID: 15a282b0-84ac-4723-b44c-7dd91e9aebf0
Ancestors: Kernel-eem.740

Reduce the number of hardcoded constants used by #raisedTo:modulo:
This will enable smooth transition to next generation of LargeIntegersPlugin.

=============== Diff against Kernel-eem.740 ===============

Item was added:
+ ----- Method: Integer>>montgomeryDigitBase (in category 'private') -----
+ montgomeryDigitBase
+ 	"Answer the base used by Montgomery algorithm."
+ 	^1 << self montgomeryDigitLength!

Item was added:
+ ----- Method: Integer>>montgomeryDigitLength (in category 'private') -----
+ montgomeryDigitLength
+ 	"Answer the number of bits composing a digit in Montgomery algorithm.
+ 	Primitive use either 8 or 32 bits digits"
+ 	<primitive: 'primMontgomeryDigitLength' module:'LargeIntegers'>
+ 	^8 "Legacy plugin which did not have this primitive did use 8 bits digits"!

Item was added:
+ ----- Method: Integer>>montgomeryDigitMax (in category 'private') -----
+ montgomeryDigitMax
+ 	"Answer the maximum value of a digit used in Montgomery algorithm."
+ 	
+ 	^1 << self montgomeryDigitLength - 1!

Item was added:
+ ----- Method: Integer>>montgomeryNumberOfDigits (in category 'private') -----
+ montgomeryNumberOfDigits
+ 	"Answer the number of montgomery digits required to represent the receiver."
+ 	^self digitLength * 8 + (self montgomeryDigitLength - 1) // self montgomeryDigitLength!

Item was changed:
  ----- Method: Integer>>montgomeryRaisedTo:times:modulo:mInvModB: (in category 'private') -----
  montgomeryRaisedTo: n times: y modulo: m mInvModB: mInv
  	"Private - do a Montgomery exponentiation of self modulo m.
  	The operation is equivalent to (self/y raisedTo: n)*y \\ m,
+ 	with y is (b raisedTo: m montgomeryNumberOfDigits),
+ 	with (m bitAnd: b-1) * mInv \\ b = (b-1)
+ 	with b = self montgomeryDigitBase (either 1<<8 or 1<<32)"
- 	with y is (256 raisedTo: m digitLength),
- 	with (m bitAnd: 255) * mInv \\ 256 = 255."
  	
  	| pow j k w index oddPowersOfSelf square |
  	
  	"Precompute powers of self for odd bit patterns xxxx1 up to length w + 1.
  	The width w is chosen with respect to the total bit length of n,
  	such that each bit pattern will on average be encoutered P times in the whole bit sequence of n.
  	This costs (2 raisedTo: w) multiplications, but more will be saved later (see below)."
  	k := n highBit.
  	w := (k highBit - 1 >> 1 min: 16) max: 1.
  	oddPowersOfSelf := Array new: 1 << w.
  	oddPowersOfSelf at: 1 put: (pow := self).
  	square := self montgomeryTimes: self modulo: m mInvModB: mInv.
  	2 to: oddPowersOfSelf size do: [:i | pow := oddPowersOfSelf at: i put: (pow montgomeryTimes: square modulo: m mInvModB: mInv)].
  	
  	"Now exponentiate by searching precomputed bit patterns with a sliding window"
  	pow := y.
  	[k > 0]
  		whileTrue:
  			[pow := pow montgomeryTimes: pow modulo: m mInvModB: mInv.
  			"Skip bits set to zero (the sliding window)"
  			(n bitAt: k) = 0
  				ifFalse:
  					["Find longest odd bit pattern up to window length (w + 1)"
  					j := k - w max: 1.
  					[j < k and: [(n bitAt: j) = 0]] whileTrue: [j := j + 1].
  					"We found a bit pattern of length k-j+1;
  					perform the square powers for each bit
  					(same cost as bitwise algorithm);
  					compute the index of this bit pattern in the precomputed powers."
  					index := 0.
  					[k > j] whileTrue:
  						[pow := pow montgomeryTimes: pow modulo: m mInvModB: mInv.
  						index := index << 1 + (n bitAt: k).
  						k := k - 1].
  					"Perform a single multiplication for the whole bit pattern.
  					This saves up to (k-j) multiplications versus a naive algorithm operating bit by bit"
  					pow := pow montgomeryTimes: (oddPowersOfSelf at: index + 1) modulo: m mInvModB: mInv].
  			k := k - 1].
  	^pow!

Item was changed:
  ----- Method: Integer>>montgomeryTimes:modulo:mInvModB: (in category 'private') -----
  montgomeryTimes: a modulo: m mInvModB: mInv
  	"Answer the result of a Montgomery multiplication
+ 	self * a * (b raisedTo: m montgomeryNumberOfDigits) inv \\ m
- 	self * a * (256 raisedTo: m digitLength) inv \\ m
  	NOTE: it is assumed that:
+ 	self montgomeryNumberOfDigits <= m montgomeryNumberOfDigits
+ 	a montgomeryNumberOfDigits <= m montgomeryNumberOfDigits
+ 	mInv * m \\ b = (-1 \\ b) = (b-1) (this implies m odd)
+ 	where b = self montgomeryDigitBase
- 	self digitLength <= m digitLength
- 	a digitLength <= m digitLength
- 	mInv * m \\ 256 = (-1 \\ 256) = 255 (this implies m odd)
  	
  	Answer nil in case of absent plugin or other failure."
  	
  	<primitive: 'primMontgomeryTimesModulo' module:'LargeIntegers'>
  	^nil!

Item was changed:
  ----- Method: Integer>>raisedTo:modulo: (in category 'mathematical functions') -----
  raisedTo: n modulo: m
  	"Answer the modular exponential.
  	Note: this implementation is optimized for case of large integers raised to large powers."
  	| a s mInv |
  	n = 0 ifTrue: [^1].
  	(self >= m or: [self < 0]) ifTrue: [^self \\ m raisedTo: n modulo: m].
  	n < 0 ifTrue: [^(self reciprocalModulo: m) raisedTo: n negated modulo: m].
  	(n < 4096 or: [m even])
  		ifTrue:
  			["Overhead of Montgomery method might cost more than naive divisions, use naive"
  			^self slidingLeftRightRaisedTo: n modulo: m].
  	
+ 	mInv := self montgomeryDigitBase - ((m bitAnd: self montgomeryDigitMax) reciprocalModulo: self montgomeryDigitBase).
- 	mInv := 256 - ((m bitAnd: 255) reciprocalModulo: 256).
   
+ 	"Initialize the result to R=self montgomeryDigitModulo raisedTo: m montgomeryNumberOfDigits"
+ 	a := (1 bitShift: m montgomeryNumberOfDigits * m montgomeryDigitLength) \\ m.
- 	"Initialize the result to R=256 raisedTo: m digitLength"
- 	a := (1 bitShift: m digitLength*8) \\ m.
  	
  	"Montgomerize self (multiply by R)"
  	(s := self montgomeryTimes: (a*a \\ m) modulo: m mInvModB: mInv)
  		ifNil:
  			["No Montgomery primitive available ? fallback to naive divisions"
  			^self slidingLeftRightRaisedTo: n modulo: m].
  
  	"Exponentiate self*R"
  	a := s montgomeryRaisedTo: n times: a modulo: m mInvModB: mInv.
  
  	"Demontgomerize the result (divide by R)"
  	^a montgomeryTimes: 1 modulo: m mInvModB: mInv!