Rosetta Code

Montgomery reduction: Difference between revisions

Page Discussion

Montgomery reduction (view source)

Revision as of 17:38, 19 June 2019

106 bytes removed ,  4 years ago
m
→‎{{header|Phix}}: replaced with gmp version
m (→‎{{header|Phix}}: replaced with gmp version)
Line 779:
=={{header|Phix}}==
{{trans|D}}
{{libheader|mpfr}}
<lang Phix>include builtins\bigintmpfr.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, bigintmpz a)
integer n = mont[BITLEN]
bigintmpz m = mont[MODULUS],
r = mpz_init_set(a)
for i=1 to n do
if bi_modmpz_odd(a,2r)=BI_ONE then -- odd
a = bi_addmpz_add(ar,r,m)
end if
{a,?} = bi_div3mpz_fdiv_q_ui(ar,r,2)
end for
if bi_comparempz_cmp(ar,m)>=0 then a = bi_submpz_sub(ar,r,m) end if
return ar
end function
 
function Montgomery(bigintmpz m)
if bi_signmpz_sign(m)=-1 then crash("must be positive") end if
if bi_modnot mpz_odd(m,2)!=BI_ONE then crash("must be odd") end if
integer n = bitLengthmpz_sizeinbase(m,2)
bigintmpz rrm = bi_modmpz_init(bi_shl(1,n*2),m)
mpz_powm_ui(rrm,rrm,n*2,m)
return {2, -- BASE
n, -- BITLEN
Line 818 ⟶ 811:
end function
bigintmpz m = bi_newmpz_init("750791094644726559640638407699"),
x1 = bi_newmpz_init("540019781128412936473322405310"),
x2 = bi_newmpz_init("515692107665463680305819378593");,
t1 = mpz_init(),
t2 = mpz_init()
sequence mont = Montgomery(m)
bigint t1 = bi_mulmpz_mul(t1,x1,mont[RRM]),
t2 = bi_mulmpz_mul(t2,x2,mont[RRM]),
mpz r1 = reduce(mont,t1),
r2 = reduce(mont,t2),
r = bi_shlmpz_init(1,mont[BITLEN])
mpz_ui_pow_ui(r,2,mont[BITLEN])
printf(1,"b : %d\n", {mont[BASE]})
printf(1,"n : %d\n", {mont[BITLEN]})
printf(1,"r : %s\n", {bi_sprintmpz_get_str(r)})
printf(1,"m : %s\n", {bi_sprintmpz_get_str(mont[MODULUS])})
printf(1,"t1: %s\n", {bi_sprintmpz_get_str(t1)})
printf(1,"t2: %s\n", {bi_sprintmpz_get_str(t2)})
printf(1,"r1: %s\n", {bi_sprintmpz_get_str(r1)})
printf(1,"r2: %s\n", {bi_sprintmpz_get_str(r2)})
printf(1,"\n")
printf(1,"Original x1 : %s\n", {bi_sprintmpz_get_str(x1)})
printf(1,"Recovered from r1 : %s\n", {bi_sprintmpz_get_str(reduce(mont,r1))})
printf(1,"Original x2 : %s\n", {bi_sprintmpz_get_str(x2)})
printf(1,"Recovered from r2 : %s\n", {bi_sprintmpz_get_str(reduce(mont,r2))})
printf(1,"\nMontgomery computation of x1 ^ x2 mod m :")
bigintmpz prod = reduce(mont,mont[RRM]),
base = reducempz_mul(montr,bi_mul(x1,mont[RRM])),
mpz base = reduce(mont,r),
expn = mpz_init_set(x2;)
while bitLength(expn)>0 do
 
if bi_mod(expn,2)=BI_ONE then -- odd▼
▲ while bi_comparempz_cmp_si(vexpn,0)!=0 do
prod = reduce(mont,bi_mul(prod,base))▼
▲ if bi_modmpz_odd(expn,2)=BI_ONE then -- odd
mpz_mul(prod,prod,base)
▲ {v,?}prod = bi_div3reduce(vmont, 2prod)
end if
{expn,?} = bi_div3mpz_fdiv_q_ui(expn,expn,2)
mpz_mul(base = reduce(mont,bi_mul(base,base))
▲ prodbase = reduce(mont,bi_mul(prod,base))
end while
printf(1,"%s\n",{bi_sprintmpz_get_str(reduce(mont,prod))})
printf(1,"\n alternate computation of x1 ^ x2 mod m :");
mpz_powm(r,x1,x2,m)
printf(1,"%s\n",{bi_sprintmpz_get_str(bi_mod_exp(x1, x2, m)r)})</lang>
{{out}}
<pre>
Petelomax
7,795

edits

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