use strict;
use warnings;
=pod

 c = a^b (mod n)

 N > 0
 R,T  such that:
    * R>N 
    * gcd(N,R)=1
    * 0 <= T < N*R
    * usually: R is some power of the base system word

    * T*R^-1(mod N)  is the montgomery reduction of T

    * 0 < R^-1 (mod N) < N
    * 0 < N' < R
    * R*R^-1 - N*N' = 1   : R^-1 is the modular inverse of R (mod N)

example: R=2^32, N=1234567=0x12D687
 ->  R'=0xFBE4C, N'=0xD5F1F0C9

 -> montgomery reduction of T mod N (wrt R) is:

 REDC(T)   // = T*R' mod m
 {
    m = (T mod R)*N' mod R   such that  0 <= m < R
    t = (T + m*N)/R
    if (t>=N) t = t-N

    return t;
 }
example: REDC(123) = (123*0xFBE4C)%1234567 = 0xEF4BA
 m = (123 % 2^32) * 0xD5F1F0C9 = 0x66CB3EB093
    m[0] = 0xCB3EB093
 t*R = (T+m[0]*N) = 123+0xCB3EB093*1234567 = 0xEF4BA00000000
 -> t = 0xEF4BA


 because m*N = T*N'*N = T*(R*R^-1  -1 ) = -T ( mod R )

 -> t = (T+m*N)/R = (T+k*R-T)/R = k  is an integer

 and t*R = T ( mod N )   ->  t = T*R^-1 (mod N)
 0 <= T+m*N < R*N+R*N  ( because T<N*R and m < R )
    -> 0 <= t < 2*N

 REDC(x*y*R)  = x*y (mod N)



--------------------------------

n' = -n^-1 (mod R)
gcd(n,R)=1

->  ( x+t*n ) / R = x*R^-1 (mod n)
 and t = x*n' (mod R)



x=sum(x_i*b^i)  and y=sum(y_i*b^i) 

calc x*y*R^-1 (mod n)

A=0
for i=0 .. n-1
    u =  ( a_0+x_i*y_0 ) * n'  ( mod b )
    A += (A + x_i*y + u * n ) / b
    A -= n  if (A > n)

=cut

package Montgomery;
use bigint lib=>'GMP';
use strict;
sub new {
    my ($class, $m, $b, $n)= @_;
    printf("r=%s^%d\nm=%s\n", $b->as_hex, $n, $m->as_hex);
    my $r= $b**$n;
    my $minv= (-$m)->bmodinv($r);
    printf("gcd(m,r)=%s\n", Math::BigInt::bgcd(-$m,$r)->as_hex);
    return bless {
        m=>$m,
        r=>$r,
        b=>$b,
        m_inv=>$minv,
    }, $class;
}

# REDC(x*R)  = x (mod N)
# REDC(REDC(x*y)) = REDC(x)*REDC(y)  (mod N)
# REDC(x*y) = R * REDC(x)*REDC(y)  (mod N)

# algorithm hac-14.32
# m={n digits base b}
# gcd(m,b)=1
# R=b^n
# m'=m^-1(mod b)
# T ={2n digits base b}  < m*R
sub redc {
    my ($self, $x)=@_;

    my $m = ( $x * $self->{m_inv} )->bmod( $self->{r} );
    my $t = ($x + $m * $self->{m}) / $self->{r};
    if ($t>=$self->{m}) { $t-=$self->{m}; }

    return $t;
}
sub mult {
    my ($self, $a, $b)= @_;
    return $self->redc($a * $b);
}
sub modpow {
    my ($self, $num, $exp)= @_;

    my $a= $self->{r};                              # a = R
    my $b= ( $num * $self->{r} ) % $self->{m};      # b = num * R  mod m
printf("modpow: num=%s  e=%s  -> A=%s, num~=%s\n", $num->as_hex, $exp->as_hex, $a->as_hex, $b->as_hex);
    while ($exp)
    {
        $a = $self->mult($a, $b) if ($exp&1);      # a = a * b * R^-1 mod m = num * R mod m

        $b = $self->mult($b, $b);                  # b = b * b * R^-1 mod m = num * num * R mod m
        $exp>>=1;
    }

    return $self->redc($a);                         # return a * R^-1  = num
}
sub modpow2 {
    my ($self, $x, $e)= @_;
    if (!defined $self->{r_mod_m} || !defined $self->{rr_mod_m}) {
        $self->{r_mod_m}= $self->{r} % $self->{m};
        $self->{rr_mod_m}= ($self->{r} * $self->{r}) % $self->{m};
    }

    my $xt= $self->mult($x, $self->{rr_mod_m});
    my $a= $self->{r_mod_m}; 
printf("modpow: num=%s  e=%s  -> A=%s, x~=%s\n", $x->as_hex, $e->as_hex, $a->as_hex, $xt->as_hex);
    while ($e)
    {
        $a = $self->mult($a, $xt) if ($e&1);
        $xt = $self->mult($xt, $xt);
        $e>>=1;
    }

    return $self->redc($a);
}
sub modpow3 {
    my ($self, $x, $e)= @_;
    if (!defined $self->{r_mod_m} || !defined $self->{rr_mod_m}) {
        $self->{r_mod_m}= $self->{r} % $self->{m};
        $self->{rr_mod_m}= ($self->{r} * $self->{r}) % $self->{m};
    }

    my $xt= $self->mult($x, $self->{rr_mod_m});
    my $a= $self->{r_mod_m}; 
printf("modpow: num=%s  e=%s  -> A=%s, x~=%s\n", $x->as_hex, $e->as_hex, $a->as_hex, $xt->as_hex);
    my $bit= 1;
    while ($bit<$e) {
        $bit <<= 1;
    }
    
    while ($bit>=1)
    {
        $a = $self->mult($a, $a);
        $a = $self->mult($a, $xt) if ($e & $bit);
        $bit >>= 1;
    }

    return $self->redc($a);
}

