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}

.




/*REXX program  displays the  modular exponentiation  of:         a**b  mod  m          */
parse arg a b m                                       /*obtain optional args from the CL*/
if a=='' | a==","  then a= 2988348162058574136915891421498819466320163312926952423791023078876139
if b=='' | b==","  then b= 2351399303373464486466122544523690094744975233415544072992656881240319
if m ='' | m =","  then m= 40                         /*Not specified? Then use default.*/
say 'a=' a               "     ("length(a)  'digits)' /*display the value of A (&length)*/
say 'b=' b               "     ("length(b)  'digits)' /*   "     "    "    " B     "    */
       do j=1  for words(m);           y= word(m, j)  /*use one of the MM powers (list).*/
       say copies('═', linesize() - 1)                /*show a nice separator fence line*/
       say 'a**b (mod 10**'y")="  powerMod(a,b,10**y) /*display the answer ───► console.*/
       end   /*j*/
exit 0                                                /*stick a fork in it, we're done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
powerMod: procedure;  parse arg x,p,mm                /*fast modular exponentiation code*/
          parse value max(x*x, p, mm)'E0'  with "E" e /*obtain the biggest of the three.*/
          numeric digits max(40, e+e)                 /*ensure big enough to handle  A².*/
          $= 1                                        /*use this for the first value.   */
                  do  until p==0                      /*perform until P is equal to zero*/
                  if p // 2  then $= $ * x  //  mm    /*is P odd?  (is ÷ remainder ≡ 1?)*/
                  p= p % 2;       x= x * x  //  mm    /*halve  P;   calculate  x² mod n */
                  end   /*until*/;         return $   /* [↑]  keep mod'ing until P=zero.*/