Rosetta Code

Montgomery reduction: Difference between revisions

Page Discussion

Montgomery reduction (view source)

Revision as of 10:50, 4 January 2024

69,825 bytes added ,  4 months ago
m
→‎{{header|Wren}}: Minor tidy
m (correction)
m (→‎{{header|Wren}}: Minor tidy)
 
(44 intermediate revisions by 18 users not shown)
Line 1:
{{clarify task}}{{draft task}}
Implement the [[wp:Montgomery reduction|Montgomery reduction]] algorithm, as explained in "Handbook of Applied Cryptography, Section 14.3.2, page 600. Montgomery reduction calculates <math>T R^{-1} \mathrm{mod} m</math>, without having to divide by <math>m</math>.
* Let <math>M</math> be a positive integer, and <math>R</math> and <math>T</math> integers such that <math>R > m</math>, <math>\mathrm{gcd}(m, R) = 1</math>, and <math>0 \le T < mR</math>.
Line 16:
A ← A - m
Return (A)
 
=={{header|11l}}==
{{trans|Python}}
 
<syntaxhighlight lang="11l">T Montgomery
BigInt m
Int n
BigInt rrm
 
F (m)
.m = m
.n = bits:length(m)
.rrm = (BigInt(2) ^ (.n * 2)) % m
 
F reduce(t)
V a = t
L 0 .< .n
I (a % 2) == 1
a += .m
a I/= 2
I a >= .m
a -= .m
R a
 
V m = BigInt(‘750791094644726559640638407699’)
V x1 = BigInt(‘540019781128412936473322405310’)
V x2 = BigInt(‘515692107665463680305819378593’)
 
V mont = Montgomery(m)
V t1 = x1 * mont.rrm
V t2 = x2 * mont.rrm
 
V r1 = mont.reduce(t1)
V r2 = mont.reduce(t2)
V r = BigInt(2) ^ mont.n
 
print(‘b : 2’)
print(‘n : ’mont.n)
print(‘r : ’r)
print(‘m : ’mont.m)
print(‘t1: ’t1)
print(‘t2: ’t2)
print(‘r1: ’r1)
print(‘r2: ’r2)
print()
print(‘Original x1 : ’x1)
print(‘Recovered from r1 : ’mont.reduce(r1))
print(‘Original x2 : ’x2)
print(‘Recovered from r2 : ’mont.reduce(r2))
 
print("\nMontgomery computation of x1 ^ x2 mod m:")
V prod = mont.reduce(mont.rrm)
V base = mont.reduce(x1 * mont.rrm)
V ex = x2
L bits:length(ex) > 0
I (ex % 2) == 1
prod = mont.reduce(prod * base)
ex I/= 2
base = mont.reduce(base * base)
print(mont.reduce(prod))
print("\nAlternate computation of x1 ^ x2 mod m:")
print(pow(x1, x2, m))</syntaxhighlight>
 
{{out}}
<pre>
b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m:
151232511393500655853002423778
 
Alternate computation of x1 ^ x2 mod m:
151232511393500655853002423778
</pre>
 
=={{header|C}}==
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
 
typedef unsigned long long ulong;
 
const int BASE = 2;
 
struct Montgomery {
ulong m;
ulong rrm;
int n;
};
 
int bitLength(ulong v) {
int result = 0;
while (v > 0) {
v >>= 1;
result++;
}
return result;
}
 
ulong modPow(ulong b, ulong e, ulong n) {
ulong result = 1;
 
if (n == 1) {
return 0;
}
 
b = b % n;
while (e > 0) {
if (e % 2 == 1) {
result = (result * b) % n;
}
e >>= 1;
b = (b * b) % n;
}
 
return result;
}
 
struct Montgomery makeMontgomery(ulong m) {
struct Montgomery result;
 
if (m == 0 || (m & 1) == 0) {
fprintf(stderr, "m must be greater than zero and odd");
exit(1);
}
 
result.m = m;
result.n = bitLength(m);
result.rrm = (1ULL << (result.n * 2)) % m;
 
return result;
}
 
ulong reduce(struct Montgomery mont, ulong t) {
ulong a = t;
int i;
 
for (i = 0; i < mont.n; i++) {
if ((a & 1) == 1) {
a += mont.m;
}
a = a >> 1;
}
 
if (a >= mont.m) {
a -= mont.m;
}
return a;
}
 
void check(ulong x1, ulong x2, ulong m) {
struct Montgomery mont = makeMontgomery(m);
ulong t1 = x1 * mont.rrm;
ulong t2 = x2 * mont.rrm;
 
ulong r1 = reduce(mont, t1);
ulong r2 = reduce(mont, t2);
ulong r = 1ULL << mont.n;
 
printf("b : %d\n", BASE);
printf("n : %d\n", mont.n);
printf("r : %llu\n", r);
printf("m : %llu\n", mont.m);
printf("t1: %llu\n", t1);
printf("t2: %llu\n", t2);
printf("r1: %llu\n", r1);
printf("r2: %llu\n", r2);
printf("\n");
printf("Original x1 : %llu\n", x1);
printf("Recovered from r1 : %llu\n", reduce(mont, r1));
printf("Original x2 : %llu\n", x2);
printf("Recovered from r2 : %llu\n", reduce(mont, r2));
 
printf("\nMontgomery computation of x1 ^ x2 mod m :\n");
ulong prod = reduce(mont, mont.rrm);
ulong base = reduce(mont, x1 * mont.rrm);
ulong exp = x2;
while (bitLength(exp) > 0) {
if ((exp & 1) == 1) {
prod = reduce(mont, prod * base);
}
exp >>= 1;
base = reduce(mont, base * base);
}
printf("%llu\n", reduce(mont, prod));
printf("\nAlternate computation of x1 ^ x2 mod m :\n");
printf("%llu\n", modPow(x1, x2, m));
}
 
int main() {
check(41, 11, 1549);
 
return 0;
}</syntaxhighlight>
{{out}}
<pre>b : 2
n : 11
r : 2048
m : 1549
t1: 47601
t2: 12771
r1: 322
r2: 842
 
Original x1 : 41
Recovered from r1 : 41
Original x2 : 11
Recovered from r2 : 11
 
Montgomery computation of x1 ^ x2 mod m :
678
 
Alternate computation of x1 ^ x2 mod m :
678</pre>
 
