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

`with Ada.Text_IO, Ada.Command_Line, Crypto.Types.Big_Numbers; 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) = ");   Ada.Text_IO.Put_Line(LN.Utils.To_String(LN.Mod_Utils.Pow((+A), (+B), M)));end Mod_Exp;`
Output:
`A**B (mod 10**40) = 1527229998585248450016808958343740453059`

## ALGOL 68

The code below uses Algol 68 Genie which provides arbitrary precision arithmetic for LONG LONG modes.

` BEGIN   PR precision=1000 PR   MODE LLI = LONG LONG INT;	CO For brevity CO   PROC mod power = (LLI base, exponent, modulus) LLI :   BEGIN      LLI result := 1, b := base, e := exponent;      IF exponent < 0      THEN	 put (stand error, (("Negative exponent", exponent, newline)))      ELSE	 WHILE e > 0 	 DO	    (ODD e | result := (result * b) MOD modulus);	    e OVERAB 2; b := (b * b) MOD modulus	 OD      FI;      result   END;   LLI a = 2988348162058574136915891421498819466320163312926952423791023078876139;   LLI b = 2351399303373464486466122544523690094744975233415544072992656881240319;   LLI m = 10000000000000000000000000000000000000000;   printf ((\$"Last 40 digits = ", 40dl\$, mod power (a, b, m)))END `
Output:
```Last 40 digits = 1527229998585248450016808958343740453059
```

## AutoHotkey

Library: MPL
`#NoEnv#SingleInstance, ForceSetBatchLines, -1#Include mpl.ahk   MP_SET(base, "2988348162058574136915891421498819466320163312926952423791023078876139"), MP_SET(exponent, "2351399303373464486466122544523690094744975233415544072992656881240319"), MP_SET(modulus, "10000000000000000000000000000000000000000") , NumGet(exponent,0,"Int") = -1 ? return : "", MP_SET(result, "1"), MP_SET(TWO, "2")while !MP_IS0(exponent)	MP_DIV(q, r, exponent, TWO)	, (MP_DEC(r) = 1		? (MP_MUL(temp, result, base)		, MP_DIV(q, result, temp, modulus))		: "")	, MP_DIV(q, r, exponent, TWO)	, MP_CPY(exponent, q)	, MP_CPY(base1, base)	, MP_MUL(base2, base1, base)	, MP_DIV(q, base, base2, modulus) msgbox % MP_DEC(result)Return`
Output:
`1527229998585248450016808958343740453059`

## Bracmat

Translation of: Icon_and_Unicon
`  ( ( 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))  )`
Output:
`last 40 digits =  1527229998585248450016808958343740453059`

## BBC BASIC

Uses the Huge Integer Math & Encryption library.

`      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)`
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
`#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;}`

Output:

```1527229998585248450016808958343740453059
```

## C#

We can use the built-in function from BigInteger.

`using System;using System.Numerics; class Program{    static void Main() {        var a = BigInteger.Parse("2988348162058574136915891421498819466320163312926952423791023078876139");        var b = BigInteger.Parse("2351399303373464486466122544523690094744975233415544072992656881240319");        var m = BigInteger.Pow(10, 40);        Console.WriteLine(BigInteger.ModPow(a, b, m));    }}`
Output:
```1527229998585248450016808958343740453059
```

## Clojure

`(defn powerMod "modular exponentiation" [b e m]  (defn m* [p q] (mod (* p q) m))  (loop [b b, e e, x 1]    (if (zero? e) x      (if (even? e) (recur (m* b b) (/ e 2) x)        (recur (m* b b) (quot e 2) (m* b x))))))`
` (defn modpow  " b^e mod m (using Java which solves some cases the pure clojure method has to be modified to tackle--i.e. with large b & e and     calculation simplications when gcd(b, m) == 1 and gcd(e, m) == 1) "  [b e m]  (.modPow (biginteger b) (biginteger e) (biginteger m)))`

