# Modular exponentiation

Modular exponentiation
You are encouraged to solve this task according to the task description, using any language you may know.

Find the last 40 decimal digits of ${\displaystyle a^{b}}$, where

• ${\displaystyle a=2988348162058574136915891421498819466320163312926952423791023078876139}$
• ${\displaystyle b=2351399303373464486466122544523690094744975233415544072992656881240319}$

A computer is too slow to find the entire value of ${\displaystyle a^{b}}$. Instead, the program must use a fast algorithm for modular exponentiation: ${\displaystyle a^{b}\mod m}$.

The algorithm must work for any integers ${\displaystyle a,b,m}$ where ${\displaystyle b\geq 0}$ and ${\displaystyle m>0}$.

Using the big integer implementation from a cryptographic library [1].

procedure Mod_Exp is

```  A: String :=
"2988348162058574136915891421498819466320163312926952423791023078876139";
B: String :=
"2351399303373464486466122544523690094744975233415544072992656881240319";
```
```  D: constant Positive := Positive'Max(Positive'Max(A'Length, B'Length), 40);
-- the number of decimals to store A, B, and result
Bits: constant Positive := (34*D)/10;
-- (slightly more than) the number of bits to store A, B, and result
package LN is new Crypto.Types.Big_Numbers (Bits + (32 - Bits mod 32));
-- the actual number of bits has to be a multiple of 32
use type LN.Big_Unsigned;
```
```  function "+"(S: String) return LN.Big_Unsigned
renames LN.Utils.To_Big_Unsigned;
```
```  M: LN.Big_Unsigned := (+"10") ** (+"40");
```

begin

```  Ada.Text_IO.Put("A**B (mod 10**40) = ");
```

end Mod_Exp;</lang>

Output:
`A**B (mod 10**40) = 1527229998585248450016808958343740453059`

## Bracmat

Translation of: Icon_and_Unicon

<lang bracmat> ( ( mod-power

```   =   base exponent modulus result
.   !arg:(?base,?exponent,?modulus)
& !exponent:~<0
& 1:?result
&   whl
' ( !exponent:>0
&     ( (   mod\$(!exponent.2):1
& mod\$(!result*!base.!modulus):?result
& -1
| 0
)
+ !exponent
)
* 1/2
: ?exponent
& mod\$(!base^2.!modulus):?base
)
& !result
)
& ( a
= 2988348162058574136915891421498819466320163312926952423791023078876139
)
& ( b
= 2351399303373464486466122544523690094744975233415544072992656881240319
)
& out\$("last 40 digits = " mod-power\$(!a,!b,10^40))
)</lang>
```

Output:

`last 40 digits =  1527229998585248450016808958343740453059`

## BBC BASIC

Uses the Huge Integer Math & Encryption library. <lang bbcbasic> INSTALL @lib\$+"HIMELIB"

```     PROC_himeinit("")

PROC_hiputdec(1, "2988348162058574136915891421498819466320163312926952423791023078876139")
PROC_hiputdec(2, "2351399303373464486466122544523690094744975233415544072992656881240319")
PROC_hiputdec(3, "10000000000000000000000000000000000000000")
h1% = 1 : h2% = 2 : h3% = 3 : h4% = 4
SYS `hi_PowMod`, ^h1%, ^h2%, ^h3%, ^h4%
PRINT FN_higetdec(4)</lang>
```

Output:

```1527229998585248450016808958343740453059
```

## C

Given numbers are too big for even 64 bit integers, so might as well take the lazy route and use GMP:

Library: 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);

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>

## Common Lisp

<lang lisp>(defun rosetta-mod-expt (base power divisor)

``` "Return BASE raised to the POWER, modulo DIVISOR.
This function is faster than (MOD (EXPT BASE POWER) DIVISOR), but
only works when POWER is a non-negative integer."
(setq base (mod base divisor))
;; Multiply product with base until power is zero.
(do ((product 1))
((zerop power) product)
;; Square base, and divide power by 2, until power becomes odd.
(do () ((oddp power))
(setq base (mod (* base base) divisor)
```

power (ash power -1)))

