=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;
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 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->redc($a * $b) if ($exp&1);      # a = a * b * R^-1 mod m = num * R mod m

        $b = $self->redc($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 (!$self->{modexpvars}) {
        $self->{r_mod_m}= $self->{r} % $self->{m};
        $self->{rr_mod_m}= ($self->{r} * $self->{r}) % $self->{m};
        $self->{modexpvars}= true;
    }

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

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

package main;
use bigint;
tst1();
tst2();
sub tst1 {
    printf("------------------\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("------------------\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);
}