# Montgomery reduction

Montgomery reduction is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Implement the Montgomery reduction algorithm, as explained in "Handbook of Applied Cryptography, Section 14.3.2, page 600. Montgomery reduction calculates TR − 1modm, without having to divide by m.

• Let M be a positive integer, and R and T integers such that R > m, gcd(m,R) = 1, and $0 \le T < mR$.
• R is usually chosen as bn, 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 − 1modb 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:
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++

`#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;}`

## 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} // constructorfunc 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 algorithmfunc (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:    fmt.Println("\nMontgomery computation of x1 ^ x2 mod m:")    // 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))    }    fmt.Println(mr.reduce(prod))     // show library-based equivalent computation as a check    fmt.Println("\nLibrary-based computation of x1 ^ x2 mod m:")    fmt.Println(new(big.Int).Exp(x1, x2, m))}`
Output:
```b:  2
n:  100
r:  1267650600228229401496703205376
m:  750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255

Original x1:        540019781128412936473322405310
Recovererd from r1: 540019781128412936473322405310
Original x2:        515692107665463680305819378593
Recovererd from r2: 515692107665463680305819378593

Montgomery computation of x1 ^ x2 mod m:
151232511393500655853002423778

Library based computation of x1 ^ x2 mod m:
151232511393500655853002423778
```

## Tcl

 This example is under development. It was marked thus on 25/06/2012. Please help complete the example.
`package require Tcl 8.5 proc montgomeryReduction {m mDash T n {b 2}} {    set A \$T    for {set i 0} {\$i < \$n} {incr i} {	# Could be simplified for cases b==2 and b==10	for {set j 0;set a \$A} {\$j < \$i} {incr j} {	    set a [expr {\$a / \$b}]	}	set ui [expr {(\$a % \$b) * \$mDash % \$b}]	incr A [expr {\$ui * \$m * \$b**\$i}]    }    set A [expr {\$A / (\$b ** \$n)}]    return [expr {\$A >= \$m ? \$A - \$m : \$A}]}`