Montgomery reduction

Implement the Montgomery reduction algorithm, as explained in "Handbook of Applied Cryptography, Section 14.3.2, page 600. Montgomery reduction calculates $TR^{-1}\mathrm {mod} m$ , without having to divide by $m$ .

Let $M$ be a positive integer, and $R$ and $T$ integers such that $R>m$ , $\mathrm {gcd} (m,R)=1$ , and $0\leq T<mR$ .
$R$ is usually chosen as $b^{n}$ , where $b$ = base (radix) in which the numbers in the calculation as represented in (so $b=10$ in ‘normal’ paper arithmetic, $b=2$ for computer implementations) and $n$ = number of digits in base $m$
The numbers $m$ ( $n$ digits long), $T$ ( $2n$ digits long), $R$ , $b$ , $n$ are known entities, a number $m'$ (often represented as m_dash in code) = $-m^{-1}\mathrm {mod} b$ 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:
   u_i ← a_i* m' mod b      // a_i is the ith digit of A, u_i is a single digit number in radix b
   A ← A + ui*m*bⁱ
A ← A/bⁿ
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