=={{header|C sharp|C#}}==
{{trans|D}}
<syntaxhighlight lang="csharp">using System;
using System.Numerics;
 
namespace MontgomeryReduction {
public static class Helper {
public static int BitLength(this BigInteger v) {
if (v < 0) {
v *= -1;
}
 
int result = 0;
while (v > 0) {
v >>= 1;
result++;
}
 
return result;
}
}
 
struct Montgomery {
public static readonly int BASE = 2;
 
public BigInteger m;
public BigInteger rrm;
public int n;
 
public Montgomery(BigInteger m) {
if (m < 0 || m.IsEven) throw new ArgumentException();
 
this.m = m;
n = m.BitLength();
rrm = (BigInteger.One << (n * 2)) % m;
}
 
public BigInteger Reduce(BigInteger t) {
var a = t;
 
for (int i = 0; i < n; i++) {
if (!a.IsEven) a += m;
a = a >> 1;
}
if (a >= m) a -= m;
return a;
}
}
 
class Program {
static void Main(string[] args) {
var m = BigInteger.Parse("750791094644726559640638407699");
var x1 = BigInteger.Parse("540019781128412936473322405310");
var x2 = BigInteger.Parse("515692107665463680305819378593");
 
var mont = new Montgomery(m);
var t1 = x1 * mont.rrm;
var t2 = x2 * mont.rrm;
 
var r1 = mont.Reduce(t1);
var r2 = mont.Reduce(t2);
var r = BigInteger.One << mont.n;
 
Console.WriteLine("b : {0}", Montgomery.BASE);
Console.WriteLine("n : {0}", mont.n);
Console.WriteLine("r : {0}", r);
Console.WriteLine("m : {0}", mont.m);
Console.WriteLine("t1: {0}", t1);
Console.WriteLine("t2: {0}", t2);
Console.WriteLine("r1: {0}", r1);
Console.WriteLine("r2: {0}", r2);
Console.WriteLine();
Console.WriteLine("Original x1 : {0}", x1);
Console.WriteLine("Recovered from r1 : {0}", mont.Reduce(r1));
Console.WriteLine("Original x2 : {0}", x2);
Console.WriteLine("Recovered from r2 : {0}", mont.Reduce(r2));
 
Console.WriteLine();
Console.WriteLine("Montgomery computation of x1 ^ x2 mod m :");
var prod = mont.Reduce(mont.rrm);
var @base = mont.Reduce(x1 * mont.rrm);
var exp = x2;
while (exp.BitLength() > 0) {
if (!exp.IsEven) prod = mont.Reduce(prod * @base);
exp >>= 1;
@base = mont.Reduce(@base * @base);
}
Console.WriteLine(mont.Reduce(prod));
Console.WriteLine();
Console.WriteLine("Alternate computation of x1 ^ x2 mod m :");
Console.WriteLine(BigInteger.ModPow(x1, x2, m));
}
}
}</syntaxhighlight>
{{out}}
<pre>b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :
151232511393500655853002423778
 
Alternate computation of x1 ^ x2 mod m :
151232511393500655853002423778</pre>
 
=={{header|C++}}==
<langsyntaxhighlight lang="cpp">#include<iostream>
#include<conio.h>
using namespace std;
Line 103 ⟶ 442:
cout<<"Montgomery domain representation = "<<e;
return 0;
}</langsyntaxhighlight>
 
=={{header|D}}==
{{trans|Kotlin}}
<syntaxhighlight lang="d">import std.bigint;
import std.stdio;
 
int bitLength(BigInt v) {
if (v < 0) {
v *= -1;
}
 
int result = 0;
while (v > 0) {
v >>= 1;
result++;
}
 
return result;
}
 
/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
BigInt modPow(BigInt b, BigInt e, BigInt n) {
if (n == 1) return BigInt(0);
BigInt result = 1;
b = b % n;
while (e > 0) {
if (e % 2 == 1) {
result = (result * b) % n;
}
e >>= 1;
b = (b*b) % n;
}
return result;
}
 
struct Montgomery {
BigInt m;
int n;
BigInt rrm;
 
this(BigInt m) in {
assert(m > 0 && (m & 1) != 0); // must be positive and odd
} body {
this.m = m;
n = m.bitLength();
rrm = (BigInt(1) << (n * 2)) % m;
}
 
BigInt reduce(BigInt t) {
auto a = t;
 
foreach(i; 0..n) {
if ((a & 1) == 1) a += m;
a = a >> 1;
}
if (a >= m) a -= m;
return a;
}
 
enum BASE = 2;
}
 
void main() {
auto m = BigInt("750791094644726559640638407699");
auto x1 = BigInt("540019781128412936473322405310");
auto x2 = BigInt("515692107665463680305819378593");
 
auto mont = Montgomery(m);
auto t1 = x1 * mont.rrm;
auto t2 = x2 * mont.rrm;
 
auto r1 = mont.reduce(t1);
auto r2 = mont.reduce(t2);
auto r = BigInt(1) << mont.n;
 
writeln("b : ", Montgomery.BASE);
writeln("n : ", mont.n);
writeln("r : ", r);
writeln("m : ", mont.m);
writeln("t1: ", t1);
writeln("t2: ", t2);
writeln("r1: ", r1);
writeln("r2: ", r2);
writeln();
writeln("Original x1 : ", x1);
writeln("Recovered from r1 : ", mont.reduce(r1));
writeln("Original x2 : ", x2);
writeln("Recovered from r2 : ", mont.reduce(r2));
 
writeln("\nMontgomery computation of x1 ^ x2 mod m :");
auto prod = mont.reduce(mont.rrm);
auto base = mont.reduce(x1 * mont.rrm);
auto exp = x2;
while (exp.bitLength() > 0) {
if ((exp & 1) == 1) prod = mont.reduce(prod * base);
exp >>= 1;
base = mont.reduce(base * base);
}
writeln(mont.reduce(prod));
writeln("\nAlternate computation of x1 ^ x2 mod m :");
writeln(x1.modPow(x2, m));
}</syntaxhighlight>
{{out}}
<pre>b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :
151232511393500655853002423778
 
Alternate computation of x1 ^ x2 mod m :
151232511393500655853002423778</pre>
 
=={{header|Factor}}==
{{trans|Sidef}}
{{works with|Factor|0.99 2020-08-14}}
<syntaxhighlight lang="factor">USING: io kernel locals math math.bitwise math.functions
prettyprint ;
 
: montgomery-reduce ( m a -- n )
over bit-length [ dup odd? [ over + ] when 2/ ] times
swap mod ;
 
CONSTANT: m 750791094644726559640638407699
CONSTANT: t1 323165824550862327179367294465482435542970161392400401329100
 
CONSTANT: r1 440160025148131680164261562101
CONSTANT: r2 435362628198191204145287283255
 
CONSTANT: x1 540019781128412936473322405310
CONSTANT: x2 515692107665463680305819378593
 
"Original x1: " write x1 .
"Recovered from r1: " write m r1 montgomery-reduce .
"Original x2: " write x2 .
"Recovered from r2: " write m r2 montgomery-reduce .
 
nl "Montgomery computation of x1^x2 mod m: " write
 
[let
m t1 x1 / montgomery-reduce :> prod!
m t1 montgomery-reduce :> base!
x2 :> exponent!
 
[ exponent zero? ] [
exponent odd?
[ m prod base * montgomery-reduce prod! ] when
m base sq montgomery-reduce base! exponent 2/ exponent!
] until
 
m prod montgomery-reduce .
"Library-based computation of x1^x2 mod m: " write
x1 x2 m ^mod .
]</syntaxhighlight>
{{out}}
<pre>
Original x1: 540019781128412936473322405310
Recovered from r1: 540019781128412936473322405310
Original x2: 515692107665463680305819378593
Recovered from r2: 515692107665463680305819378593
 
Montgomery computation of x1^x2 mod m: 151232511393500655853002423778
Library-based computation of x1^x2 mod m: 151232511393500655853002423778
</pre>
 
=={{header|Go}}==
<langsyntaxhighlight 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 {
Line 145 ⟶ 659:
}
return a
}
 