## Common 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))))`
Works with: CLISP
`;; CLISP provides EXT:MOD-EXPT(let ((a 2988348162058574136915891421498819466320163312926952423791023078876139)      (b 2351399303373464486466122544523690094744975233415544072992656881240319))  (format t "~A~%" (mod-expt a b (expt 10 40))))`

### Implementation with LOOP

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

## D

Translation of: Icon_and_Unicon

Compile this module with `-version=modular_exponentiation` to see the output.

`module modular_exponentiation; private import std.bigint; BigInt powMod(BigInt base, BigInt exponent, in BigInt modulus)pure nothrow /*@safe*/ in {   assert(exponent >= 0);} body {    BigInt result = 1;     while (exponent) {        if (exponent & 1)            result = (result * base) % modulus;        exponent /= 2;        base = base ^^ 2 % modulus;    }     return result;} version (modular_exponentiation)    void main() {        import std.stdio;         powMod(BigInt("29883481620585741369158914214988194" ~                      "66320163312926952423791023078876139"),               BigInt("235139930337346448646612254452369009" ~                      "4744975233415544072992656881240319"),               BigInt(10) ^^ 40).writeln;    }`
Output:
`1527229998585248450016808958343740453059`

## Dc

`2988348162058574136915891421498819466320163312926952423791023078876139 2351399303373464486466122544523690094744975233415544072992656881240319 10 40^|p`

## EchoLisp

` (lib 'bigint) (define a 2988348162058574136915891421498819466320163312926952423791023078876139)(define b 2351399303373464486466122544523690094744975233415544072992656881240319)(define m 1e40) (powmod a b m)    → 1527229998585248450016808958343740453059 ;; powmod is a native function;; it could be defined as follows : (define (xpowmod base exp mod)    (define result 1)    (while ( !zero? exp)        (when (odd? exp) (set! result (% (* result base) mod)))    (/= exp 2)    (set! base (% (* base base) mod)))result) (xpowmod a b m)    → 1527229998585248450016808958343740453059 `

## Emacs Lisp

Library: Calc
`(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)))`

## F#

`let expMod a b n =    let rec loop a b c =        if b = 0I then c else            loop (a*a%n) (b>>>1) (if b&&&1I = 0I then c else c*a%n)    loop a b 1I [<EntryPoint>]let main argv =    let a = 2988348162058574136915891421498819466320163312926952423791023078876139I    let b = 2351399303373464486466122544523690094744975233415544072992656881240319I    printfn "%A" (expMod a b (10I**40))    // -> 1527229998585248450016808958343740453059    0`

## Factor

`! Built-in2988348162058574136915891421498819466320163312926952423791023078876139235139930337346448646612254452369009474497523341554407299265688124031910 40 ^^mod .`
Output:
```1527229998585248450016808958343740453059
```

## GAP

`# Built-ina := 2988348162058574136915891421498819466320163312926952423791023078876139;b := 2351399303373464486466122544523690094744975233415544072992656881240319;PowerModInt(a, b, 10^40);1527229998585248450016808958343740453059 # ImplementationPowerModAlt := function(a, n, m)    local r;    r := 1;    while n > 0 do        if IsOddInt(n) then            r := RemInt(r*a, m);        fi;        n := QuoInt(n, 2);        a := RemInt(a*a, m);    od;    return r;end; PowerModAlt(a, b, 10^40);`

## gnuplot

`_powm(b, e, m, r) = (e == 0 ? r : (e % 2 == 1 ? _powm(b * b % m, e / 2, m, r * b % m) : _powm(b * b % m, e / 2, m, r)))powm(b, e, m) = _powm(b, e, m, 1)# Usageprint powm(2, 3453, 131)# Where b is the base, e is the exponent, m is the modulus, i.e.: b^e mod m`

## 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)}`
Output:
```1527229998585248450016808958343740453059
```

