Montgomery reduction algorithm with Python
[Primes Home][Home]
The Montgomery reduction algorithm is a fast way to perform multiplications in the form of \(x y \pmod N\) or \(x^y \pmod N\). In this case we will take a prime number of \(N\) and determine \(x \times y \pmod N\) and \(x^y \pmod N\) [Golang version]:
|
Outline
With the RSA and the Diffie-Hellman method we perform large exponential calculations, such as:
\(C = M^e \pmod N\)
and where we will continually multiply large integers by an exponent to get a result. If we were to just calculate \(x^y\) and then take \(\pmod n\) it would take a while to produce the result. Thus we use Montgomery modular multiplication, and which significantly reduces the time to compute the result. In a traditional multiplication of two value (\(x\) and \(y\)) for a modulus of \(N\), we multiply \(x\) times \(y\) and then divide by \(N\) to find the remainder. The number of bits in the multiplication with this be the number of bits in \(x\) added to the number of bits in \(y\). In Montgomery reduction we add multiples of \(N\) in order to simplify the multiplication.
An example of this is here, and a sample run for \(x=10\), \(y=5\) and \(N=29\) is:
a=10, b=5, p=29 Result: 10*5 (mod 29) = 21 Traditional method result = 21 Result: 10^5 (mod 29) = 8 Traditional method result = 8
In this case we get \(50 \pmod {29}\) which is 21, and \(10^5 \pmod {29}\) which is 8.
The sample code is:
import sys from montred import MontgomeryReducer a=99 b=5 p=83 if (len(sys.argv)>1): a=int(sys.argv[1]) if (len(sys.argv)>2): b=int(sys.argv[2]) if (len(sys.argv)>3): p=int(sys.argv[3]) print (f"x={a}, y={b}, p={p}") mr = MontgomeryReducer(p) aval=mr.convert_in(a) bval=mr.convert_in(b) res=mr.convert_out(mr.multiply(aval,bval)) print (f"\n({a}x{b}) (mod {p})={res}") aval=mr.convert_in(a) res=mr.convert_out(mr.pow(aval,b)) print (f"\n({a}^{b}) (mod {p})={res}") try: inv_b=mr.reciprocal_mod(b,p) print (f"\nInv y (mod {p})={inv_b}") inv_b=mr.convert_in(inv_b) res=mr.convert_out(mr.multiply(aval,inv_b)) print (f"\n({a}/{b}) (mod {p})={res}") except: print(f"\nReciprocal of {b} (mod {p}) does not exist, so we cannot divide.")
and Montgomery reduction:
# https://www.nayuki.io/res/montgomery-reduction-algorithm/montgomery-reducer.py import math class MontgomeryReducer: def __init__(self, mod): # Modulus if mod < 3 or mod % 2 == 0: raise ValueError("Modulus must be an odd number at least 3") self.modulus = mod # Reducer self.reducerbits = (mod.bit_length() // 8 + 1) * 8 # This is a multiple of 8 self.reducer = 1 << self.reducerbits # This is a power of 256 self.mask = self.reducer - 1 assert self.reducer > mod and math.gcd(self.reducer, mod) == 1 # Other computed numbers self.reciprocal = MontgomeryReducer.reciprocal_mod(self.reducer % mod, mod) self.factor = (self.reducer * self.reciprocal - 1) // mod self.convertedone = self.reducer % mod # The range of x is unlimited def convert_in(self, x): return (x << self.reducerbits) % self.modulus # The range of x is unlimited def convert_out(self, x): return (x * self.reciprocal) % self.modulus # Inputs and output are in Montgomery form and in the range [0, modulus) def multiply(self, x, y): mod = self.modulus assert 0 <= x < mod and 0 <= y < mod product = x * y temp = ((product & self.mask) * self.factor) & self.mask reduced = (product + temp * mod) >> self.reducerbits result = reduced if (reduced < mod) else (reduced - mod) assert 0 <= result < mod return result # Input x (base) and output (power) are in Montgomery form and in the range [0, modulus); input y (exponent) is in standard form def pow(self, x, y): assert 0 <= x < self.modulus if y < 0: raise ValueError("Negative exponent") z = self.convertedone while y != 0: if y & 1 != 0: z = self.multiply(z, x) x = self.multiply(x, x) y >>= 1 return z @staticmethod def reciprocal_mod(x, mod): # Based on a simplification of the extended Euclidean algorithm assert mod > 0 and 0 <= x < mod y = x x = mod a = 0 b = 1 while y != 0: a, b = b, a - x // y * b x, y = y, x % y if x == 1: return a % mod else: raise ValueError("Reciprocal does not exist")
A sample run of \(30^{19} \pmod {53}\) [Try!] is:
x=30, y=19, p=53 (30x19) (mod 53)=40 (30^19) (mod 53)=23 Inv y (mod 53)=14 (30/19) (mod 53)=49
In this case:
- Multiply: \(30 \times 19 \pmod {53} = 570 \pmod {53} = 40\)
- Exponent: \(30^{19} \pmod {53} = 11622614670000000000000000000 \pmod {53} = 23\)
- Divide: \(19^{-1} \pmod {53} = 14\). \(30 \times 14 \pmod {53} = 420 \pmod {53} = 49\). We can check this with \(49 \times 19 \pmod {53}=30\)
Another example is:
x=99, y=5, p=997 (99x5) (mod 997)=495 (99^5) (mod 997)=47 Inv y (mod 997)=399 (99/5) (mod 997)=618
In this case:
- Multiply: \(99 \times 5 \pmod {997} = 495 \pmod {997} = 495\)
- Exponent: \(99^{5} \pmod {997} = 9509900499 \pmod {997} = 47\)
- Divide: \(5^{-1} \pmod {997} = 399\). \(99 \times 399 \pmod {997} = 39501 \pmod {997} = 618\). We can check this with \(618 \times 5 \pmod {997}=99\)
Presentation
The following is a presentation on the method [slides]:
References
[1] Montgomery, Peter L. "Modular multiplication without trial division." Mathematics of computation 44.170 (1985): 519-521.