// example use:
func main() {
Line 194 ⟶ 708:
 
// and demonstrate a use:
fmt.Println("\nMontgomery computation of x1 ^ x2 mod m:")
show("\nLibrary:", func() *big.Int {
// this is the modular exponentiation algorithm, except we call
return new(big.Int).Exp(x1, x2, m)
// mont.reduce instead of using a mod function.
})
prod := mr.reduce(mr.r2) // 1
show("\nMontgomery:", func() *big.Int {
base := mr.reduce(t1.Mul(x1, mr.r2)) // x1^1
// this is the modular exponentiation algorithm, except we call
exp := new(big.Int).Set(x2) // not reduced
// mont.reduce instead of using a mod function.
for prod := mrexp.reduceBitLen(mr.r2) > 0 // 1{
baseif := mrexp.reduceBit(t1.Mul(x1, mr.r2)0) //== x1^1 {
prod = mr.reduce(prod.Mul(prod, base))
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 mrexp.reduceRsh(prodexp, 1)
base = mr.reduce(base.Mul(base, base))
})
} }
fmt.Println(mr.reduce(prod))
 
// show library-based equivalent computation as a check
func show(heading string, f func() *big.Int) {
fmt.Println(heading"\nLibrary-based computation of x1 ^ x2 mod m:")
fmt.Println(new(big.Int).Exp(x1, x2, m))
t0 := time.Now()
}</syntaxhighlight>
fmt.Println("x1 ^ x2 mod m:", f())
{{out}}
fmt.Println(time.Now().Sub(t0))
}</lang>
Output:
<pre>
b: 2
n: 100
r: 1267650600228229401496703205376
m: 750791094644726559640638407699
m: 1191151171693032142151966564621
t1: 323165824550862327179367294465482435542970161392400401329100
t1: 138602318824179275477121611740471794618999274340657947731554
t2: 308607334419945011411837686695175944083084270671482464168730
t2: 397645077922552596057001924720561733079744309909810064024787
r1: 440160025148131680164261562101
r1: 1073769066977456952449715239458
r2: 435362628198191204145287283255
r2: 460620928979319347925360938242
 
Original x1: 141982853055250888109454150154540019781128412936473322405310
Recovererd from r1: 141982853055250888109454150154540019781128412936473322405310
Original x2: 407343709295665106703621643087515692107665463680305819378593
Recovererd from r2: 407343709295665106703621643087515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m:
Library:
151232511393500655853002423778
x1 ^ x2 mod m: 599840174322403511400105423249
 
231us
Library based computation of x1 ^ x2 mod m:
151232511393500655853002423778
</pre>
 
=={{header|Java}}==
{{trans|Kotlin}}
<syntaxhighlight lang="java">import java.math.BigInteger;
 
public class MontgomeryReduction {
private static final BigInteger ZERO = BigInteger.ZERO;
private static final BigInteger ONE = BigInteger.ONE;
private static final BigInteger TWO = BigInteger.valueOf(2);
 
public static class Montgomery {
public static final int BASE = 2;
 
BigInteger m;
BigInteger rrm;
int n;
 
public Montgomery(BigInteger m) {
if (m.compareTo(BigInteger.ZERO) <= 0 || !m.testBit(0)) {
throw new IllegalArgumentException();
}
this.m = m;
this.n = m.bitLength();
this.rrm = ONE.shiftLeft(n * 2).mod(m);
}
 
public BigInteger reduce(BigInteger t) {
BigInteger a = t;
for (int i = 0; i < n; i++) {
if (a.testBit(0)) a = a.add(this.m);
a = a.shiftRight(1);
}
if (a.compareTo(m) >= 0) a = a.subtract(this.m);
return a;
}
}
 
public static void main(String[] args) {
BigInteger m = new BigInteger("750791094644726559640638407699");
BigInteger x1 = new BigInteger("540019781128412936473322405310");
BigInteger x2 = new BigInteger("515692107665463680305819378593");
 
Montgomery mont = new Montgomery(m);
BigInteger t1 = x1.multiply(mont.rrm);
BigInteger t2 = x2.multiply(mont.rrm);
 
BigInteger r1 = mont.reduce(t1);
BigInteger r2 = mont.reduce(t2);
BigInteger r = ONE.shiftLeft(mont.n);
 
System.out.printf("b : %s\n", Montgomery.BASE);
System.out.printf("n : %s\n", mont.n);
System.out.printf("r : %s\n", r);
System.out.printf("m : %s\n", mont.m);
System.out.printf("t1: %s\n", t1);
System.out.printf("t2: %s\n", t2);
System.out.printf("r1: %s\n", r1);
System.out.printf("r2: %s\n", r2);
System.out.println();
System.out.printf("Original x1 : %s\n", x1);
System.out.printf("Recovered from r1 : %s\n", mont.reduce(r1));
System.out.printf("Original x2 : %s\n", x2);
System.out.printf("Recovered from r2 : %s\n", mont.reduce(r2));
 
System.out.println();
System.out.println("Montgomery computation of x1 ^ x2 mod m :");
BigInteger prod = mont.reduce(mont.rrm);
BigInteger base = mont.reduce(x1.multiply(mont.rrm));
BigInteger exp = x2;
while (exp.bitLength()>0) {
if (exp.testBit(0)) prod=mont.reduce(prod.multiply(base));
exp = exp.shiftRight(1);
base = mont.reduce(base.multiply(base));
}
System.out.println(mont.reduce(prod));
 
System.out.println();
System.out.println("Library-based computation of x1 ^ x2 mod m :");
System.out.println(x1.modPow(x2, m));
}
}</syntaxhighlight>
{{out}}
<pre>b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :
151232511393500655853002423778
 