## Groovy

`println 2988348162058574136915891421498819466320163312926952423791023078876139.modPow(            2351399303373464486466122544523690094744975233415544072992656881240319,            10000000000000000000000000000000000000000)`

Ouput:

`1527229998585248450016808958343740453059`

`powm :: Integer -> Integer -> Integer -> Integer -> Integerpowm b 0 m r = rpowm b e m r  | e `mod` 2 == 1 = powm (b * b `mod` m) (e `div` 2) m (r * b `mod` m)powm b e m r = powm (b * b `mod` m) (e `div` 2) m r main :: IO ()main =  print \$  powm    2988348162058574136915891421498819466320163312926952423791023078876139    2351399303373464486466122544523690094744975233415544072992656881240319    (10 ^ 40)    1`
Output:
`1527229998585248450016808958343740453059`

## Icon and Unicon

This uses the exponentiation procedure from RSA Code an example of the right to left binary method.

`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 resultend`
Output:
`last 40 digits = 1527229998585248450016808958343740453059`

## J

Solution:
`   m&|@^`
Example:
`   a =: 2988348162058574136915891421498819466320163312926952423791023078876139x   b =: 2351399303373464486466122544523690094744975233415544072992656881240319x   m =: 10^40x    a m&|@^ b1527229998585248450016808958343740453059`

Discussion: The phrase a m&|@^ b 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.

`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));    }}`
Output:
`1527229998585248450016808958343740453059`

## Julia

Works with: Julia version 1.0

We can use the built-in `powermod` function with the built-in `BigInt` type (based on GMP):

`a = 2988348162058574136915891421498819466320163312926952423791023078876139b = 2351399303373464486466122544523690094744975233415544072992656881240319m = big(10) ^ 40@show powermod(a, b, m)`
Output:
`powermod(a, b, m) = 1527229998585248450016808958343740453059`

## Kotlin

`// version 1.0.6 import java.math.BigInteger fun main(args: Array<String>) {    val a = BigInteger("2988348162058574136915891421498819466320163312926952423791023078876139")    val b = BigInteger("2351399303373464486466122544523690094744975233415544072992656881240319")    val m = BigInteger.TEN.pow(40)    println(a.modPow(b, m))}`
Output:
```1527229998585248450016808958343740453059
```

## Maple

`a := 2988348162058574136915891421498819466320163312926952423791023078876139:b := 2351399303373464486466122544523690094744975233415544072992656881240319:a &^ b mod 10^40;`
Output:
`1527229998585248450016808958343740453059`

## Mathematica

`a = 2988348162058574136915891421498819466320163312926952423791023078876139;b = 2351399303373464486466122544523690094744975233415544072992656881240319;m = 10^40;PowerMod[a, b, m]-> 1527229998585248450016808958343740453059`

## Maxima

`a: 2988348162058574136915891421498819466320163312926952423791023078876139\$b: 2351399303373464486466122544523690094744975233415544072992656881240319\$power_mod(a, b, 10^40);/* 1527229998585248450016808958343740453059 */`

## Nim

Library: bigints
`import bigints proc powmod(b, e, m: BigInt): BigInt =  assert e >= 0  var e = e  var b = b  result = initBigInt(1)  while e > 0:    if e mod 2 == 1:      result = (result * b) mod m    e = e div 2    b = (b.pow 2) mod m var  a = initBigInt("2988348162058574136915891421498819466320163312926952423791023078876139")  b = initBigInt("2351399303373464486466122544523690094744975233415544072992656881240319") echo powmod(a, b, 10.pow 40)`
Output:
`1527229998585248450016808958343740453059`

## Oforth

Usage : a b mod powmod

