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}

.




with javascript_semantics 
include mpfr.e
procedure mpz_mod_exp(mpz base, exponent, modulus, result)
    if mpz_cmp_si(exponent,1)=0 then
        mpz_set(result,base)
    else
        mpz _exp = mpz_init_set(exponent) -- (use a copy)
        bool odd = mpz_odd(_exp)
        if odd then
            mpz_sub_ui(_exp,_exp,1)
        end if
        mpz_fdiv_q_2exp(_exp,_exp,1)
        mpz_mod_exp(base,_exp,modulus,result)
        _exp = mpz_free(_exp)
        mpz_mul(result,result,result)
        if odd then
            mpz_mul(result,result,base)
        end if
    end if
    mpz_mod(result,result,modulus)
end procedure
 
mpz base     = mpz_init("2988348162058574136915891421498819466320163312926952423791023078876139"),
    exponent = mpz_init("2351399303373464486466122544523690094744975233415544072992656881240319"),
    modulus  = mpz_init("1"&repeat('0',40)),
    result   = mpz_init()
 
mpz_mod_exp(base,exponent,modulus,result)
?mpz_get_str(result)
 
-- check against the builtin:
mpz_powm(result,base,exponent,modulus)
?mpz_get_str(result)
