Montgomery reduction: Difference between revisions
→{{header|Phix}}
(→{{header|Perl}}: Formatting, add libheader) |
|||
Line 775:
Montgomery computation x1**x2 mod m: 151232511393500655853002423778
Built-in op computation x1**x2 mod m: 151232511393500655853002423778
</pre>
=={{header|Phix}}==
{{trans|D}}
<lang Phix>include builtins\bigint.e
function bitLength(bigint v)
if bi_sign(v)=-1 then v = bi_neg(v) end if
integer result = 0
while bi_compare(v,0)!=0 do
{v,?} = bi_div3(v, 2)
result += 1
end while
return result
end function
enum BASE, BITLEN, MODULUS, RRM
function reduce(sequence mont, bigint a)
integer n = mont[BITLEN]
bigint m = mont[MODULUS]
for i=1 to n do
if bi_mod(a,2)=BI_ONE then -- odd
a = bi_add(a,m)
end if
{a,?} = bi_div3(a,2)
end for
if bi_compare(a,m)>=0 then a = bi_sub(a,m) end if
return a
end function
function Montgomery(bigint m)
if bi_sign(m)=-1 then crash("must be positive") end if
if bi_mod(m,2)!=BI_ONE then crash("must be odd") end if
integer n = bitLength(m)
bigint rrm = bi_mod(bi_shl(1,n*2),m)
return {2, -- BASE
n, -- BITLEN
m, -- MODULUS
rrm -- 1<<(n*2) % m
}
end function
bigint m = bi_new("750791094644726559640638407699"),
x1 = bi_new("540019781128412936473322405310"),
x2 = bi_new("515692107665463680305819378593");
sequence mont = Montgomery(m)
bigint t1 = bi_mul(x1,mont[RRM]),
t2 = bi_mul(x2,mont[RRM]),
r1 = reduce(mont,t1),
r2 = reduce(mont,t2),
r = bi_shl(1,mont[BITLEN])
printf(1,"b : %d\n", {mont[BASE]})
printf(1,"n : %d\n", {mont[BITLEN]})
printf(1,"r : %s\n", {bi_sprint(r)})
printf(1,"m : %s\n", {bi_sprint(mont[MODULUS])})
printf(1,"t1: %s\n", {bi_sprint(t1)})
printf(1,"t2: %s\n", {bi_sprint(t2)})
printf(1,"r1: %s\n", {bi_sprint(r1)})
printf(1,"r2: %s\n", {bi_sprint(r2)})
printf(1,"\n")
printf(1,"Original x1 : %s\n", {bi_sprint(x1)})
printf(1,"Recovered from r1 : %s\n", {bi_sprint(reduce(mont,r1))})
printf(1,"Original x2 : %s\n", {bi_sprint(x2)})
printf(1,"Recovered from r2 : %s\n", {bi_sprint(reduce(mont,r2))})
printf(1,"\nMontgomery computation of x1 ^ x2 mod m :")
bigint prod = reduce(mont,mont[RRM]),
base = reduce(mont,bi_mul(x1,mont[RRM])),
expn = x2;
while bitLength(expn)>0 do
if bi_mod(expn,2)=BI_ONE then -- odd
prod = reduce(mont,bi_mul(prod,base))
end if
{expn,?} = bi_div3(expn,2)
base = reduce(mont,bi_mul(base,base))
end while
printf(1,"%s\n",{bi_sprint(reduce(mont,prod))})
printf(1,"\n alternate computation of x1 ^ x2 mod m :");
printf(1,"%s\n",{bi_sprint(bi_mod_exp(x1, x2, m))})</lang>
{{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>
| |||