`: powmod(base, exponent, modulus)   1 exponent dup ifZero: [ return ]    while ( dup 0 > ) [       dup isEven ifFalse: [ swap base * modulus mod swap ]       2 / base sq modulus mod ->base      ] drop ;`
Output:
```>2988348162058574136915891421498819466320163312926952423791023078876139
ok
>2351399303373464486466122544523690094744975233415544072992656881240319
ok
>10 40 pow
ok
>powmod println
1527229998585248450016808958343740453059
ok
```

## ooRexx

`/* Modular exponentiation */ numeric digits 100say powerMod(, 2988348162058574136915891421498819466320163312926952423791023078876139,, 2351399303373464486466122544523690094744975233415544072992656881240319,, 1e40)exit powerMod: procedure use strict arg base, exponent, modulus exponent=exponent~d2x~x2b~strip('L','0')result=1base = base // modulusdo exponentPos=exponent~length to 1 by -1  if (exponent~subChar(exponentPos) == '1')    then result = (result * base) // modulus  base = (base * base) // modulusendreturn result`
Output:
```1527229998585248450016808958343740453059
```

## PARI/GP

`a=2988348162058574136915891421498819466320163312926952423791023078876139;b=2351399303373464486466122544523690094744975233415544072992656881240319;lift(Mod(a,10^40)^b)`

## Pascal

Works with: Free_Pascal
Library: GMP

A port of the C example using gmp.

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

## Perl

Calculating the result both with an explicit algorithm (borrowed from Perl 6 example) and with a built-in available when the `use bigint` pragma is invoked. Note that `bmodpow` modifies the base value (here `\$a`) while `expmod` does not.

`use bigint; sub expmod {    my(\$a, \$b, \$n) = @_;    my \$c = 1;    do {        (\$c *= \$a) %= \$n if \$b % 2;        (\$a *= \$a) %= \$n;    } while (\$b = int \$b/2);    \$c;} my \$a = 2988348162058574136915891421498819466320163312926952423791023078876139;my \$b = 2351399303373464486466122544523690094744975233415544072992656881240319;my \$m = 10 ** 40; print expmod(\$a, \$b, \$m), "\n";print \$a->bmodpow(\$b, \$m), "\n";`
Output:
```1527229998585248450016808958343740453059
1527229998585248450016808958343740453059```

## Perl 6

This is specced as a built-in, but here's an explicit version:

`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;`
Output:
`1527229998585248450016808958343740453059`

## Phix

One I already had lying around. Note this routine is fairly likely to end up being added to bigatom.e in a future release.

`include bigatom.e function ba_mod_exp(object base, exponent, modulus)-- base/exponent/modulus can be integer/string/bigatom-- returns mod(power(base,exponent),modulus), but in bigatoms and faster.    bigatom res    if ba_compare(exponent,1)=0 then        res = ba_new(base)    else        bool odd = (ba_compare(ba_mod(exponent,2),0)!=0)        if odd then            exponent = ba_sub(exponent,1)        end if        exponent = ba_divide(exponent,2)        res = ba_mod_exp(base,exponent,modulus)        res = ba_multiply(res,res)        if odd then            res = ba_multiply(res,base)        end if    end if    res = ba_mod(res,modulus)    return resend function string s = ba_sprint(ba_mod_exp("2988348162058574136915891421498819466320163312926952423791023078876139",                                "2351399303373464486466122544523690094744975233415544072992656881240319",                                1e40))?s`
Output:
```"1527229998585248450016808958343740453059"
```

## PHP

`<?php\$a = '2988348162058574136915891421498819466320163312926952423791023078876139';\$b = '2351399303373464486466122544523690094744975233415544072992656881240319';\$m = '1' . str_repeat('0', 40);echo bcpowmod(\$a, \$b, \$m), "\n";`
Output:
`1527229998585248450016808958343740453059`

## PicoLisp

The following function is taken from "lib/rsa.l":

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

Test:

`: (**Mod   2988348162058574136915891421498819466320163312926952423791023078876139   2351399303373464486466122544523690094744975233415544072992656881240319   10000000000000000000000000000000000000000 )-> 1527229998585248450016808958343740453059`

