Montgomery reduction: Difference between revisions
(Added back online viewing book link (shouldn't expect readers to fork over $70 just to work on this task. Unless the Amazon link carries commissions for this site, it needs not to exist); restore subscript i as in the book) |
(Go solution) |
||
Line 103:
return 0;
}</lang>
=={{header|Go}}==
<lang go>package main
import (
"fmt"
"math/big"
"math/rand"
"time"
)
// mont holds numbers useful for working in Mongomery representation.
type mont struct {
n uint // m.BitLen()
m *big.Int // modulus, must be odd
r2 *big.Int // (1<<2n) mod m
}
// constructor
func newMont(m *big.Int) *mont {
if m.Bit(0) != 1 {
return nil
}
n := uint(m.BitLen())
x := big.NewInt(1)
x.Sub(x.Lsh(x, n), m)
return &mont{n, new(big.Int).Set(m), x.Mod(x.Mul(x, x), m)}
}
// Montgomery reduction algorithm
func (m mont) reduce(t *big.Int) *big.Int {
a := new(big.Int).Set(t)
for i := uint(0); i < m.n; i++ {
if a.Bit(0) == 1 {
a.Add(a, m.m)
}
a.Rsh(a, 1)
}
if a.Cmp(m.m) >= 0 {
a.Sub(a, m.m)
}
return a
}
// example use:
func main() {
const n = 100 // bit length for numbers in example
// generate random n-bit odd number for modulus m
rnd := rand.New(rand.NewSource(time.Now().UnixNano()))
one := big.NewInt(1)
r1 := new(big.Int).Lsh(one, n-1)
r2 := new(big.Int).Lsh(one, n-2)
m := new(big.Int)
m.Or(r1, m.Or(m.Lsh(m.Rand(rnd, r2), 1), one))
// make Montgomery reduction object around m
mr := newMont(m)
// generate a couple more numbers in the range 0..m.
// these are numbers we will do some computations on, mod m.
x1 := new(big.Int).Rand(rnd, m)
x2 := new(big.Int).Rand(rnd, m)
// t1, t2 are examples of T, from the task description.
// Generated this way, they will be in the range 0..m^2, and so < mR.
t1 := new(big.Int).Mul(x1, mr.r2)
t2 := new(big.Int).Mul(x2, mr.r2)
// reduce. r1 and r2 are now montgomery representations of x1 and x2.
r1 = mr.reduce(t1)
r2 = mr.reduce(t2)
// this is the end of what is described in the task so far.
fmt.Println("b: 2")
fmt.Println("n: ", mr.n)
fmt.Println("r: ", new(big.Int).Lsh(one, mr.n))
fmt.Println("m: ", mr.m)
fmt.Println("t1:", t1)
fmt.Println("t2:", t2)
fmt.Println("r1:", r1)
fmt.Println("r2:", r2)
// but now demonstrate that it works:
fmt.Println()
fmt.Println("Original x1: ", x1)
fmt.Println("Recovererd from r1:", mr.reduce(r1))
fmt.Println("Original x2: ", x2)
fmt.Println("Recovererd from r2:", mr.reduce(r2))
// and demonstrate a use:
show("\nLibrary:", func() *big.Int {
return new(big.Int).Exp(x1, x2, m)
})
show("\nMontgomery:", func() *big.Int {
// this is the modular exponentiation algorithm, except we call
// mont.reduce instead of using a mod function.
prod := mr.reduce(mr.r2) // 1
base := mr.reduce(t1.Mul(x1, mr.r2)) // x1^1
exp := new(big.Int).Set(x2) // not reduced
for exp.BitLen() > 0 {
if exp.Bit(0) == 1 {
prod = mr.reduce(prod.Mul(prod, base))
}
exp.Rsh(exp, 1)
base = mr.reduce(base.Mul(base, base))
}
return mr.reduce(prod)
})
}
func show(heading string, f func() *big.Int) {
fmt.Println(heading)
t0 := time.Now()
fmt.Println("x1 ^ x2 mod m:", f())
fmt.Println(time.Now().Sub(t0))
}</lang>
Output:
<pre>
b: 2
n: 100
r: 1267650600228229401496703205376
m: 1191151171693032142151966564621
t1: 138602318824179275477121611740471794618999274340657947731554
t2: 397645077922552596057001924720561733079744309909810064024787
r1: 1073769066977456952449715239458
r2: 460620928979319347925360938242
Original x1: 141982853055250888109454150154
Recovererd from r1: 141982853055250888109454150154
Original x2: 407343709295665106703621643087
Recovererd from r2: 407343709295665106703621643087
Library:
x1 ^ x2 mod m: 599840174322403511400105423249
231us
Montgomery:
x1 ^ x2 mod m: 599840174322403511400105423249
2.316ms
</pre>
|
Revision as of 22:24, 23 December 2011
Implement the Montgomery reduction algorithm, as explained in "Handbook of Applied Cryptography, Section 14.3.2, page 600. Montgomery reduction calculates , without having to divide by .
- Let be a positive integer, and and integers such that , , and .
- is usually chosen as , where = base (radix) in which the numbers in the calculation as represented in (so in ‘normal’ paper arithmetic, for computer implementations) and = number of digits in base
- The numbers ( digits long), ( digits long), , , are known entities, a number (often represented as
m_dash
in code) = is precomputed.
See <amazon id=0849385237>[%buy% Handbook of Applied Cryptography]</amazon> for brief introduction to theory and numerical example in radix 10. Individual chapters of the book can be viewed online as provided by the authors.
Algorithm:
A ← T (temporary variable) For i from 0 to (n-1) do the following: ui ← ai* m' mod b // ai is the ith digit of A, ui is a single digit number in radix b A ← A + ui*m*bi A ← A/bn if A >= m, A ← A - m Return (A)
C++
<lang cpp>#include<iostream>
- include<conio.h>
using namespace std; typedef unsigned long ulong;
int ith_digit_finder(long long n, long b, long i){
/** n = number whose digits we need to extract b = radix in which the number if represented i = the ith bit (ie, index of the bit that needs to be extracted) **/ while(i>0){ n/=b; i--; } return (n%b);
}
long eeuclid(long m, long b, long *inverse){ /// eeuclid( modulus, num whose inv is to be found, variable to put inverse )
/// Algorithm used from Stallings book long A1 = 1, A2 = 0, A3 = m, B1 = 0, B2 = 1, B3 = b, T1, T2, T3, Q;
cout<<endl<<"eeuclid() started"<<endl;
while(1){ if(B3 == 0){ *inverse = 0; return A3; // A3 = gcd(m,b) }
if(B3 == 1){ *inverse = B2; // B2 = b^-1 mod m return B3; // A3 = gcd(m,b) }
Q = A3/B3;
T1 = A1 - Q*B1; T2 = A2 - Q*B2; T3 = A3 - Q*B3;
A1 = B1; A2 = B2; A3 = B3; B1 = T1; B2 = T2; B3 = T3;
} cout<<endl<<"ending eeuclid() "<<endl;
}
long long mon_red(long m, long m_dash, long T, int n, long b = 2){ /**
m = modulus m_dash = m' = -m^-1 mod b T = number whose modular reduction is needed, the o/p of the function is TR^-1 mod m n = number of bits in m (2n is the number of bits in T) b = radix used (for practical implementations, is equal to 2, which is the default value)
- /
long long A,ui, temp, Ai; // Ai is the ith bit of A, need not be llong long probably if( m_dash < 0 ) m_dash = m_dash + b; A = T; for(int i = 0; i<n; i++){ /// ui = ( (A%b)*m_dash ) % b; // step 2.1; A%b gives ai (MISTAKE -- A%b will always give the last digit of A if A is represented in base b); hence we need the function ith_digit_finder() Ai = ith_digit_finder(A, b, i); ui = ( ( Ai % b) * m_dash ) % b; temp = ui*m*power(b, i); A = A + temp; } A = A/power(b, n); if(A >= m) A = A - m; return A;
}
int main(){
long a, b, c, d=0, e, inverse = 0; cout<<"m >> "; cin >> a; cout<<"T >> "; cin>>b; cout<<"Radix b >> "; cin>>c; eeuclid(c, a, &d); // eeuclid( modulus, num whose inverse is to be found, address of variable which is to store inverse) e = mon_red(a, -d, b, length_finder(a, c), c); cout<<"Montgomery domain representation = "<<e; return 0;
}</lang>
Go
<lang go>package main
import (
"fmt" "math/big" "math/rand" "time"
)
// mont holds numbers useful for working in Mongomery representation. type mont struct {
n uint // m.BitLen() m *big.Int // modulus, must be odd r2 *big.Int // (1<<2n) mod m
}
// constructor func newMont(m *big.Int) *mont {
if m.Bit(0) != 1 { return nil } n := uint(m.BitLen()) x := big.NewInt(1) x.Sub(x.Lsh(x, n), m) return &mont{n, new(big.Int).Set(m), x.Mod(x.Mul(x, x), m)}
}
// Montgomery reduction algorithm func (m mont) reduce(t *big.Int) *big.Int {
a := new(big.Int).Set(t) for i := uint(0); i < m.n; i++ { if a.Bit(0) == 1 { a.Add(a, m.m) } a.Rsh(a, 1) } if a.Cmp(m.m) >= 0 { a.Sub(a, m.m) } return a
}
// example use: func main() {
const n = 100 // bit length for numbers in example
// generate random n-bit odd number for modulus m rnd := rand.New(rand.NewSource(time.Now().UnixNano())) one := big.NewInt(1) r1 := new(big.Int).Lsh(one, n-1) r2 := new(big.Int).Lsh(one, n-2) m := new(big.Int) m.Or(r1, m.Or(m.Lsh(m.Rand(rnd, r2), 1), one))
// make Montgomery reduction object around m mr := newMont(m)
// generate a couple more numbers in the range 0..m. // these are numbers we will do some computations on, mod m. x1 := new(big.Int).Rand(rnd, m) x2 := new(big.Int).Rand(rnd, m)
// t1, t2 are examples of T, from the task description. // Generated this way, they will be in the range 0..m^2, and so < mR. t1 := new(big.Int).Mul(x1, mr.r2) t2 := new(big.Int).Mul(x2, mr.r2)
// reduce. r1 and r2 are now montgomery representations of x1 and x2. r1 = mr.reduce(t1) r2 = mr.reduce(t2)
// this is the end of what is described in the task so far. fmt.Println("b: 2") fmt.Println("n: ", mr.n) fmt.Println("r: ", new(big.Int).Lsh(one, mr.n)) fmt.Println("m: ", mr.m) fmt.Println("t1:", t1) fmt.Println("t2:", t2) fmt.Println("r1:", r1) fmt.Println("r2:", r2)
// but now demonstrate that it works: fmt.Println() fmt.Println("Original x1: ", x1) fmt.Println("Recovererd from r1:", mr.reduce(r1)) fmt.Println("Original x2: ", x2) fmt.Println("Recovererd from r2:", mr.reduce(r2))
// and demonstrate a use: show("\nLibrary:", func() *big.Int { return new(big.Int).Exp(x1, x2, m) }) show("\nMontgomery:", func() *big.Int { // this is the modular exponentiation algorithm, except we call // mont.reduce instead of using a mod function. prod := mr.reduce(mr.r2) // 1 base := mr.reduce(t1.Mul(x1, mr.r2)) // x1^1 exp := new(big.Int).Set(x2) // not reduced for exp.BitLen() > 0 { if exp.Bit(0) == 1 { prod = mr.reduce(prod.Mul(prod, base)) } exp.Rsh(exp, 1) base = mr.reduce(base.Mul(base, base)) } return mr.reduce(prod) })
}
func show(heading string, f func() *big.Int) {
fmt.Println(heading) t0 := time.Now() fmt.Println("x1 ^ x2 mod m:", f()) fmt.Println(time.Now().Sub(t0))
}</lang> Output:
b: 2 n: 100 r: 1267650600228229401496703205376 m: 1191151171693032142151966564621 t1: 138602318824179275477121611740471794618999274340657947731554 t2: 397645077922552596057001924720561733079744309909810064024787 r1: 1073769066977456952449715239458 r2: 460620928979319347925360938242 Original x1: 141982853055250888109454150154 Recovererd from r1: 141982853055250888109454150154 Original x2: 407343709295665106703621643087 Recovererd from r2: 407343709295665106703621643087 Library: x1 ^ x2 mod m: 599840174322403511400105423249 231us Montgomery: x1 ^ x2 mod m: 599840174322403511400105423249 2.316ms