```   (setq product (mod (* product base) divisor)
```

power (1- power))))

(let ((a 2988348162058574136915891421498819466320163312926952423791023078876139)

```     (b 2351399303373464486466122544523690094744975233415544072992656881240319))
(format t "~A~%" (rosetta-mod-expt a b (expt 10 40))))</lang>
```
Works with: CLISP

<lang lisp>;; CLISP provides EXT:MOD-EXPT (let ((a 2988348162058574136915891421498819466320163312926952423791023078876139)

```     (b 2351399303373464486466122544523690094744975233415544072992656881240319))
(format t "~A~%" (mod-expt a b (expt 10 40))))</lang>
```

### Implementation with LOOP

<lang lisp>(defun mod-expt (a n m)

```  (loop with c = 1 while (plusp n) do
(if (oddp n) (setf c (mod (* a c) m)))
(setf n (ash n -1))
(setf a (mod (* a a) m))
finally (return c)))</lang>
```

## D

Translation of: Icon_and_Unicon

<lang d>import std.stdio, std.bigint;

BigInt powMod(BigInt base, BigInt exponent, BigInt modulus) in {

```  assert(exponent >= 0);
```

} body {

```   BigInt result = 1;
while (exponent > 0) {
if (exponent % 2 == 1)
result = (result * base) % modulus;
exponent /= 2;
base = base ^^ 2 % modulus;
}
return result;
```

}

void main() {

```   powMod(BigInt("29883481620585741369158914214988194" ~
"66320163312926952423791023078876139"),
BigInt("235139930337346448646612254452369009" ~
"4744975233415544072992656881240319"),
BigInt(10) ^^ 40).writeln();
```

}</lang>

Output:
`1527229998585248450016808958343740453059`

## Dc

<lang Dc>2988348162058574136915891421498819466320163312926952423791023078876139 2351399303373464486466122544523690094744975233415544072992656881240319 10 40^|p</lang>

## Emacs Lisp

Library: Calc

<lang lisp>(let ((a "2988348162058574136915891421498819466320163312926952423791023078876139")

```     (b "2351399303373464486466122544523690094744975233415544072992656881240319"))
;; "\$ ^ \$\$ mod (10 ^ 40)" performs modular exponentiation.
;; "unpack(-5, x)_1" unpacks the integer from the modulo form.
(message "%s" (calc-eval "unpack(-5, \$ ^ \$\$ mod (10 ^ 40))_1" nil a b)))</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 := big.NewInt(10)
r := big.NewInt(40)
m.Exp(m, r, nil)
```
```   r.Exp(a, b, m)
fmt.Println(r)
```

}</lang> Output:

```1527229998585248450016808958343740453059
```

Kind of a hack. We partially implement a "modular arithmetic" instance of `Num`, so that we can take advantage of the efficient built-in exponentiation-by-squaring operation without implementing it ourselves. Since there are no "local" instances, we must keep the modulo base around with us in the type, which makes the code inelegant. <lang haskell>-- Private type. Do not use outside of the modPow function newtype ModN = ModN (Integer, Integer) deriving (Eq, Show) instance Num ModN where

``` -- actually only multiplication needs to be implemented
-- but we do some of the other ones too for good measure
ModN (x, m) + ModN (y, m') | m == m' = ModN ((x + y) `mod` m, m)
| otherwise = undefined
ModN (x, m) * ModN (y, m') | m == m' = ModN ((x * y) `mod` m, m)
| otherwise = undefined
negate (ModN (x, m)) = ModN ((- x) `mod` m, m)
abs _ = undefined
signum _ = undefined
fromInteger _ = undefined
```

modPow :: Integer -> Integer -> Integer -> Integer modPow _ 0 m = 1 `mod` m modPow a b m = c

``` where a' = ModN (a, m)
ModN (c, _) = a' ^ b
```

main :: IO () main = print \$ modPow a b m

``` where a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ^ 40</lang>
```

Output:

`1527229998585248450016808958343740453059`

## Icon and Unicon

This uses the exponentiation procedure from RSA Code an example of the right to left binary method. <lang Icon>procedure main()

```   a := 2988348162058574136915891421498819466320163312926952423791023078876139
b := 2351399303373464486466122544523690094744975233415544072992656881240319
write("last 40 digits = ",mod_power(a,b,(10^40))
```

end

procedure mod_power(base, exponent, modulus) # fast modular exponentation

```  if exponent < 0 then runerr(205,m)          # added for this task
result := 1
while exponent > 0 do {
if exponent % 2 = 1 then
result := (result * base) % modulus
exponent /:= 2
base := base ^ 2 % modulus
}
return result
```

end</lang>

Output:
`last 40 digits = 1527229998585248450016808958343740453059`

## 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

`java.math.BigInteger.modPow` solves this task. Inside OpenJDK, BigInteger.java implements `BigInteger.modPow` with a fast algorithm from Colin Plumb's bnlib. This "window algorithm" caches odd powers of the base, to decrease the number of squares and multiplications. It also exploits both the Chinese remainder theorem and the Montgomery reduction.

<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`

## Maple

<lang Maple>a := 2988348162058574136915891421498819466320163312926952423791023078876139: b := 2351399303373464486466122544523690094744975233415544072992656881240319: a &^ b mod 10^40;</lang> Output:

`1527229998585248450016808958343740453059`

## Mathematica

<lang Mathematica>a = 2988348162058574136915891421498819466320163312926952423791023078876139; b = 2351399303373464486466122544523690094744975233415544072992656881240319; m = 10^40; PowerMod[a, b, m] -> 1527229998585248450016808958343740453059</lang>

## Maxima

<lang maxima>a: 2988348162058574136915891421498819466320163312926952423791023078876139\$ b: 2351399303373464486466122544523690094744975233415544072992656881240319\$ power_mod(a, b, 10^40); /* 1527229998585248450016808958343740453059 */</lang>

## PARI/GP

<lang parigp>a=2988348162058574136915891421498819466320163312926952423791023078876139; b=2351399303373464486466122544523690094744975233415544072992656881240319; lift(Mod(a,10^40)^b)</lang>

## Pascal

Works with: Free_Pascal
Library: GMP

A port of the C example using gmp. <lang pascal>Program ModularExponentiation(output);

uses

``` gmp;

```

var

``` a, b, m, r: mpz_t;
fmt: pchar;
```

begin

``` mpz_init_set_str(a, '2988348162058574136915891421498819466320163312926952423791023078876139', 10);
mpz_init_set_str(b, '2351399303373464486466122544523690094744975233415544072992656881240319', 10);
mpz_init(m);
mpz_ui_pow_ui(m, 10, 40);
```
``` mpz_init(r);
mpz_powm(r, a, b, m);
```
``` fmt := '%Zd' + chr(13) + chr(10);
mp_printf(fmt, @r); (* ...16808958343740453059 *)

mpz_clear(a);
mpz_clear(b);
mpz_clear(m);
mpz_clear(r);
```

end.</lang> Output:

```% ./ModularExponentiation
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`

## Perl 6

This is specced as a built-in, but here's an explicit version: <lang perl6>sub expmod(Int \$a is copy, Int \$b is copy, \$n) {

```   my \$c = 1;
repeat while \$b div= 2 {
(\$c *= \$a) %= \$n if \$b % 2;
(\$a *= \$a) %= \$n;
}
\$c;
```

}

say expmod

```   2988348162058574136915891421498819466320163312926952423791023078876139,
2351399303373464486466122544523690094744975233415544072992656881240319,
10**40;</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`

## PicoLisp

The following function is taken from "lib/rsa.l": <lang PicoLisp>(de **Mod (X Y N)

```  (let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )</lang>
```

Test: <lang PicoLisp>: (**Mod

```  2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
10000000000000000000000000000000000000000 )
```

-> 1527229998585248450016808958343740453059</lang>

## Python

<lang python>a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 print(pow(a, b, m))</lang> Output:

`1527229998585248450016808958343740453059`

## REXX

This REXX program attempts to handle any a,b, or m, but there are limits for any computer language.
For REXX, it's around eight million digits, unless ${\displaystyle a^{2}}$ or ${\displaystyle 10^{m}}$ exceeds that. <lang rexx>/*REXX program to show modular exponentation: a**b mod M */ parse arg a b mm /*get the arguments (maybe).*/ if a== | a==',' then a=,

```   2988348162058574136915891421498819466320163312926952423791023078876139
```

if b== | b==',' then b=,

```   2351399303373464486466122544523690094744975233415544072992656881240319
```

if mm== then mm=40 say 'a=' a; say ' ('length(a) "digits)" say 'b=' b; say ' ('length(b) "digits)"

```     do j=1 for words(mm);   m=word(mm,j);   say copies('─',linesize()-1)
say 'a**b (mod 10**'m")=" powerModulated(a,b,10**m)
end   /*j*/
```

exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────────POWERMODULATED subroutine───────*/ powerModulated: procedure; parse arg x,p,n /*fast modular exponentation*/ if p==0 then return 1 /*special case. */ if p==1 then return x /*special case. */ if p<0 then do; say '***error!*** power is negative:' p; exit 13; end parse value max(x**2,p,n)'E0' with "E" e /*pick biggest of the three.*/ numeric digits max(20,e*2) /*big enough to handle A² */ _=1

```         do while p\==0;   if p//2==1 then _=_*x//n
p=p%2;     x=x*x // n
end    /*while*/
```

return _</lang> output when using the input of: 40 80 180 888
Note the REXX program was executing within a window of 600 bytes wide.

```a= 2988348162058574136915891421498819466320163312926952423791023078876139
(70 digits)
b= 2351399303373464486466122544523690094744975233415544072992656881240319
(70 digits)
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
a**b (mod 10**40)= 1527229998585248450016808958343740453059
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
a**b (mod 10**80)= 53259517041910225328867076245698908287781527229998585248450016808958343740453059
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
a**b (mod 10**180)= 31857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
a**b (mod 10**888)= 2612849643808365153970307063634422265713972370574889513136845452410856423299436762487557161242604471887885300171829510516527484255607339748359444160694661767131561827274483018385170003434853270016569482853811730383390737793312301323406698998964489388587853627711904603124125798753498716559994462054260496622614506334484689315735068762556447491553489235236807309998697854727791160093566968169527719659307289405305177993299425901141782840092602984267350865792542825912897568403588118221513074793528568569833937153488707152390200379629380198479929609788498528506130631774711751914442
62586321233906926671000476591123695550566585083205841790404069511972417770392822283604206143472509425391114072344402850867571806031857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059
```

## 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.

Library: OpenSSL

<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. (Program needs Ruby 1.8.7 for Integer#odd?.)

Works with: Ruby version 1.8.7

<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.odd? and 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>

## Scala

<lang scala>import scala.math.BigInt

val a = BigInt(

``` "2988348162058574136915891421498819466320163312926952423791023078876139")
```

val b = BigInt(

``` "2351399303373464486466122544523690094744975233415544072992656881240319")
```

println(a.modPow(b, BigInt(10).pow(40)))</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

1. Algorithm from http://introcs.cs.princeton.edu/java/78crypto/ModExp.java.html
2. but Tcl has arbitrary-width integers and an exponentiation operator, which
3. 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} {

```   for {set c 1} {\$b} {set a [expr {\$a*\$a % \$n}]} {
```

if {\$b & 1} { set c [expr {\$c*\$a % \$n}] } set b [expr {\$b >> 1}]

```   }
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"```