Library-based computation of x1 ^ x2 mod m :
151232511393500655853002423778</pre>
 
=={{header|Julia}}==
{{trans|Python}}
<syntaxhighlight lang="julia">""" base 2 type Montgomery numbers """
struct Montgomery2
m::BigInt
n::Int64
rrm::BigInt
end
 
function Montgomery2(x::BigInt)
bitlen = length(string(x, base=2))
r = (x == 0) ? 0 : (BigInt(1) << (bitlen * 2)) % x
Montgomery2(x, bitlen, r)
end
Montgomery2(n) = Montgomery2(BigInt(n))
 
function reduce(mm::Montgomery2, t)
a = BigInt(t)
for i in 1:mm.n
if isodd(a)
a += mm.m
end
a >>= 1
end
return a >= mm.m ? a - mm.m : a
end
 
BASE(::Montgomery2) = 2
 
const mmm = Montgomery2(20)
 
function testmontgomery2()
m = big"750791094644726559640638407699"
x1 = big"540019781128412936473322405310"
x2 = big"515692107665463680305819378593"
 
mont = Montgomery2(m)
t1 = x1 * mont.rrm
t2 = x2 * mont.rrm
r1 = reduce(mont, t1)
r2 = reduce(mont, t2)
r = 1 << mont.n
println("b : ", BASE(mont))
println("n : ", mont.n)
println("r : ", r)
println("m : ", mont.m)
println("t1: ", t1)
println("t2: ", t2)
println("r1: ", r1)
println("r2: ", r2)
println()
println("Original x1 :", x1)
println("Recovered from r1 :", reduce(mont, r1))
println("Original x2 :", x2)
println("Recovered from r2 :", reduce(mont, r2))
println("\nMontgomery computation of x1 ^ x2 mod m:")
prod = reduce(mont, mont.rrm)
base = reduce(mont, x1 * mont.rrm)
pow = x2
while pow > 0
if isodd(pow)
prod = reduce(mont, prod * base)
end
pow >>= 1
base = reduce(mont, base * base)
end
println(reduce(mont, prod))
println("\nAlternate computation of x1 ^ x2 mod m :")
println(powermod(x1, x2, m))
end
 
testmontgomery2()
</syntaxhighlight>{{out}}
<pre>
b : 2
n : 100
r : 0
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 :540019781128412936473322405310
Recovered from r1 :540019781128412936473322405310
Original x2 :515692107665463680305819378593
Recovered from r2 :515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m:
151232511393500655853002423778
 
Alternate computation of x1 ^ x2 mod m :
151232511393500655853002423778
</pre>
 
=={{header|Kotlin}}==
{{trans|Go}}
<syntaxhighlight lang="scala">// version 1.1.3
 
import java.math.BigInteger
 
val bigZero = BigInteger.ZERO
val bigOne = BigInteger.ONE
val bigTwo = BigInteger.valueOf(2L)
 
class Montgomery(val m: BigInteger) {
val n: Int
val rrm: BigInteger
 
init {
require(m > bigZero && m.testBit(0)) // must be positive and odd
n = m.bitLength()
rrm = bigOne.shiftLeft(n * 2).mod(m)
}
 
fun reduce(t: BigInteger): BigInteger {
var a = t
for (i in 0 until n) {
if (a.testBit(0)) a += m
a = a.shiftRight(1)
}
if (a >= m) a -= m
return a
}
 
companion object {
const val BASE = 2
}
}
 
fun main(args: Array<String>) {
val m = BigInteger("750791094644726559640638407699")
val x1 = BigInteger("540019781128412936473322405310")
val x2 = BigInteger("515692107665463680305819378593")
 
val mont = Montgomery(m)
val t1 = x1 * mont.rrm
val t2 = x2 * mont.rrm
 
val r1 = mont.reduce(t1)
val r2 = mont.reduce(t2)
val r = bigOne.shiftLeft(mont.n)
 
println("b : ${Montgomery.BASE}")
println("n : ${mont.n}")
println("r : $r")
println("m : ${mont.m}")
println("t1: $t1")
println("t2: $t2")
println("r1: $r1")
println("r2: $r2")
println()
println("Original x1 : $x1")
println("Recovered from r1 : ${mont.reduce(r1)}")
println("Original x2 : $x2")
println("Recovered from r2 : ${mont.reduce(r2)}")
 
println("\nMontgomery computation of x1 ^ x2 mod m :")
var prod = mont.reduce(mont.rrm)
var base = mont.reduce(x1 * mont.rrm)
var exp = x2
while (exp.bitLength() > 0) {
if (exp.testBit(0)) prod = mont.reduce(prod * base)
exp = exp.shiftRight(1)
base = mont.reduce(base * base)
}
println(mont.reduce(prod))
println("\nLibrary-based computation of x1 ^ x2 mod m :")
println(x1.modPow(x2, m))
}</syntaxhighlight>
 
{{out}}
<pre>
b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :
151232511393500655853002423778
 
Library-based computation of x1 ^ x2 mod m :
151232511393500655853002423778
</pre>
 
=={{header|Nim}}==
{{trans|D}}
{{libheader|bignum}}
<syntaxhighlight lang="nim">import bignum
 
# Missing functions in "bignum".
 
template isOdd(val: Int): bool =
## Needed as bignum "odd" function crashes.
(val and 1) != 0
 
func exp(x, y, m: Int): Int =
## Missing "exp" function in "bignum".
if m == 1: return newInt(0)
result = newInt(1)
var x = x mod m
var y = y
while y > 0:
if y.isOdd:
result = (result * x) mod m
y = y shr 1
x = (x * x) mod m
 
 
type Montgomery = object
m: Int # Modulus; must be odd.
n: int # m.bitLen().
rrm: Int # (1<<2n) mod m.
 
const Base = 2
 
 
func initMontgomery(m: Int): Montgomery =
## Initialize a Mongtgomery object.
doAssert m > 0 and m.isOdd, "argument must be positive and odd."
result.m = m
result.n = m.bitLen
result.rrm = newInt(1) shl culong(result.n * 2) mod m
 
 
func reduce(mont: Montgomery; t: Int): Int =
## Montgomery reduction algorithm.
result = t
for i in 0..<mont.n:
if result.isOdd: inc result, mont.m
result = result shr 1
if result >= mont.m: dec result, mont.m
 
 
when isMainModule:
 
