The Trunk: Kernel-nice.741.mcz

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!