package main;
use bigint lib=>'GMP';
use strict;
tst1();
tst2();
tst3();
tst4();
sub tst1 {
    printf("------------------ tst1\n");
    my $m= Montgomery->new(1234567, 2**32, 1);
    printf("r= %s\n", $m->{r}->as_hex);
    printf("m_inv= %s\n", $m->{m_inv}->as_hex);
    printf("r(x*R) mod m= %s   x mod m= %s\n", 
        ($m->redc(123456789*2**32) %1234567)->as_hex, 
        (123456789 % 1234567)->as_hex);

    printf("123= %s  456= %s  123*456= %s\n",
        $m->redc(123)->as_hex, $m->redc(456)->as_hex, $m->redc(123*456)->as_hex);

    printf("r(123)*r(456)*R= %s    r(123*456) (modN)= %s\n", 
        (($m->redc(123)*$m->redc(456)*2**32)%1234567)->as_hex, 
        ($m->redc(123*456)%1234567)->as_hex);

    printf("montg= %s  %s normal= %s\n", 
        $m->modpow(123, 31)->as_hex, $m->modpow2(123, 31)->as_hex, 123->bmodpow(31, 1234567)->as_hex);
}
sub tst2 {
    printf("------------------ tst2\n");

    my $m= Montgomery->new(0xF765A3A0C9C291D81A56FE73794A746B8DA23DBE155D0D495B49D581B5C6545F449A10FDF1C26A92FBD1F43A0687044927A6A21B69A73999E6083D03ACDAFFA6409F1BC71D810628F6E18F76231ED6E22D54ED2502E66F8A33D0D5F07B3EB605F7418110E2EF9A5EE77B070F4EADFCF3D70C53E870F29C9D4F229F2CB6C25383, 2**32, 32);
    printf("r=     %s\n", $m->{r}->as_hex);
    printf("m_inv= %s\n", $m->{m_inv}->as_hex);
    printf("montg= %s\n", 
        $m->modpow(0x0857bde36f680cda9535784bc4aec5f4344131071b419f732ac9c74d0e61db49dd958c7344236e0279df009c6e66aec6ba574c2820d4aeb0c4d814c8e184c6ea7e6d8aa3e15d1c251c78c5364ea2b3edb3c19e90739afa765506242e78fcdc71a87efdfe2df6ce6039fc62cb3b360cb77cd5574292282df352886cbc3fcfbff2, 0x10001)->as_hex);
    printf("mont2= %s\n", 
        $m->modpow2(0x0857bde36f680cda9535784bc4aec5f4344131071b419f732ac9c74d0e61db49dd958c7344236e0279df009c6e66aec6ba574c2820d4aeb0c4d814c8e184c6ea7e6d8aa3e15d1c251c78c5364ea2b3edb3c19e90739afa765506242e78fcdc71a87efdfe2df6ce6039fc62cb3b360cb77cd5574292282df352886cbc3fcfbff2, 0x10001)->as_hex);
}