let
m = newInt("750791094644726559640638407699")
x1 = newInt("540019781128412936473322405310")
x2 = newInt("515692107665463680305819378593")
 
mont = initMontgomery(m)
t1 = x1 * mont.rrm
t2 = x2 * mont.rrm
 
r1 = mont.reduce(t1)
r2 = mont.reduce(t2)
r = newInt(1) shl culong(mont.n)
 
echo "b: ", Base
echo "n: ", mont.n
echo "r: ", r
echo "m: ", mont.m
echo "t1: ", t1
echo "t2: ", t2
echo "r1: ", r1
echo "r2: ", r2
echo()
echo "Original x1: ", x1
echo "Recovered from r1: ", mont.reduce(r1)
echo "Original x2: ", x2
echo "Recovered from r2: ", mont.reduce(r2)
 
echo "\nMontgomery computation of x1^x2 mod m:"
var
prod = mont.reduce(mont.rrm)
base = mont.reduce(x1 * mont.rrm)
e = x2
while e > 0:
if e.isOdd: prod = mont.reduce(prod * base)
e = e shr 1
base = mont.reduce(base * base)
echo mont.reduce(prod)
echo "\nAlternate computation of x1^x2 mod m:"
echo x1.exp(x2, m)</syntaxhighlight>
 
{{out}}
<pre>b: 2
n: 100
r: 1267650600228229401496703205376
m: 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1: 540019781128412936473322405310
Recovered from r1: 540019781128412936473322405310
Original x2: 515692107665463680305819378593
Recovered from r2: 515692107665463680305819378593
 
Montgomery computation of x1^x2 mod m:
151232511393500655853002423778
 
Alternate computation of x1^x2 mod m:
151232511393500655853002423778</pre>
 
=={{header|Perl}}==
{{trans|Raku}}
{{libheader|ntheory}}
<syntaxhighlight lang="perl">use bigint;
use ntheory qw(powmod);
 
sub msb {
my ($n, $base) = (shift, 0);
$base++ while $n >>= 1;
$base;
}
 
sub montgomery_reduce {
my($m, $a) = @_;
for (0 .. msb($m)) {
$a += $m if $a & 1;
$a >>= 1
}
$a % $m
}
 
my $m = 750791094644726559640638407699;
my $t1 = 323165824550862327179367294465482435542970161392400401329100;
 
my $r1 = 440160025148131680164261562101;
my $r2 = 435362628198191204145287283255;
 
my $x1 = 540019781128412936473322405310;
my $x2 = 515692107665463680305819378593;
 
printf "Original x1: %s\n", $x1;
printf "Recovered from r1: %s\n", montgomery_reduce($m, $r1);
printf "Original x2: %s\n", $x2;
printf "Recovered from r2: %s\n", montgomery_reduce($m, $r2);
 
print "\nMontgomery computation x1**x2 mod m: ";
my $prod = montgomery_reduce($m, $t1/$x1);
my $base = montgomery_reduce($m, $t1);
 
for (my $exponent = $x2; $exponent >= 0; $exponent >>= 1) {
$prod = montgomery_reduce($m, $prod * $base) if $exponent & 1;
$base = montgomery_reduce($m, $base * $base);
last if $exponent == 0;
}
 
print montgomery_reduce($m, $prod) . "\n";
printf "Built-in op computation x1**x2 mod m: %s\n", powmod($x1, $x2, $m);</syntaxhighlight>
{{out}}
<pre>Original x1: 540019781128412936473322405310
Recovered from r1: 540019781128412936473322405310
Original x2: 515692107665463680305819378593
Recovered from r2: 515692107665463680305819378593
 
Montgomery computation x1**x2 mod m: 151232511393500655853002423778
Built-in op computation x1**x2 mod m: 151232511393500655853002423778</pre>
 
=={{header|Phix}}==
{{trans|D}}
{{libheader|Phix/mpfr}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">BASE</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">BITLEN</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">MODULUS</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">RRM</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">mpz</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">BITLEN</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">MODULUS</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_odd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_fdiv_q_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)>=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">mpz_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">r</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">Montgomery</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_sign</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)=-</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"must be positive"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #7060A8;">mpz_odd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"must be odd"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_sizeinbase</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">rrm</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_powm_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rrm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rrm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">*</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- BASE</span>
<span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- BITLEN</span>
<span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- MODULUS</span>
<span style="color: #000000;">rrm</span> <span style="color: #000080;font-style:italic;">-- 1&lt;&lt;(n*2) % m</span>
<span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"750791094644726559640638407699"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">x1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"540019781128412936473322405310"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">x2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"515692107665463680305819378593"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">t2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">mont</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">Montgomery</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">RRM</span><span style="color: #0000FF;">])</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">RRM</span><span style="color: #0000FF;">])</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">r1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t1</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">r2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">mpz_ui_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">BITLEN</span><span style="color: #0000FF;">])</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"b : %d\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">BASE</span><span style="color: #0000FF;">]})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"n : %d\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">BITLEN</span><span style="color: #0000FF;">]})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"r : %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"m : %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">MODULUS</span><span style="color: #0000FF;">])})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"t1: %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t1</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"t2: %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t2</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"r1: %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r1</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"r2: %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Original x1 : %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Recovered from r1 : %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r1</span><span style="color: #0000FF;">))})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Original x2 : %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Recovered from r2 : %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r2</span><span style="color: #0000FF;">))})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\nMontgomery computation of x1 ^ x2 mod m :"</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">prod</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">RRM</span><span style="color: #0000FF;">])</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">[</span><span style="color: #000000;">RRM</span><span style="color: #0000FF;">])</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">base</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">expn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">while</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">expn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_odd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">expn</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">prod</span><span style="color: #0000FF;">,</span><span style="color: #000000;">prod</span><span style="color: #0000FF;">,</span><span style="color: #000000;">base</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">prod</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">prod</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_fdiv_q_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">expn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">expn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">base</span><span style="color: #0000FF;">,</span><span style="color: #000000;">base</span><span style="color: #0000FF;">,</span><span style="color: #000000;">base</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">base</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">base</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mont</span><span style="color: #0000FF;">,</span><span style="color: #000000;">prod</span><span style="color: #0000FF;">))})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" alternate computation of x1 ^ x2 mod m :"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_powm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :151232511393500655853002423778
alternate computation of x1 ^ x2 mod m :151232511393500655853002423778
</pre>
 