## Python

`a = 2988348162058574136915891421498819466320163312926952423791023078876139b = 2351399303373464486466122544523690094744975233415544072992656881240319m = 10 ** 40print(pow(a, b, m))`
Output:
`1527229998585248450016808958343740453059`

## OCaml

Using the zarith library:

` let a = Z.of_string "2988348162058574136915891421498819466320163312926952423791023078876139" inlet b = Z.of_string "2351399303373464486466122544523690094744975233415544072992656881240319" inlet m = Z.pow (Z.of_int 10) 40 inZ.powm a b m|> Z.to_string|> print_endline`
Output:
`1527229998585248450016808958343740453059`

## Racket

` #lang racket(require math)(define a 2988348162058574136915891421498819466320163312926952423791023078876139)(define b 2351399303373464486466122544523690094744975233415544072992656881240319)(define m (expt 10 40))(modular-expt a b m) `
Output:
```1527229998585248450016808958343740453059
```

## REXX

### version 1

This REXX program attempts to handle   any   a, b, or m,   but there are limits for any computer language.
For some REXXes, it's around eight million digits for any arithmetic expression or value, which puts limitations on the
values of   a2   or   10m.

There is REXX code (below) to automatically adjust the number of digits to accommodate huge numbers which are
computed when raising large numbers to some arbitrary power.

`/*REXX program  displays  modular exponentiation of:             a**b  mod  M           */parse arg a b mm                                      /*obtain optional args from the CL*/if a=='' | a==","  then a=2988348162058574136915891421498819466320163312926952423791023078876139if b=='' | b==","  then b=2351399303373464486466122544523690094744975233415544072992656881240319if mm='' | mm=","  then mm=40                         /*MM not specified?   Use default.*/say 'a=' a;   say "        ("length(a) 'digits)'      /*display the  value of  A.       */say 'b=' b;   say "        ("length(b) 'digits)'      /*   "     "     "    "  B.       */      do j=1  for words(mm);   m=word(mm,j)            /*use one of the MM powers (list).*/     say copies('─', linesize()-1)                    /*show a nice separator fence line*/     say 'a**b (mod 10**'m")="   powerMod(a,b,10**m)  /*display the answer ───► console.*/     end   /*j*/exit                                                  /*stick a fork in it, we're done. *//*──────────────────────────────────────────────────────────────────────────────────────*/powerMod: procedure;  parse arg x,p,n                 /*fast modular exponentiation code*/          if p==0  then return 1                      /*special case of  P  being zero. */          if p==1  then return x                      /*   "      "   "  "    "   unity.*/          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  /*obtain the biggest of the three.*/          numeric digits max(20, e*2)                 /*big enough to handle  A².       */          \$=1                                         /*use this for the first value.   */                     do  while p\==0                  /*perform  while   P   isn't zero.*/                     if p//2  then \$=\$ * x  // n      /*is  P  odd?  (is ÷ remainder≡1).*/                     p=p%2;        x=x * x  // n      /*halve  P;   calculate  x² mod n */                     end   /*while*/                  /* [↑]  keep mod'ing 'til equal 0.*/          return \$`

This REXX program makes use of   LINESIZE   REXX program (or BIF) which is used to determine the screen width (or linesize) of the terminal (console).
The   LINESIZE.REX   REXX program is included here   ──►   LINESIZE.REX.

output   when using the input of:     ,   ,   40   80   180   888

Note the REXX program was executing within a window of 600 characters wide, and even with that, the last result wrapped.

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

### version 2

This REXX version handles only up to 100 decimal digits.