sub tst3 {
    printf("------------------ tst3\n");
    my $m= Montgomery->new(0xc386528327edd25ad50f3458b272d0ee6375276d56b029ca8b56baa1d8237929, 2**32, 8);
    printf("montg= %s\n", $m->modpow(0x8d0071d6b068fd5da6e5137b26cfcbea19abda46ce6f258dd539bf0c9102e5c6, 0x10001)->as_hex);
    printf("mont2= %s\n", $m->modpow2(0x8d0071d6b068fd5da6e5137b26cfcbea19abda46ce6f258dd539bf0c9102e5c6, 0x10001)->as_hex);
    printf("mont3= %s\n", $m->modpow3(0x8d0071d6b068fd5da6e5137b26cfcbea19abda46ce6f258dd539bf0c9102e5c6, 0x10001)->as_hex);
    printf("rr=    %s\n", $m->{rr_mod_m}->as_hex);
    printf("m_inv= %s\n", $m->{m_inv}->as_hex);
}
sub tst4 {
    printf("------------------ tst4\n");
    my $m= Montgomery->new(0xFFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F24117C4B1FE649286651ECE45B3DC2007CB8A163BF0598DA48361C55D39A69163FA8FD24CF5F83655D23DCA3AD961C62F356208552BB9ED529077096966D670C354E4ABC9804F1746C08CA18217C32905E462E36CE3BE39E772C180E86039B2783A2EC07A28FB5C55DF06F4C52C9DE2BCBF6955817183995497CEA956AE515D2261898FA051015728E5A8AAAC42DAD33170D04507A33A85521ABDF1CBA64ECFB850458DBEF0A8AEA71575D060C7DB3970F85A6E1E4C7ABF5AE8CDB0933D71E8C94E04A25619DCEE3D2261AD2EE6BF12FFA06D98A0864D87602733EC86A64521F2B18177B200CBBE117577A615D6C770988C0BAD946E208E24FA074E5AB3143DB5BFCE0FD108E4B82D120A92108011A723C12A787E6D788719A10BDBA5B2699C327186AF4E23C1A946834B6150BDA2583E9CA2AD44CE8DBBBC2DB04DE8EF92E8EFC141FBECAA6287C59474E6BC05D99B2964FA090C3A2233BA186515BE7ED1F612970CEE2D7AFB81BDD762170481CD0069127D5B05AA993B4EA988D8FDDC186FFB7DC90A6C08F4DF435C934063199FFFFFFFFFFFFFFFF, 2**32, 128);
    printf("montg= %s\n", $m->modpow(0xf2b27ddf8fbb4d435a07b659406374ad99636ed953187f1a9fa1462bc9df8d4634fbaa1a4a5133e1991bdcfb150bb3da0389806f7fc98768446f5e426b9711cac8426bd79fb80bb626fdb73198e57984804730579f76e4fec7e1561141bae3dc73fb7d76b6b885ee79c69ac5716947eabf82d4b7c5d08d2ae91aaeb43f3405d619826e8ec3cab341088c88bdecd9fa8d6ae240d0c5953d2b6800554fad4798087a78b878fe91dbd65fb19497deb4a89a9360083f0571fd68ae5d882c7517ac43a5cecccdb00c144914edf49a4382766a8af87c11a2714fde33bb1bb410723673eed6e0827a84494561f0768ef15bc2fe73ef3582928fc09ea6bd7b8fa53fac77470ab11135f9214d6dc2b6ca8eee9d68569e4e898eaad62b5fe952dd6a25731d546ab84e718a8bd3b81b172e8f0984232b6725ce45032513641ecc752d744aa6747a673abd9dd3b462d8551c0a3299943a082000431a1bca79cf6ef22644a9a21a4a7a39135775c99338ca458fff4e9cbf45aba630396c6af6aa0b006c5cbef8d6ba3929e29c2b1200ba5e3082f5d5939f7fd1d44b6a8d43f3c721eb418b47226e79c4f68b30741a91dd625a83587228346ae7453fd7be74d7fbb0358beec6a91b1cde28b62ed33dad9f1e1d2aca384ab2a0363121b93f0b605ac6ed07e24ea879e7adec07341b54edd7729adbb5dcee8688fa77a253813af297063181741f2a, 0x3f819317dc1b4562be66998481db09a1c20e0a55f56ff819db7d565b08799cd4)->as_hex);
    printf("mont2= %s\n", $m->modpow2(0xf2b27ddf8fbb4d435a07b659406374ad99636ed953187f1a9fa1462bc9df8d4634fbaa1a4a5133e1991bdcfb150bb3da0389806f7fc98768446f5e426b9711cac8426bd79fb80bb626fdb73198e57984804730579f76e4fec7e1561141bae3dc73fb7d76b6b885ee79c69ac5716947eabf82d4b7c5d08d2ae91aaeb43f3405d619826e8ec3cab341088c88bdecd9fa8d6ae240d0c5953d2b6800554fad4798087a78b878fe91dbd65fb19497deb4a89a9360083f0571fd68ae5d882c7517ac43a5cecccdb00c144914edf49a4382766a8af87c11a2714fde33bb1bb410723673eed6e0827a84494561f0768ef15bc2fe73ef3582928fc09ea6bd7b8fa53fac77470ab11135f9214d6dc2b6ca8eee9d68569e4e898eaad62b5fe952dd6a25731d546ab84e718a8bd3b81b172e8f0984232b6725ce45032513641ecc752d744aa6747a673abd9dd3b462d8551c0a3299943a082000431a1bca79cf6ef22644a9a21a4a7a39135775c99338ca458fff4e9cbf45aba630396c6af6aa0b006c5cbef8d6ba3929e29c2b1200ba5e3082f5d5939f7fd1d44b6a8d43f3c721eb418b47226e79c4f68b30741a91dd625a83587228346ae7453fd7be74d7fbb0358beec6a91b1cde28b62ed33dad9f1e1d2aca384ab2a0363121b93f0b605ac6ed07e24ea879e7adec07341b54edd7729adbb5dcee8688fa77a253813af297063181741f2a, 0x3f819317dc1b4562be66998481db09a1c20e0a55f56ff819db7d565b08799cd4)->as_hex);
    printf("rr=    %s\n", $m->{rr_mod_m}->as_hex);
    printf("m_inv= %s\n", $m->{m_inv}->as_hex);
}