Find the last   40   decimal digits of   




a

b




{\displaystyle a^{b}}

,   where

   



a
=
2988348162058574136915891421498819466320163312926952423791023078876139


{\displaystyle a=2988348162058574136915891421498819466320163312926952423791023078876139}


   



b
=
2351399303373464486466122544523690094744975233415544072992656881240319


{\displaystyle b=2351399303373464486466122544523690094744975233415544072992656881240319}




A computer is too slow to find the entire value of   




a

b




{\displaystyle a^{b}}

.

Instead, the program must use a fast algorithm for modular exponentiation:   




a

b



mod


m


{\displaystyle a^{b}\mod m}

.

The algorithm must work for any integers   



a
,
b
,
m


{\displaystyle a,b,m}

,     where   



b
≥
0


{\displaystyle b\geq 0}

   and   



m
>
0


{\displaystyle m>0}

.




 
%%% For details of the algorithms used see
%%% https://en.wikipedia.org/wiki/Modular_exponentiation
 
-module modexp_rosetta.
-export [mod_exp/3,binary_exp/2,test/0].
 
mod_exp(Base,Exp,Mod) when
      is_integer(Base),    
      is_integer(Exp),
      is_integer(Mod),
      Base > 0,
      Exp > 0,
      Mod > 0 ->
    binary_exp_mod(Base,Exp,Mod).
 
binary_exp(Base,Exponent) ->
    binary_exp(Base,Exponent,1).
binary_exp(_,0,Result) ->
    Result;
binary_exp(Base,Exponent,Acc) ->
    binary_exp(Base*Base,Exponent bsr 1,Acc * exp_factor(Base,Exponent)).
 
 
binary_exp_mod(Base,Exponent,Mod) ->
    binary_exp_mod(Base rem Mod,Exponent,Mod,1).
binary_exp_mod(_,0,_,Result) ->
   Result;
binary_exp_mod(Base,Exponent,Mod,Acc) ->
    binary_exp_mod((Base*Base) rem Mod,
		   Exponent bsr 1,Mod,(Acc * exp_factor(Base,Exponent))rem Mod).
 
exp_factor(_,0) ->
    1;
exp_factor(Base,1) ->
    Base;
exp_factor(Base,Exponent) ->
    exp_factor(Base,Exponent band 1).
 
test() ->
    445 = mod_exp(4,13,497),
    %% Rosetta code example:
    mod_exp(2988348162058574136915891421498819466320163312926952423791023078876139,
	    2351399303373464486466122544523690094744975233415544072992656881240319,
	    binary_exp(10,40)).
 
 