Modular exponentiation
Find the last 40 decimal digits of , where
A computer is too slow to find the entire value of . Instead, the program must use a fast algorithm for modular exponentiation: .
The algorithm must work for any integers where and .
C
Given numbers are too big for even 64 bit integers, so might as well take the lazy route and use GMP: <lang c>#include <gmp.h>
int main() { mpz_t a, b, m, r;
mpz_init_set_str(a, "2988348162058574136915891421498819466320" "163312926952423791023078876139", 0); mpz_init_set_str(b, "2351399303373464486466122544523690094744" "975233415544072992656881240319", 0); mpz_init(m); mpz_ui_pow_ui(m, 10, 40);
- if 0 /* Don't actually need these two lines: hopefully GMP does it. */
mpz_mod(a, a, m); /* obvious */ mpz_mod(b, b, m); /* Fermat's little theorem */
- endif
mpz_init(r); mpz_powm(r, a, b, m);
gmp_printf("%Zd\n", r); /* 16808958343740453059 */
mpz_clear(a); mpz_clear(b); mpz_clear(m); mpz_clear(r); return 0; }</lang>
Go
<lang go>package main
import (
"fmt" "math/big"
)
func main() {
a, _ := new(big.Int).SetString( "2988348162058574136915891421498819466320163312926952423791023078876139", 10) b, _ := new(big.Int).SetString( "2351399303373464486466122544523690094744975233415544072992656881240319", 10) m, _ := new(big.Int).SetString("100000000000000000000", 10) r := new(big.Int).Exp(a, b, m) fmt.Println(r)
}</lang> Output:
16808958343740453059
J
Solution:<lang j> m&|@^</lang> Example:<lang j> a =: 2988348162058574136915891421498819466320163312926952423791023078876139x
b =: 2351399303373464486466122544523690094744975233415544072992656881240319x m =: 10^40x
a m&|@^ b
1527229998585248450016808958343740453059</lang> Discussion: The phrase m&|@^ is the natural expression of a^b mod m in J, and is recognized by the interpreter as an opportunity for optimization, by avoiding the exponentiation.
Java
<lang java>import java.math.BigInteger;
public class PowMod {
public static void main(String[] args){ BigInteger a = new BigInteger( "2988348162058574136915891421498819466320163312926952423791023078876139"); BigInteger b = new BigInteger( "2351399303373464486466122544523690094744975233415544072992656881240319"); BigInteger m = new BigInteger("10000000000000000000000000000000000000000"); System.out.println(a.modPow(b, m)); }
}</lang> Output:
1527229998585248450016808958343740453059
Perl
<lang perl>use bigint;
my $a = 2988348162058574136915891421498819466320163312926952423791023078876139; my $b = 2351399303373464486466122544523690094744975233415544072992656881240319; my $m = 10 ** 40; print $a->bmodpow($b, $m), "\n";</lang> Output:
1527229998585248450016808958343740453059
PHP
<lang php><?php $a = '2988348162058574136915891421498819466320163312926952423791023078876139'; $b = '2351399303373464486466122544523690094744975233415544072992656881240319'; $m = '1' . str_repeat('0', 40); echo bcpowmod($a, $b, $m), "\n"; ?></lang> Output:
1527229998585248450016808958343740453059
Python
<lang python>a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 print(pow(a, b, m))</lang> Output:
1527229998585248450016808958343740453059
Ruby
Ruby's core library has no modular exponentiation. OpenSSL, in Ruby's standard library, provides OpenSSL::BN#mod_exp. To reach this method, we call Integer#to_bn to convert a from Integer to OpenSSL::BN. The method implicitly converts b and m.
<lang ruby>require 'openssl'
a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 puts a.to_bn.mod_exp(b, m)</lang>
Or we can implement a custom method, Integer#rosetta_mod_exp, to calculate the same result. This method does exponentiation by successive squaring, but replaces each intermediate product with a congruent value.
<lang ruby>class Integer
def rosetta_mod_exp(exp, mod) exp < 0 and raise ArgumentError, "negative exponent" prod = 1 base = self % mod until exp.zero? exp[0].zero? or prod = (prod * base) % mod exp >>= 1 base = (base * base) % mod end prod end
end
a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 puts a.rosetta_mod_exp(b, m)</lang>
Tcl
While Tcl does have arbitrary-precision arithmetic (from 8.5 onwards), it doesn't expose a modular exponentiation function. Thus we implement one ourselves.
Recursive
<lang tcl>package require Tcl 8.5
- Algorithm from http://introcs.cs.princeton.edu/java/78crypto/ModExp.java.html
- but Tcl has arbitrary-width integers and an exponentiation operator, which
- helps simplify the code.
proc tcl::mathfunc::modexp {a b n} {
if {$b == 0} {return 1} set c [expr {modexp($a, $b / 2, $n)**2 % $n}] if {$b & 1} {
set c [expr {($c * $a) % $n}]
} return $c
}</lang> Demonstrating: <lang tcl>set a 2988348162058574136915891421498819466320163312926952423791023078876139 set b 2351399303373464486466122544523690094744975233415544072992656881240319 set n [expr {10**40}] puts [expr {modexp($a,$b,$n)}]</lang> Output:
1527229998585248450016808958343740453059
Iterative
<lang tcl>package require Tcl 8.5 proc modexp {a b n} {
set c 1 for {for {set t 1} {$t<$b} {set t [expr {$t<<1}]} {}} {$t} {set t [expr {$t>>1}]} {
if {$b & $t} { set c [expr {($c**2 * $a) % $n}] } else { set c [expr {$c**2 % $n}] }
} return $c
}</lang> Demonstrating: <lang tcl>set a 2988348162058574136915891421498819466320163312926952423791023078876139 set b 2351399303373464486466122544523690094744975233415544072992656881240319 set n [expr {10**40}] puts [modexp $a $b $n]</lang> Output:
1527229998585248450016808958343740453059
TXR
<lang txr>@(bind result @(exptmod 2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319 (expt 10 40)))</lang>
$ ./txr rosetta/modexp.txr result="1527229998585248450016808958343740453059"