`/* Modular exponentiation */  numeric digits 100say powerMod(,  2988348162058574136915891421498819466320163312926952423791023078876139,,  2351399303373464486466122544523690094744975233415544072992656881240319,,  1e40)exit  powerMod: procedure parse arg base, exponent, modulus  exponent = strip(x2b(d2x(exponent)), 'L', '0')result = 1base = base // modulusdo exponentPos=length(exponent) to 1 by -1    if substr(exponent, exponentPos, 1) = '1'        then result = (result * base) // modulus    base = (base * base) // modulusendreturn result`
Output:
```1527229998585248450016808958343740453059
```

## Ruby

Ruby has modular exponentiation built in since version 2.5.

`a = 2988348162058574136915891421498819466320163312926952423791023078876139b = 2351399303373464486466122544523690094744975233415544072992656881240319m = 10**40puts a.pow(b, m)`

Earlier versions required OpenSSL, in Ruby's standard library, which 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
`require 'openssl' a = 2988348162058574136915891421498819466320163312926952423791023078876139b = 2351399303373464486466122544523690094744975233415544072992656881240319m = 10 ** 40puts a.to_bn.mod_exp(b, m)`

## Rust

`/* Add this line to the [dependencies] section of your Cargo.toml file:num = "0.2.0"*/  use num::bigint::BigInt;use num::bigint::ToBigInt;  // The modular_exponentiation() function takes three identical types// (which get cast to BigInt), and returns a BigInt:fn modular_exponentiation<T: ToBigInt>(n: &T, e: &T, m: &T) -> BigInt {    // Convert n, e, and m to BigInt:    let n = n.to_bigint().unwrap();    let e = e.to_bigint().unwrap();    let m = m.to_bigint().unwrap();     // Sanity check:  Verify that the exponent is not negative:    assert!(e >= Zero::zero());     use num::traits::{Zero, One};     // As most modular exponentiations do, return 1 if the exponent is 0:    if e == Zero::zero() {        return One::one()    }     // Now do the modular exponentiation algorithm:    let mut result: BigInt = One::one();    let mut base = n % &m;    let mut exp = e;     // Loop until we can return out result:    loop {        if &exp % 2 == One::one() {            result *= &base;            result %= &m;        }         if exp == One::one() {            return result        }         exp /= 2;        base *= base.clone();        base %= &m;    }}`

Test code:

`fn main() {    let (a, b, num_digits) = (  "2988348162058574136915891421498819466320163312926952423791023078876139",  "2351399303373464486466122544523690094744975233415544072992656881240319",  "40",                    );     // Covert a, b, and num_digits to numbers:    let a = BigInt::parse_bytes(a.as_bytes(), 10).unwrap();    let b = BigInt::parse_bytes(b.as_bytes(), 10).unwrap();    let num_digits = num_digits.parse().unwrap();     // Calculate m from num_digits:    let m = num::pow::pow(10.to_bigint().unwrap(), num_digits);     // Get the result and print it:    let result = modular_exponentiation(&a, &b, &m);    println!("The last {} digits of\n{}\nto the power of\n{}\nare:\n{}",             num_digits, a, b, result);}`
Output:
```The last 40 digits of
2988348162058574136915891421498819466320163312926952423791023078876139
to the power of
2351399303373464486466122544523690094744975233415544072992656881240319
are:
1527229998585248450016808958343740453059```

## Scala

`import scala.math.BigInt val a = BigInt(  "2988348162058574136915891421498819466320163312926952423791023078876139")val b = BigInt(  "2351399303373464486466122544523690094744975233415544072992656881240319") println(a.modPow(b, BigInt(10).pow(40)))`

## Scheme

` (define (square n)  (* n n)) (define (mod-exp a n mod)  (cond ((= n 0) 1)        ((even? n)          (remainder (square (mod-exp a (/ n 2) mod))                     mod))        (else (remainder (* a (mod-exp a (- n 1) mod))                          mod)))) (define result  (mod-exp 2988348162058574136915891421498819466320163312926952423791023078876139            2351399303373464486466122544523690094744975233415544072992656881240319            (expt 10 40)))`
Output:
```> result
1527229998585248450016808958343740453059
```