=={{header|PicoLisp}}==
<syntaxhighlight lang="picolisp">(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y))) M)
(setq X (% (* X X) N)) ) ) )
(de rrm (M)
(% (>> (- (* 2 Mbins)) 1) M) )
(de reduce (A)
(do Mbins
(and (bit? 1 A) (inc 'A M))
(setq A (>> 1 A)) )
(and (>= A M) (dec 'A M))
A )
(let
(M 750791094644726559640638407699
Mbins (length (bin M))
RRM (rrm M)
X1 540019781128412936473322405310
X2 515692107665463680305819378593
T1 (* X1 RRM)
T2 (* X2 RRM)
R1 (reduce T1)
R2 (reduce T2)
R (>> (- Mbins) 1)
Prod (reduce RRM)
Base (reduce (* X1 RRM))
Exp X2 )
(println 'b ': 2)
(println 'n ': Mbins)
(println 'r ': R)
(println 'm ': M)
(println 't1 ': T1)
(println 't2 ': T2)
(println 'r1 ': R1)
(println 'r2 ': R2)
(prinl)
(prinl "Original x1 : " X1)
(prinl "Recovered from r1 : " (reduce R1))
(prinl "Original x2 : " X2)
(prinl "Recovered from r2 : " (reduce R2))
(prinl)
(prin "Montgomery computation of x1 \^ x2 mod m : ")
(while (gt0 Exp)
(and
(bit? 1 Exp)
(setq Prod (reduce (* Prod Base))) )
(setq
Exp (>> 1 Exp)
Base (reduce (* Base Base)) ) )
(prinl (reduce Prod))
(prinl "Montgomery computation of x1 \^ x2 mod m : " (**Mod X1 X2 M)) )</syntaxhighlight>
{{out}}
<pre>b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1 : 323165824550862327179367294465482435542970161392400401329100
t2 : 308607334419945011411837686695175944083084270671482464168730
r1 : 440160025148131680164261562101
r2 : 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m : 151232511393500655853002423778
Montgomery computation of x1 ^ x2 mod m : 151232511393500655853002423778</pre>
 
=={{header|Python}}==
{{todo|python|Update the output}}
{{trans|D}}
<syntaxhighlight lang="python">
class Montgomery:
BASE = 2
 
def __init__(self, m):
self.m = m
self.n = m.bit_length()
self.rrm = (1 << (self.n * 2)) % m
 
def reduce(self, t):
a = t
for i in xrange(self.n):
if (a & 1) == 1:
a = a + self.m
a = a >> 1
if a >= self.m:
a = a - self.m
return a
 
# Main
m = 750791094644726559640638407699
x1 = 540019781128412936473322405310
x2 = 515692107665463680305819378593
 
mont = Montgomery(m)
t1 = x1 * mont.rrm
t2 = x2 * mont.rrm
 
r1 = mont.reduce(t1)
r2 = mont.reduce(t2)
r = 1 << mont.n
 
print(
f"b: {Montgomery.BASE}\n"
f"n: {mont.n}\n"
f"r: {r}\n"
f"m: {mont.m}\n"
f"t1: {t1}\n"
f"t2: {t2}\n"
f"r1: {r1}\n"
f"r2: {r2}\n"
f"Original x1: {x1}\n"
f"Recovered from r1: {mont.reduce(r1)}\n"
f"Original x2: {x2}\n"
f"Recovered from r2: {mont.reduce(r2)}\n"
)
 
print("Montgomery computation of x1 ^ x2 mod m:")
prod = mont.reduce(mont.rrm)
base = mont.reduce(x1 * mont.rrm)
exp = x2
while exp.bit_length() > 0:
if (exp & 1) == 1:
prod = mont.reduce(prod * base)
exp = exp >> 1
base = mont.reduce(base * base)
print(mont.reduce(prod))
print(f"\nAlternate computation of x1 ^ x2 mod m: {pow(x1, x2, m)}")
</syntaxhighlight>
{{out}}
<pre>
b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m:
151232511393500655853002423778
 
Alternate computation of x1 ^ x2 mod m :
151232511393500655853002423778
</pre>
 
=={{header|Quackery}}==
 
{{trans|Factor}}
 
<code>**mod</code> is defined at [[Modular exponentiation#Quackery]].
 
<syntaxhighlight lang="Quackery"> [ 0 swap [ dup while dip 1+ 1 >> again ] drop ] is bits ( n --> n )
 
[ 1 & ] is odd ( n --> b )
 
[ over bits times [ dup odd if [ over + ] 1 >> ] swap mod ] is monred ( n n --> n )
 
[ 750791094644726559640638407699 ] is m ( --> n )
[ 323165824550862327179367294465482435542970161392400401329100 ] is t1 ( --> n )
 
[ 440160025148131680164261562101 ] is r1 ( --> n )
[ 435362628198191204145287283255 ] is r2 ( --> n )
 
[ 540019781128412936473322405310 ] is x1 ( --> n )
[ 515692107665463680305819378593 ] is x2 ( --> n )
 
[ unrot dip
[ dip dup
over t1 rot / monred temp put
t1 monred ]
[ dup 0 != while
dup odd if
[ over temp take *
m swap monred
temp put ]
dip [ dup * m swap monred ]
1 >> again ]
2drop temp take monred ] is **mon ( n n n --> n )
 
say "Original x1: " x1 echo cr
say "Recovered from r1: " m r1 monred echo cr
cr
say "Original x2: " x2 echo cr
say "Recovered from r2: " m r2 monred echo cr
cr
say "Montgomery computation of x1^x2 mod m: " x1 x2 m **mon echo cr
say "Modular exponentiation of x1^x2 mod m: " x1 x2 m **mod echo cr</syntaxhighlight>
 
{{out}}
 
<pre>Original x1: 540019781128412936473322405310
Recovered from r1: 540019781128412936473322405310
 
Original x2: 515692107665463680305819378593
Recovered from r2: 515692107665463680305819378593
 
Montgomery computation of x1^x2 mod m: 151232511393500655853002423778
Modular exponentiation of x1^x2 mod m: 151232511393500655853002423778</pre>
 
=={{header|Racket}}==
 
<syntaxhighlight lang="racket">#lang typed/racket
(require math/number-theory)
 
(: montgomery-reduce-fn
(-> Positive-Integer Natural [#:m-dash Natural]
(Nonnegative-Integer Natural -> Integer)))
 
(: ith-digit (Integer Nonnegative-Integer Natural -> Nonnegative-Integer))
(define (ith-digit a i b)
(modulo (quotient a (expt b i)) b))
 
(: m-dash (Integer Integer -> Natural))
(define (m-dash m b) ; for if you want to precompute it yourself
(modular-inverse (- m) b))
 
(define ((montgomery-reduce-fn m b #:m-dash (m′ (m-dash m b))) T n)
(define A
(for/fold : Nonnegative-Integer
((A : Nonnegative-Integer T))
((i : Nonnegative-Integer (in-range n)))
(let* ((a_i (ith-digit A i b))
(u_i (modulo (* a_i m′) b)))
(+ A (* u_i m (expt b i))))))
(define A/b^n (quotient A (expt b n)))
(if (>= A/b^n m)
(- A/b^n m)
A/b^n))
 
; ---------------------------------------------------------------------------------------------------
(module+ test
(require typed/rackunit)
(check-equal? (ith-digit 1234 0 10) 4)
(check-equal? (ith-digit 1234 3 10) 1)
;; e.g. ripped off from {{trans|Go}}
(let ((b 2)
(n 100)
(r 1267650600228229401496703205376)
(m 750791094644726559640638407699)
(T1 323165824550862327179367294465482435542970161392400401329100)
(T2 308607334419945011411837686695175944083084270671482464168730)
(R1 440160025148131680164261562101)
(R2 435362628198191204145287283255)
(x1 540019781128412936473322405310)
(x2 515692107665463680305819378593))
(define mr (montgomery-reduce-fn m b))
(check-equal? (mr R1 n) x1)
(check-equal? (mr R2 n) x2)))</syntaxhighlight>
 
Tests, which are courtesy of #Go implementation, all pass.
 
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2018.03}}
{{trans|Sidef}}
 
<syntaxhighlight lang="raku" line>sub montgomery-reduce($m, $a is copy) {
for 0..$m.msb {
$a += $m if $a +& 1;
$a +>= 1
}
$a % $m
}
 
my $m = 750791094644726559640638407699;
my $t1 = 323165824550862327179367294465482435542970161392400401329100;
 
my $r1 = 440160025148131680164261562101;
my $r2 = 435362628198191204145287283255;
 
my $x1 = 540019781128412936473322405310;
my $x2 = 515692107665463680305819378593;
 
say "Original x1: ", $x1;
say "Recovered from r1: ", montgomery-reduce($m, $r1);
say "Original x2: ", $x2;
say "Recovered from r2: ", montgomery-reduce($m, $r2);
 
print "\nMontgomery computation x1**x2 mod m: ";
my $prod = montgomery-reduce($m, $t1/$x1);
my $base = montgomery-reduce($m, $t1);
 
for $x2, {$_ +> 1} ... 0 -> $exponent {
$prod = montgomery-reduce($m, $prod * $base) if $exponent +& 1;
$base = montgomery-reduce($m, $base * $base);
}
 
say montgomery-reduce($m, $prod);
say "Built-in op computation x1**x2 mod m: ", $x1.expmod($x2, $m);</syntaxhighlight>
{{out}}
<pre>Original x1: 540019781128412936473322405310
Recovered from r1: 540019781128412936473322405310
Original x2: 515692107665463680305819378593
Recovered from r2: 515692107665463680305819378593
 
Montgomery computation x1**x2 mod m: 151232511393500655853002423778
Built-in op computation x1**x2 mod m: 151232511393500655853002423778
</pre>
 
=={{header|Sidef}}==
{{trans|zkl}}
<syntaxhighlight lang="ruby">func montgomeryReduce(m, a) {
{
a += m if a.is_odd
a >>= 1
} * m.as_bin.len
 
a % m
}
 
var m = 750791094644726559640638407699
var t1 = 323165824550862327179367294465482435542970161392400401329100
 
var r1 = 440160025148131680164261562101
var r2 = 435362628198191204145287283255
 
var x1 = 540019781128412936473322405310
var x2 = 515692107665463680305819378593
 
say("Original x1: ", x1)
say("Recovererd from r1: ", montgomeryReduce(m, r1))
say("Original x2: ", x2)
say("Recovererd from r2: ", montgomeryReduce(m, r2))
 
print("\nMontgomery computation of x1^x2 mod m: ")
var prod = montgomeryReduce(m, t1/x1)
var base = montgomeryReduce(m, t1)
 
for (var exponent = x2; exponent ; exponent >>= 1) {
prod = montgomeryReduce(m, prod * base) if exponent.is_odd
base = montgomeryReduce(m, base * base)
}
 
say(montgomeryReduce(m, prod))
say("Library-based computation of x1^x2 mod m: ", x1.powmod(x2, m))</syntaxhighlight>
{{out}}
<pre>
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
</pre>
 
=={{header|Tcl}}==
{{in progress|lang=Tcl|day=25|month=06|year=2012}}
<syntaxhighlight lang="tcl">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}]
}</syntaxhighlight>
<!-- Not quite sure how to demonstrate this working; examples above aren't very clear… -->
 
=={{header|Visual Basic .NET}}==
{{trans|C#}}
<syntaxhighlight lang="vbnet">Imports System.Numerics
Imports System.Runtime.CompilerServices
 
Module Module1
 
<Extension()>
Function BitLength(v As BigInteger) As Integer
If v < 0 Then
v *= -1
End If
 
Dim result = 0
While v > 0
v >>= 1
result += 1
End While
Return result
End Function
 
Structure Montgomery
Shared ReadOnly BASE = 2
Dim m As BigInteger
Dim rrm As BigInteger
Dim n As Integer
 
Sub New(m As BigInteger)
If m < 0 OrElse m.IsEven Then
Throw New ArgumentException()
End If
 
Me.m = m
n = m.BitLength
rrm = (BigInteger.One << (n * 2)) Mod m
End Sub
 
Function Reduce(t As BigInteger) As BigInteger
Dim a = t
For i = 1 To n
If Not a.IsEven Then
a += m
End If
a >>= 1
Next
If a >= m Then
a -= m
End If
Return a
End Function
End Structure
 
Sub Main()
Dim m = BigInteger.Parse("750791094644726559640638407699")
Dim x1 = BigInteger.Parse("540019781128412936473322405310")
Dim x2 = BigInteger.Parse("515692107665463680305819378593")
 
Dim mont As New Montgomery(m)
Dim t1 = x1 * mont.rrm
Dim t2 = x2 * mont.rrm
 
Dim r1 = mont.Reduce(t1)
Dim r2 = mont.Reduce(t2)
Dim r = BigInteger.One << mont.n
 
Console.WriteLine("b : {0}", Montgomery.BASE)
Console.WriteLine("n : {0}", mont.n)
Console.WriteLine("r : {0}", r)
Console.WriteLine("m : {0}", mont.m)
Console.WriteLine("t1: {0}", t1)
Console.WriteLine("t2: {0}", t2)
Console.WriteLine("r1: {0}", r1)
Console.WriteLine("r2: {0}", r2)
Console.WriteLine()
Console.WriteLine("Original x1 : {0}", x1)
Console.WriteLine("Recovered from r1 : {0}", mont.Reduce(r1))
Console.WriteLine("Original x2 : {0}", x2)
Console.WriteLine("Recovered from r2 : {0}", mont.Reduce(r2))
 
Console.WriteLine()
Console.WriteLine("Montgomery computation of x1 ^ x2 mod m :")
Dim prod = mont.Reduce(mont.rrm)
Dim base = mont.Reduce(x1 * mont.rrm)
Dim exp = x2
While exp.BitLength > 0
If Not exp.IsEven Then
prod = mont.Reduce(prod * base)
End If
exp >>= 1
base = mont.Reduce(base * base)
End While
Console.WriteLine(mont.Reduce(prod))
Console.WriteLine()
Console.WriteLine("Alternate computation of x1 ^ x2 mod m :")
Console.WriteLine(BigInteger.ModPow(x1, x2, m))
End Sub
 
End Module</syntaxhighlight>
{{out}}
<pre>b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :
151232511393500655853002423778
 
Alternate computation of x1 ^ x2 mod m :
151232511393500655853002423778</pre>
 
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-big}}
<syntaxhighlight lang="wren">import "./big" for BigInt
 
class Montgomery {
static base { 2 }
 
construct new(m) {
if (m <= BigInt.zero || !m.testBit(0)) {
Fiber.abort("Argument must be a positive, odd big integer.")
}
_m = m
_n = m.bitLength.toSmall
_rrm = (BigInt.one << (n*2)) % m
}
 
m { _m }
n { _n }
rrm { _rrm }
 
reduce(t) {
var a = t.copy()
for (i in 0..._n) {
if (a.testBit(0)) a = a + _m
a = a >> 1
}
if (a >= _m) a = a - _m
return a
}
}
 
var m = BigInt.new("750791094644726559640638407699")
var x1 = BigInt.new("540019781128412936473322405310")
var x2 = BigInt.new("515692107665463680305819378593")
 
var mont = Montgomery.new(m)
var t1 = x1 * mont.rrm
var t2 = x2 * mont.rrm
 
var r1 = mont.reduce(t1)
var r2 = mont.reduce(t2)
var r = BigInt.one << (mont.n)
 
System.print("b : %(Montgomery.base)")
System.print("n : %(mont.n)")
System.print("r : %(r)")
System.print("m : %(mont.m)")
System.print("t1: %(t1)")
System.print("t2: %(t2)")
System.print("r1: %(r1)")
System.print("r2: %(r2)")
System.print()
System.print("Original x1 : %(x1)")
System.print("Recovered from r1 : %(mont.reduce(r1))")
System.print("Original x2 : %(x2)")
System.print("Recovered from r2 : %(mont.reduce(r2))")
 
System.print("\nMontgomery computation of x1 ^ x2 mod m :")
var prod = mont.reduce(mont.rrm)
var base = mont.reduce(x1 * mont.rrm)
var exp = x2
while (exp.bitLength > 0) {
if (exp.testBit(0)) prod = mont.reduce(prod * base)
exp = exp >> 1
base = mont.reduce(base * base)
}
System.print(mont.reduce(prod))
System.print("\nLibrary-based computation of x1 ^ x2 mod m :")
System.print(x1.modPow(x2, m))</syntaxhighlight>
 
{{out}}
<pre>
b : 2
n : 100
r : 1267650600228229401496703205376
m : 750791094644726559640638407699
t1: 323165824550862327179367294465482435542970161392400401329100
t2: 308607334419945011411837686695175944083084270671482464168730
r1: 440160025148131680164261562101
r2: 435362628198191204145287283255
 
Original x1 : 540019781128412936473322405310
Recovered from r1 : 540019781128412936473322405310
Original x2 : 515692107665463680305819378593
Recovered from r2 : 515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m :
151232511393500655853002423778
 
Library-based computation of x1 ^ x2 mod m :
151232511393500655853002423778
</pre>
 
=={{header|zkl}}==
{{Trans|Go}}
Uses GMP (GNU Multi Precision library).
<syntaxhighlight lang="zkl">var [const] BN=Import("zklBigNum"); // libGMP
 
fcn montgomeryReduce(modulus,T){
_assert_(modulus.isOdd);
a:=BN(T); // we'll do in place math
do(modulus.len(2)){ // bits needed to hold modulus
if(a.isOdd) a.add(modulus);
a.div(2); // a>>=1
}
if(a>=modulus) a.sub(modulus);
a
}</syntaxhighlight>
<syntaxhighlight lang="zkl"> // magic numbers from the Go solution
//b:= 2;
//n:= 100;
//r:= BN("1267650600228229401496703205376");
m:= BN("750791094644726559640638407699");
 
t1:=BN("323165824550862327179367294465482435542970161392400401329100");
t2:=BN("308607334419945011411837686695175944083084270671482464168730");
 
r1:=BN("440160025148131680164261562101");
r2:=BN("435362628198191204145287283255");
 
x1:=BN("540019781128412936473322405310");
x2:=BN("515692107665463680305819378593");
 
// now demonstrate that it works:
println("Original x1: ", x1);
println("Recovererd from r1:",montgomeryReduce(m,r1));
println("Original x2: ", x2);
println("Recovererd from r2:", montgomeryReduce(m,r2));
 
// and demonstrate a use:
print("\nMontgomery computation of x1 ^ x2 mod m: ");
// this is the modular exponentiation algorithm, except we call
// montgomeryReduce instead of using a mod function.
prod:=montgomeryReduce(m,t1/x1); // 1
base:=montgomeryReduce(m,t1); // x1^1
exp :=BN(x2); // not reduced
while(exp){
if(exp.isOdd) prod=montgomeryReduce(m,prod.mul(base));
exp.div(2); // exp>>=1
base=montgomeryReduce(m,base.mul(base));
}
println(montgomeryReduce(m,prod));
println("Library-based computation of x1 ^ x2 mod m: ",x1.powm(x2,m));</syntaxhighlight>
{{out}}
<pre>
Original x1: 540019781128412936473322405310
Recovererd from r1:540019781128412936473322405310
Original x2: 515692107665463680305819378593
Recovererd from r2:515692107665463680305819378593
 
Montgomery computation of x1 ^ x2 mod m: 151232511393500655853002423778
Montgomery:
Library-based computation of x1 ^ x2 mod m: 599840174322403511400105423249151232511393500655853002423778
2.316ms
</pre>
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/210417...351381"