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 the Handbook of Applied Cryptography 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. The said algorithm can be found at [1] at page 600 (page 11 of pdf file)

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 + u_i*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