## Seed7

The library bigint.s7i defines the function modPow, which does modular exponentiation.

`\$ include "seed7_05.s7i";  include "bigint.s7i"; const proc: main is func  begin    writeln(modPow(2988348162058574136915891421498819466320163312926952423791023078876139_,                   2351399303373464486466122544523690094744975233415544072992656881240319_,                   10_ ** 40));  end func;`
Output:
```1527229998585248450016808958343740453059
```

The library bigint.s7i defines modPow with:

`const func bigInteger: modPow (in var bigInteger: base,    in var bigInteger: exponent, in bigInteger: modulus) is func  result    var bigInteger: power is 1_;  begin    if exponent < 0_ or modulus < 0_ then      raise RANGE_ERROR;    else      while exponent > 0_ do        if odd(exponent) then          power := (power * base) mod modulus;        end if;        exponent >>:= 1;        base := base ** 2 mod modulus;      end while;    end if;  end func;`

Original source: [2]

## Sidef

Built-in:

`say expmod(    2988348162058574136915891421498819466320163312926952423791023078876139,    2351399303373464486466122544523690094744975233415544072992656881240319,    10**40)`

User-defined:

`func expmod(a, b, n) {    var c = 1    do {        (c *= a) %= n if b.is_odd        (a *= a) %= n    } while (b //= 2)    c}`
Output:
```1527229998585248450016808958343740453059
```

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

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

Demonstrating:

`set a 2988348162058574136915891421498819466320163312926952423791023078876139set b 2351399303373464486466122544523690094744975233415544072992656881240319set n [expr {10**40}]puts [expr {modexp(\$a,\$b,\$n)}]`
Output:
``` 1527229998585248450016808958343740453059
```

### Iterative

`package require Tcl 8.5proc 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 }`

Demonstrating:

`set a 2988348162058574136915891421498819466320163312926952423791023078876139set b 2351399303373464486466122544523690094744975233415544072992656881240319set n [expr {10**40}]puts [modexp \$a \$b \$n]`
Output:
``` 1527229998585248450016808958343740453059
```

## TXR

`\$ txr -p '(exptmod 2988348162058574136915891421498819466320163312926952423791023078876139                   2351399303373464486466122544523690094744975233415544072992656881240319                   (expt 10 40)))'1527229998585248450016808958343740453059`

## Visual Basic .NET

Works with: Visual Basic .NET version 2011
`' Modular exponentiation - VB.Net - 21/01/2019    Imports System.Numerics     Private Sub Main()        Dim a, b, m, x As BigInteger        a = BigInteger.Parse("2988348162058574136915891421498819466320163312926952423791023078876139")        b = BigInteger.Parse("2351399303373464486466122544523690094744975233415544072992656881240319")        m = BigInteger.Pow(10, 40)   '=10^40        x = ModPowBig(a, b, m)        Debug.Print("x=" & x.ToString)    End Sub 'Main     Function ModPowBig(ByVal base As BigInteger, ByVal exponent As BigInteger, ByVal modulus As BigInteger) As BigInteger        Dim result As BigInteger        result = 1        Do While exponent > 0            If (exponent Mod 2) = 1 Then                result = (result * base) Mod modulus            End If            exponent = exponent / 2            base = (base * base) Mod modulus        Loop        Return result    End Function 'ModPowBig`
Output:
```x=1527229998585248450016808958343740453059
```

## zkl

Using the GMP big num library:

`var BN=Import("zklBigNum");a:=BN("2988348162058574136915891421498819466320163312926952423791023078876139");b:=BN("2351399303373464486466122544523690094744975233415544072992656881240319");m:=BN(10).pow(40);a.powm(b,m).println();a.powm(b,m) : "%,d".fmt(_).println();`
Output:
```1527229998585248450016808958343740453059
1,527,229,998,585,248,450,016,808,958,343,740,453,059
```