# Modular exponentiation

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

Modular exponentiation
You are encouraged to solve this task according to the task description, using any language you may know.
•   ${\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}$.

## 11l

Translation of: D
```F pow_mod(BigInt =base, BigInt =exponent, BigInt modulus)
BigInt result = 1

L exponent != 0
I exponent % 2 != 0
result = (result * base) % modulus
exponent I/= 2
base = (base * base) % modulus

R result

print(pow_mod(BigInt(‘2988348162058574136915891421498819466320163312926952423791023078876139’),
BigInt(‘2351399303373464486466122544523690094744975233415544072992656881240319’),
BigInt(10) ^ 40))```
Output:
```1527229998585248450016808958343740453059
```

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

## Arturo

```a: 2988348162058574136915891421498819466320163312926952423791023078876139
b: 2351399303373464486466122544523690094744975233415544072992656881240319

loop [40 80 180 888] 'm ->
print ["(a ^ b) % 10 ^" m "=" powmod a b 10^m]
```
Output:
```(a ^ b) % 10 ^ 40 = 1527229998585248450016808958343740453059
(a ^ b) % 10 ^ 80 = 53259517041910225328867076245698908287781527229998585248450016808958343740453059
(a ^ b) % 10 ^ 180 = 31857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059
(a ^ b) % 10 ^ 888 = 261284964380836515397030706363442226571397237057488951313684545241085642329943676248755716124260447188788530017182951051652748425560733974835944416069466176713156182727448301838517000343485327001656948285381173038339073779331230132340669899896448938858785362771190460312412579875349871655999446205426049662261450633448468931573506876255644749155348923523680730999869785472779116009356696816952771965930728940530517799329942590114178284009260298426735086579254282591289756840358811822151307479352856856983393715348870715239020037962938019847992960978849852850613063177471175191444262586321233906926671000476591123695550566585083205841790404069511972417770392822283604206143472509425391114072344402850867571806031857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059```

## ATS

Library: ats2-xprelude
Library: GMP

For its multiple precision support, you will need ats2-xprelude with GNU MP support enabled.

There is GNU MP support that comes with the ATS2 compiler, however, so perhaps someone will write a demonstration using that. Unlike ats2-xprelude (which assumes there is garbage collection), that support represents multi-precision numbers as linear types.

```(* You will need
https://sourceforge.net/p/chemoelectric/ats2-xprelude/ *)

#include "xprelude/HATS/xprelude.hats"

val a = exrat_make_string_exn "2988348162058574136915891421498819466320163312926952423791023078876139"
val b = exrat_make_string_exn "2351399303373464486466122544523690094744975233415544072992656881240319"

val modulus = exrat_make (10, 1) ** 40

(* xprelude/SATS/exrat.sats includes the "exrat_numerator_modular_pow"
function, based on GMP's mpz_powm. *)
val result1 = exrat_numerator_modular_pow (a, b, modulus)

(* But that was too easy. Here is the right-to-left binary method,
https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
*)
val result2 =
(lam (base     : exrat,
exponent : exrat,
modulus  : exrat) : exrat =>
let
val zero = exrat_make (0, 1)
and one = exrat_make (1, 1)
and two = exrat_make (2, 1)
macdef divrem = exrat_numerator_euclid_division
macdef rem = exrat_numerator_euclid_remainder
in
if modulus = one then
zero
else
let
fun
loop (result   : exrat,
base     : exrat,
exponent : exrat) : exrat =
if iseqz exponent then
result
else
let
val @(exponent, remainder) = exponent \divrem two
val result =
if remainder = one then
(result * base) \rem modulus
else
result
val base = (base * base) \rem modulus
in
loop (result, base, exponent)
end
in
loop (one, base \rem modulus, exponent)
end
end) (a, b, modulus)

implement
main0 () =
begin
println! result1;
println! result2
end```
Output:
```\$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW \$(pkg-config --cflags ats2-xprelude) \$(pkg-config --variable=PATSCCFLAGS ats2-xprelude) modular-exponentiation.dats \$(pkg-config --libs ats2-xprelude) -lgc && ./a.out
1527229998585248450016808958343740453059
1527229998585248450016808958343740453059
```

## AutoHotkey

Library: MPL
```#NoEnv
#SingleInstance, Force
SetBatchLines, -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`

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

## bc

```define p(n, e, m) {
auto r
for (r = 1; e > 0; e /= 2) {
if (e % 2 == 1) r = n * r % m
n = n * n % m
}
return(r)
}

a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
p(a, b, 10 ^ 40)
```
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`

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

## C++

Library: Boost
```#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/integer.hpp>

int main() {
using boost::multiprecision::cpp_int;
using boost::multiprecision::pow;
using boost::multiprecision::powm;
cpp_int a("2988348162058574136915891421498819466320163312926952423791023078876139");
cpp_int b("2351399303373464486466122544523690094744975233415544072992656881240319");
std::cout << powm(a, b, pow(cpp_int(10), 40)) << '\n';
return 0;
}
```
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)))
```

## Crystal

```require "big"

module Integer
module Powmod

# Compute self**e mod m
def powmod(e, m)
r, b = 1, self.to_big_i
while e > 0
r = (b * r) % m if e.odd?
b = (b * b) % m
e >>= 1
end
r
end
end
end

struct Int; include Integer::Powmod end

a = "2988348162058574136915891421498819466320163312926952423791023078876139".to_big_i
b = "2351399303373464486466122544523690094744975233415544072992656881240319".to_big_i
m = 10.to_big_i ** 40

puts a.powmod(b, m)
```
Output:
`1527229998585248450016808958343740453059`

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

## Delphi

Translation of: C#

Thanks for Rudy Velthuis, BigIntegers library [2].

```program Modular_exponentiation;

{\$APPTYPE CONSOLE}

uses
System.SysUtils,
Velthuis.BigIntegers;

var
a, b, m: BigInteger;

begin
a := BigInteger.Parse('2988348162058574136915891421498819466320163312926952423791023078876139');
b := BigInteger.Parse('2351399303373464486466122544523690094744975233415544072992656881240319');
m := BigInteger.Pow(10, 40);
Writeln(BigInteger.ModPow(a, b, m).ToString);
end.
```
Output:
`1527229998585248450016808958343740453059`

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

## Erlang

```%%% For details of the algorithms used see
%%% https://en.wikipedia.org/wiki/Modular_exponentiation

-module modexp_rosetta.
-export [mod_exp/3,binary_exp/2,test/0].

mod_exp(Base,Exp,Mod) when
is_integer(Base),
is_integer(Exp),
is_integer(Mod),
Base > 0,
Exp > 0,
Mod > 0 ->
binary_exp_mod(Base,Exp,Mod).

binary_exp(Base,Exponent) ->
binary_exp(Base,Exponent,1).
binary_exp(_,0,Result) ->
Result;
binary_exp(Base,Exponent,Acc) ->
binary_exp(Base*Base,Exponent bsr 1,Acc * exp_factor(Base,Exponent)).

binary_exp_mod(Base,Exponent,Mod) ->
binary_exp_mod(Base rem Mod,Exponent,Mod,1).
binary_exp_mod(_,0,_,Result) ->
Result;
binary_exp_mod(Base,Exponent,Mod,Acc) ->
binary_exp_mod((Base*Base) rem Mod,
Exponent bsr 1,Mod,(Acc * exp_factor(Base,Exponent))rem Mod).

exp_factor(_,0) ->
1;
exp_factor(Base,1) ->
Base;
exp_factor(Base,Exponent) ->
exp_factor(Base,Exponent band 1).

test() ->
445 = mod_exp(4,13,497),
%% Rosetta code example:
mod_exp(2988348162058574136915891421498819466320163312926952423791023078876139,
2351399303373464486466122544523690094744975233415544072992656881240319,
binary_exp(10,40)).
```
```34> modexp_rosetta:test().
modexp_rosetta:test().
1527229998585248450016808958343740453059
35>
```

## 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-in
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
10 40 ^
^mod .
```
Output:
```1527229998585248450016808958343740453059
```

## Fortran

```module big_integers             ! Big (but not very big) integers.

!
! A very primitive multiple precision module, using the classical
! algorithms (long multiplication, long division) and a mere 8-bit
! "digit" size. Some procedures might assume integers are in a
! two's-complement representation. This module is good enough for us
!

! NOTE: I assume that iachar and achar do not alter the most
!       significant bit.

use, intrinsic :: iso_fortran_env, only: int16
implicit none
private

public :: big_integer
public :: integer2big
public :: string2big
public :: big2string
public :: big_sgn
public :: big_cmp, big_cmpabs
public :: big_neg, big_abs
public :: big_subabs, big_sub
public :: big_mul      ! One might also include a big_muladd.
public :: big_divrem   ! One could also include big_div and big_rem.
public :: big_pow
public :: operator(+)
public :: operator(-)
public :: operator(*)
public :: operator(**)

type :: big_integer
! The representation is sign-magnitude. The radix is 256, which
! is not speed-efficient, but which seemed relatively easy to
! work with if one were writing in standard Fortran (and assuming
! iachar and achar were "8-bit clean").
logical :: sign = .false.  ! .false. => +sign, .true. => -sign.
character, allocatable :: bytes(:)
end type big_integer

character, parameter :: zero = achar (0)
character, parameter :: one = achar (1)

! An integer type capable of holding an unsigned 8-bit value.
integer, parameter :: bytekind = int16

interface operator(+)
end interface

interface operator(-)
module procedure big_neg
module procedure big_sub
end interface

interface operator(*)
module procedure big_mul
end interface

interface operator(**)
module procedure big_pow
end interface

contains

elemental function logical2byte (bool) result (byte)
logical, intent(in) :: bool
character :: byte
if (bool) then
byte = one
else
byte = zero
end if
end function logical2byte

elemental function logical2i (bool) result (i)
logical, intent(in) :: bool
integer :: i
if (bool) then
i = 1
else
i = 0
end if
end function logical2i

elemental function byte2i (c) result (i)
character, intent(in) :: c
integer :: i
i = iachar (c)
end function byte2i

elemental function i2byte (i) result (c)
integer, intent(in) :: i
character :: c
c = achar (i)
end function i2byte

elemental function byte2bk (c) result (i)
character, intent(in) :: c
integer(bytekind) :: i
i = iachar (c, kind = bytekind)
end function byte2bk

elemental function bk2byte (i) result (c)
integer(bytekind), intent(in) :: i
character :: c
c = achar (i)
end function bk2byte

elemental function bk2i (i) result (j)
integer(bytekind), intent(in) :: i
integer :: j
j = int (i)
end function bk2i

elemental function i2bk (i) result (j)
integer, intent(in) :: i
integer(bytekind) :: j
j = int (iand (i, 255), kind = bytekind)
end function i2bk

! Left shift of the least significant 8 bits of a bytekind integer.
elemental function lshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ishft (ibits (a, 0, 8 - i), i)
end function lshftbk

! Right shift of the least significant 8 bits of a bytekind integer.
elemental function rshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ibits (a, i, 8 - i)
end function rshftbk

! Left shift an integer.
elemental function lshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, i)
end function lshfti

! Right shift an integer.
elemental function rshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, -i)
end function rshfti

function integer2big (i) result (a)
integer, intent(in) :: i
type(big_integer), allocatable :: a

!
! To write a more efficient implementation of this procedure is
! left as an exercise for the reader.
!

character(len = 100) :: buffer

write (buffer, '(I0)') i
a = string2big (trim (buffer))
end function integer2big

function string2big (s) result (a)
character(len = *), intent(in) :: s
type(big_integer), allocatable :: a

integer :: n, i, istart, iend
integer :: digit

if ((s(1:1) == '-') .or. s(1:1) == '+') then
istart = 2
else
istart = 1
end if

iend = len (s)

n = (iend - istart + 2) / 2

allocate (a)
allocate (a%bytes(n))

a%bytes = zero
do i = istart, iend
digit = ichar (s(i:i)) - ichar ('0')
if (digit < 0 .or. 9 < digit) error stop
a = short_multiplication (a, 10)
end do
a%sign = (s(1:1) == '-')
call normalize (a)
end function string2big

function big2string (a) result (s)
type(big_integer), intent(in) :: a
character(len = :), allocatable :: s

type(big_integer), allocatable :: q
integer :: r
integer :: sgn

sgn = big_sgn (a)
if (sgn == 0) then
s = '0'
else
q = a
s = ''
do while (big_sgn (q) /= 0)
call short_division (q, 10, q, r)
s = achar (r + ichar ('0')) // s
end do
if (sgn < 0) s = '-' // s
end if
end function big2string

function big_sgn (a) result (sgn)
type(big_integer), intent(in) :: a
integer :: sgn

integer :: n, i

n = size (a%bytes)
i = 1
sgn = 1234
do while (sgn == 1234)
if (i == n + 1) then
sgn = 0
else if (a%bytes(i) /= zero) then
if (a%sign) then
sgn = -1
else
sgn = 1
end if
else
i = i + 1
end if
end do
end function big_sgn

function big_cmp (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp

if (a%sign) then
if (b%sign) then
cmp = -big_cmpabs (a, b)
else
cmp = -1
end if
else
if (b%sign) then
cmp = 1
else
cmp = big_cmpabs (a, b)
end if
end if
end function big_cmp

function big_cmpabs (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp

integer :: n, i
integer :: ia, ib

cmp = 1234
n = max (size (a%bytes), size (b%bytes))
i = n
do while (cmp == 1234)
if (i == 0) then
cmp = 0
else
ia = byteval (a, i)
ib = byteval (b, i)
if (ia < ib) then
cmp = -1
else if (ia > ib) then
cmp = 1
else
i = i - 1
end if
end if
end do
end function big_cmpabs

function big_neg (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .not. c%sign
end function big_neg

function big_abs (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .false.
end function big_abs

function big_add (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c

logical :: sign

if (a%sign) then
if (b%sign) then      ! a <= 0, b <= 0
sign = .true.
else                  ! a <= 0, b >= 0
c = big_subabs (a, b)
sign = .not. c%sign
end if
else
if (b%sign) then      ! a >= 0, b <= 0
c = big_subabs (a, b)
sign = c%sign
else                  ! a >= 0, b >= 0
sign = .false.
end if
end if
c%sign = sign

function big_sub (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c

logical :: sign

if (a%sign) then
if (b%sign) then      ! a <= 0, b <= 0
c = big_subabs (a, b)
sign = .not. c%sign
else                  ! a <= 0, b >= 0
sign = .true.
end if
else
if (b%sign) then      ! a >= 0, b <= 0
sign = .false.
else                  ! a >= 0, b >= 0
c = big_subabs (a, b)
sign = c%sign
end if
end if
c%sign = sign
end function big_sub

function big_addabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c

! Compute abs(a) + abs(b).

integer :: n, nc, i
logical :: carry
type(big_integer), allocatable :: tmp

n = max (size (a%bytes), size (b%bytes))
nc = n + 1

allocate(tmp)
allocate(tmp%bytes(nc))

call add_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call add_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp

function big_subabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c

! Compute abs(a) - abs(b). The result is signed.

integer :: n, i
logical :: carry
type(big_integer), allocatable :: tmp

n = max (size (a%bytes), size (b%bytes))
allocate(tmp)
allocate(tmp%bytes(n))

if (big_cmpabs (a, b) >= 0) then
tmp%sign = .false.
call sub_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
else
tmp%sign = .true.
call sub_bytes (get_byte (b, 1), get_byte (a, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (b, i), get_byte (a, i), carry, tmp%bytes(i), carry)
end do
end if
call normalize (tmp)
c = tmp
end function big_subabs

function big_mul (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c

!
! This is Knuth, Volume 2, Algorithm 4.3.1M.
!

integer :: na, nb, nc
integer :: i, j
integer :: ia, ib, ic
integer :: carry
type(big_integer), allocatable :: tmp

na = size (a%bytes)
nb = size (b%bytes)
nc = na + nb + 1

allocate (tmp)
allocate (tmp%bytes(nc))

tmp%bytes = zero
j = 1
do j = 1, nb
ib = byte2i (b%bytes(j))
if (ib /= 0) then
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = byte2i (tmp%bytes(i + j - 1))
ic = (ia * ib) + ic + carry
tmp%bytes(i + j - 1) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(na + j) = i2byte (carry)
end if
end do
tmp%sign = (a%sign .neqv. b%sign)
call normalize (tmp)
c = tmp
end function big_mul

subroutine big_divrem (a, b, q, r)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable, intent(inout) :: q, r

!
! Division with a remainder that is never negative. Equivalently,
! this is floor division if the divisor is positive, and ceiling
! division if the divisor is negative.
!
! See Raymond T. Boute, "The Euclidean definition of the functions
! div and mod", ACM Transactions on Programming Languages and
! Systems, Volume 14, Issue 2, pp. 127-144.
! https://doi.org/10.1145/128861.128862
!

call nonnegative_division (a, b, .true., .true., q, r)
if (a%sign) then
if (big_sgn (r) /= 0) then
r = big_sub (big_abs (b), r)
end if
q%sign = .not. b%sign
else
q%sign = b%sign
end if
end subroutine big_divrem

function short_addition (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c

! Compute abs(a) + b.

integer :: na, nc, i
logical :: carry
type(big_integer), allocatable :: tmp

na = size (a%bytes)
nc = na + 1

allocate(tmp)
allocate(tmp%bytes(nc))

call add_bytes (a%bytes(1), i2byte (b), .false., tmp%bytes(1), carry)
do i = 2, na
call add_bytes (a%bytes(i), zero, carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp

function short_multiplication (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c

integer :: i, na, nc
integer :: ia, ic
integer :: carry
type(big_integer), allocatable :: tmp

na = size (a%bytes)
nc = na + 1

allocate (tmp)
allocate (tmp%bytes(nc))

tmp%sign = a%sign
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = (ia * b) + carry
tmp%bytes(i) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(nc) = i2byte (carry)
call normalize (tmp)
c = tmp
end function short_multiplication

! Division without regard to signs.
subroutine nonnegative_division (a, b, want_q, want_r, q, r)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_q, want_r
type(big_integer), intent(inout), allocatable :: q, r

integer :: na, nb
integer :: remainder

na = size (a%bytes)
nb = size (b%bytes)

! It is an error if b has "significant" zero-bytes or is equal to
! zero.
if (b%bytes(nb) == zero) error stop

if (nb == 1) then
if (want_q) then
call short_division (a, byte2i (b%bytes(1)), q, remainder)
else
block
type(big_integer), allocatable :: bit_bucket
call short_division (a, byte2i (b%bytes(1)), bit_bucket, remainder)
end block
end if
if (want_r) then
if (allocated (r)) deallocate (r)
allocate (r)
allocate (r%bytes(1))
r%bytes(1) = i2byte (remainder)
end if
else
if (na >= nb) then
call long_division (a, b, want_q, want_r, q, r)
else
if (want_q) q = string2big ("0")
if (want_r) r = a
end if
end if
end subroutine nonnegative_division

subroutine short_division (a, b, q, r)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), intent(inout), allocatable :: q
integer, intent(inout) :: r

!
! This is Knuth, Volume 2, Exercise 4.3.1.16.
!
! The divisor is assumed to be positive.
!

integer :: n, i
integer :: ia, ib, iq
type(big_integer), allocatable :: tmp

ib = b
n = size (a%bytes)

allocate (tmp)
allocate (tmp%bytes(n))

r = 0
do i = n, 1, -1
ia = (256 * r) + byte2i (a%bytes(i))
iq = ia / ib
r = mod (ia, ib)
tmp%bytes(i) = i2byte (iq)
end do
tmp%sign = a%sign
call normalize (tmp)
q = tmp
end subroutine short_division

subroutine long_division (a, b, want_quotient, want_remainder, quotient, remainder)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_quotient, want_remainder
type(big_integer), intent(inout), allocatable :: quotient
type(big_integer), intent(inout), allocatable :: remainder

!
! This is Knuth, Volume 2, Algorithm 4.3.1D.
!
! We do not deal here with the signs of the inputs and outputs.
!
! It is assumed size(a%bytes) >= size(b%bytes), and that b has no
! leading zero-bytes and is at least two bytes long. If b is one
! byte long and nonzero, use short division.
!

integer :: na, nb, m, n
integer :: num_lz, num_nonlz
integer :: j
integer :: qhat
logical :: carry

!
! We will NOT be working with VERY large numbers, and so it will
! be safe to put temporary storage on the stack. (Note: your
! Fortran might put this storage in a heap instead of the stack.)
!
!    v = b, normalized to put its most significant 1-bit all the
!           way left.
!
!    u = a, shifted left by the same amount as b.
!
!    q = the quotient.
!
! The remainder, although shifted left, will end up in u.
!
integer(bytekind) :: u(0:size (a%bytes) + size (b%bytes))
integer(bytekind) :: v(0:size (b%bytes) - 1)
integer(bytekind) :: q(0:size (a%bytes) - size (b%bytes))

na = size (a%bytes)
nb = size (b%bytes)

n = nb
m = na - nb

! In the most significant byte of the divisor, find the number of
! leading zero bits, and the number of bits after that.
block
integer(bytekind) :: tmp
tmp = byte2bk (b%bytes(n))
num_nonlz = bit_size (tmp) - leadz (tmp)
num_lz = 8 - num_nonlz
end block

call normalize_v (b%bytes) ! Make the most significant bit of v be one.
call normalize_u (a%bytes) ! Shifted by the same amount as v.

! Assure ourselves that the most significant bit of v is a one.
if (.not. btest (v(n - 1), 7)) error stop

do j = m, 0, -1
call calculate_qhat (qhat)
call multiply_and_subtract (carry)
q(j) = i2bk (qhat)
end do

if (want_quotient) then
if (allocated (quotient)) deallocate (quotient)
allocate (quotient)
allocate (quotient%bytes(m + 1))
quotient%bytes = bk2byte (q)
call normalize (quotient)
end if

if (want_remainder) then
if (allocated (remainder)) deallocate (remainder)
allocate (remainder)
allocate (remainder%bytes(n))
call unnormalize_u (remainder%bytes)
call normalize (remainder)
end if

contains

subroutine normalize_v (b_bytes)
character, intent(in) :: b_bytes(n)

!
! Normalize v so its most significant bit is a one. Any
! normalization factor that achieves this goal will suffice; we
! choose 2**num_lz. (Knuth uses (2**32) div (y[n-1] + 1).)
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!

integer :: i
integer(bytekind) :: btmp(0:n - 1)

btmp = byte2bk (b_bytes)

v(0) = lshftbk (btmp(0), num_lz)
do i = 1, n - 1
v(i) = ior (lshftbk (btmp(i), num_lz), &
&      rshftbk (btmp(i - 1), num_nonlz))
end do
end subroutine normalize_v

subroutine normalize_u (a_bytes)
character, intent(in) :: a_bytes(m + n)

!
! Shift a leftwards to get u. Shift by as much as b was shifted
! to get v.
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!

integer :: i
integer(bytekind) :: atmp(0:m + n - 1)

atmp = byte2bk (a_bytes)

u(0) = lshftbk (atmp(0), num_lz)
do i = 1, m + n - 1
u(i) = ior (lshftbk (atmp(i), num_lz), &
&      rshftbk (atmp(i - 1), num_nonlz))
end do
u(m + n) = rshftbk (atmp(m + n - 1), num_nonlz)
end subroutine normalize_u

subroutine unnormalize_u (r_bytes)
character, intent(out) :: r_bytes(n)

!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!

integer :: i
integer(bytekind) :: rtmp(0:n - 1)

do i = 0, n - 1
rtmp(i) = ior (rshftbk (u(i), num_lz), &
&         lshftbk (u(i + 1), num_nonlz))
end do
rtmp(n - 1) = rshftbk (u(n - 1), num_lz)

r_bytes = bk2byte (rtmp)
end subroutine unnormalize_u

subroutine calculate_qhat (qhat)
integer, intent(out) :: qhat

integer :: itmp, rhat

itmp = ior (lshfti (bk2i (u(j + n)), 8), &
&      bk2i (u(j + n - 1)))
qhat = itmp / bk2i (v(n - 1))
rhat = mod (itmp, bk2i (v(n - 1)))
if (rshfti (qhat, 8) /= 0) then
continue
else if (qhat * bk2i (v(n - 2)) &
&     > ior (lshfti (rhat, 8), &
&            bk2i (u(j + n - 2)))) then
continue
else
end if
qhat = qhat - 1
rhat = rhat + bk2i (v(n - 1))
if (rshfti (rhat, 8) == 0) then
end if
end if
end do
end subroutine calculate_qhat

subroutine multiply_and_subtract (carry)
logical, intent(out) :: carry

integer :: i
integer :: qhat_v
integer :: mul_carry, sub_carry
integer :: diff

mul_carry = 0
sub_carry = 0
do i = 0, n
! Multiplication.
qhat_v = mul_carry
if (i /= n) qhat_v = qhat_v + (qhat * bk2i (v(i)))
mul_carry = rshfti (qhat_v, 8)
qhat_v = iand (qhat_v, 255)

! Subtraction.
diff = bk2i (u(j + i)) - qhat_v + sub_carry
sub_carry = -(logical2i (diff < 0)) ! Carry 0 or -1.
u(j + i) = i2bk (diff)
end do
carry = (sub_carry /= 0)
end subroutine multiply_and_subtract

integer :: i, carry, sum

q(j) = q(j) - 1_bytekind
carry = 0
do i = 0, n - 1
sum = bk2i (u(j + i)) + bk2i (v(i)) + carry
carry = ishft (sum, -8)
u(j + i) = i2bk (sum)
end do

end subroutine long_division

function big_pow (a, i) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
type(big_integer), allocatable :: c

type(big_integer), allocatable :: base
integer :: exponent, exponent_halved
integer :: j, last_set

if (i < 0) error stop

if (i == 0) then
c = integer2big (1)
else
last_set = bit_size (i) - leadz (i)
base = a
exponent = i
c = integer2big (1)
do j = 0, last_set - 1
exponent_halved = exponent / 2
if (2 * exponent_halved /= exponent) then
c = c * base
end if
exponent = exponent_halved
base = base * base
end do
end if
end function big_pow

subroutine add_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out

integer :: ia, ib, ic

ia = byte2i (a)
if (carry_in) ia = ia + 1
ib = byte2i (b)
ic = ia + ib
c = i2byte (iand (ic, 255))
carry_out = (ic >= 256)

subroutine sub_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out

integer :: ia, ib, ic

ia = byte2i (a)
ib = byte2i (b)
if (carry_in) ib = ib + 1
ic = ia - ib
carry_out = (ic < 0)
if (carry_out) ic = ic + 256
c = i2byte (iand (ic, 255))
end subroutine sub_bytes

function get_byte (a, i) result (byte)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
character :: byte

if (size (a%bytes) < i) then
byte = zero
else
byte = a%bytes(i)
end if
end function get_byte

function byteval (a, i) result (v)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
integer :: v

if (size (a%bytes) < i) then
v = 0
else
v = byte2i (a%bytes(i))
end if
end function byteval

subroutine normalize (a)
type(big_integer), intent(inout) :: a

logical :: done
integer :: i
character, allocatable :: fewer_bytes(:)

! Shorten to the minimum number of bytes.
i = size (a%bytes)
done = .false.
do while (.not. done)
if (i == 1) then
done = .true.
else if (a%bytes(i) /= zero) then
done = .true.
else
i = i - 1
end if
end do
if (i /= size (a%bytes)) then
allocate (fewer_bytes (i))
fewer_bytes = a%bytes(1:i)
call move_alloc (fewer_bytes, a%bytes)
end if

! If the magnitude is zero, then clear the sign bit.
if (size (a%bytes) == 1) then
if (a%bytes(1) == zero) then
a%sign = .false.
end if
end if
end subroutine normalize

end module big_integers

use, non_intrinsic :: big_integers
implicit none

type(big_integer), allocatable :: zero, one, two
type(big_integer), allocatable :: a, b, modulus

zero = integer2big (0)
one = integer2big (1)
two = integer2big (2)

a = string2big ("2988348162058574136915891421498819466320163312926952423791023078876139")
b = string2big ("2351399303373464486466122544523690094744975233415544072992656881240319")
modulus = string2big ("10") ** 40
write (*,*) big2string (modular_pow (a, b, modulus))

contains

! The right-to-left binary method,
! https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
function modular_pow (base, exponent, modulus) result (retval)
type(big_integer), intent(in) :: base, exponent, modulus
type(big_integer), allocatable :: retval

type(big_integer), allocatable :: bas, expnt, remainder, bit_bucket

if (big_sgn (modulus - one) == 0) then
retval = zero
else
retval = one
bas = base
expnt = exponent
do while (big_sgn (expnt) /= 0)
call big_divrem (expnt, two, expnt, remainder)
if (big_sgn (remainder) /= 0) then
call big_divrem (retval * bas, modulus, bit_bucket, retval)
end if
call big_divrem (bas * bas, modulus, bit_bucket, bas)
end do
end if
end function modular_pow

```
Output:
```\$ gfortran -g -fbounds-check -Wall -Wextra modular_exponentiation.f90 && ./a.out
1527229998585248450016808958343740453059
```

## FreeBASIC

```
'From first principles (No external library)
Function _divide(n1 As String,n2 As String,decimal_places As Integer=10,dpflag As String="s") As String
Dim As String number=n1,divisor=n2
dpflag=Lcase(dpflag)
'For MOD
Dim As Integer modstop
If dpflag="mod" Then
If Len(n1)<Len(n2) Then Return n1
If Len(n1)=Len(n2) Then
If n1<n2 Then Return n1
End If
modstop=Len(n1)-Len(n2)+1
End If
If dpflag<>"mod" Then
If dpflag<>"s"  Then dpflag="raw"
End If
Dim runcount As Integer
'_______  LOOK UP TABLES ______________
Dim Qmod(0 To 19) As Ubyte
Dim bool(0 To 19) As Ubyte
For z As Integer=0 To 19
Qmod(z)=(z Mod 10+48)
bool(z)=(-(10>z))
Next z

'_______ SET THE DECIMAL WHERE IT SHOULD BE AT _______
Dim As String part1,part2
#macro set(decimal)
#macro insert(s,char,position)
If position > 0 And position <=Len(s) Then
part1=Mid(s,1,position-1)
part2=Mid(s,position)
s=part1+char+part2
End If
#endmacro
If dpflag="raw" Then
End If
#endmacro
'______________________________________________
'__________ SPLIT A STRING ABOUT A CHARACTRR __________
Dim As String var1,var2
Dim pst As Integer
#macro split(stri,char,var1,var2)
pst=Instr(stri,char)
var1="":var2=""
If pst<>0 Then
var1=Rtrim(Mid(stri,1,pst),".")
var2=Ltrim(Mid(stri,pst),".")
Else
var1=stri
End If
#endmacro

#macro Removepoint(s)
split(s,".",var1,var2)
#endmacro
'__________ GET THE SIGN AND CLEAR THE -ve __________________
Dim sign As String
If Left(number,1)="-" Xor Left (divisor,1)="-" Then sign="-"
If Left(number,1)="-" Then  number=Ltrim(number,"-")
If Left (divisor,1)="-" Then divisor=Ltrim(divisor,"-")

'DETERMINE THE DECIMAL POSITION BEFORE THE DIVISION
Dim As Integer lennint,lenddec,lend,lenn,difflen
split(number,".",var1,var2)
lennint=Len(var1)
split(divisor,".",var1,var2)
lenddec=Len(var2)

If Instr(number,".") Then
Removepoint(number)
number=var1+var2
End If
If Instr(divisor,".") Then
Removepoint(divisor)
divisor=var1+var2
End If
Dim As Integer numzeros
numzeros=Len(number)
number=Ltrim(number,"0"):divisor=Ltrim (divisor,"0")
numzeros=numzeros-Len(number)
lend=Len(divisor):lenn=Len(number)
If lend>lenn Then difflen=lend-lenn
Dim decpos As Integer=lenddec+lennint-lend+2-numzeros 'THE POSITION INDICATOR
Dim _sgn As Byte=-Sgn(decpos)
If _sgn=0 Then _sgn=1
Dim As String thepoint=String(_sgn,".") 'DECIMAL AT START (IF)
Dim As String zeros=String(-decpos+1,"0")'ZEROS AT START (IF) e.g. .0009
If dpflag<>"mod" Then
If Len(zeros) =0 Then dpflag="s"
End If
Dim As Integer runlength
If Len(zeros) Then
runlength=decimal_places
If dpflag="raw" Then
runlength=1
If decimal_places>Len(zeros) Then
runlength=runlength+(decimal_places-Len(zeros))
End If
End If

Else
decimal_places=decimal_places+decpos
runlength=decimal_places
End If
'___________DECIMAL POSITION DETERMINED  _____________

'SET UP THE VARIABLES AND START UP CONDITIONS
number=number+String(difflen+decimal_places,"0")
Dim count As Integer
Dim temp As String
Dim copytemp As String
Dim topstring As String
Dim copytopstring As String
Dim As Integer lenf,lens
Dim As Ubyte takeaway,subtractcarry
Dim As Integer n3,diff
If Ltrim(divisor,"0")="" Then Return "Error :division by zero"
lens=Len(divisor)
topstring=Left(number,lend)
copytopstring=topstring
Do
count=0
Do
count=count+1
copytemp=temp

Do
'___________________ QUICK SUBTRACTION loop _________________

lenf=Len(topstring)
If  lens<lenf=0 Then 'not
If Lens>lenf Then
temp= "done"
Exit Do
End If
If divisor>topstring Then
temp= "done"
Exit Do
End If
End If

diff=lenf-lens
temp=topstring
subtractcarry=0

For n3=lenf-1 To diff Step -1
takeaway= topstring[n3]-divisor[n3-diff]+10-subtractcarry
temp[n3]=Qmod(takeaway)
subtractcarry=bool(takeaway)
Next n3
If subtractcarry=0 Then Exit Do
If n3=-1 Then Exit Do
For n3=n3 To 0 Step -1
takeaway= topstring[n3]-38-subtractcarry
temp[n3]=Qmod(takeaway)
subtractcarry=bool(takeaway)
If subtractcarry=0 Then Exit Do
Next n3
Exit Do

Loop 'single run
temp=Ltrim(temp,"0")
If temp="" Then temp= "0"
topstring=temp
Loop Until temp="done"
' INDIVIDUAL CHARACTERS CARVED OFF ________________
runcount=runcount+1
If count=1 Then
topstring=copytopstring+Mid(number,lend+runcount,1)
Else
topstring=copytemp+Mid(number,lend+runcount,1)
End If
copytopstring=topstring
topstring=Ltrim(topstring,"0")
If dpflag="mod" Then
If runcount=modstop Then
If topstring="" Then Return "0"
Return Mid(topstring,1,Len(topstring)-1)
End If
End If
If topstring="" And runcount>Len(n1)+1 Then
Exit Do
End If
Loop Until runcount=runlength+1

' END OF RUN TO REQUIRED DECIMAL PLACES
set(decimal) 'PUT IN THE DECIMAL POINT
'THERE IS ALWAYS A DECIMAL POINT SOMEWHERE IN THE ANSWER
'NOW GET RID OF IT IF IT IS REDUNDANT
End Function

Dim Shared As Integer _Mod(0 To 99),_Div(0 To 99)
For z As Integer=0 To 99:_Mod(z)=(z Mod 10+48):_Div(z)=z\10:Next

Function Qmult(a As String,b As String) As String
Var flag=0,la = Len(a),lb = Len(b)
If Len(b)>Len(a) Then flag=1:Swap a,b:Swap la,lb
Dim As Ubyte n,carry,ai
Var c =String(la+lb,"0")
For i As Integer =la-1 To 0 Step -1
carry=0:ai=a[i]-48
For j As Integer =lb-1 To 0 Step -1
n = ai * (b[j]-48) + (c[i+j+1]-48) + carry
carry =_Div(n):c[i+j+1]=_Mod(n)
Next j
c[i]+=carry
Next i
If flag Then Swap a,b
Return Ltrim(c,"0")
End Function
'=======================================================================
#define mod_(a,b) _divide((a),(b),,"mod")
#define div_(a,b) iif(len((a))<len((b)),"0",_divide((a),(b),-2))

Function Modular_Exponentiation(base_num As String, exponent As String, modulus As String) As String
Dim b1 As String = base_num
Dim e1 As String = exponent
Dim As String result="1"
b1 = mod_(b1,modulus)
Do While e1 <> "0"
Var L=Len(e1)-1
If e1[L] And 1 Then
result=mod_(Qmult(result,b1),modulus)
End If
b1=mod_(qmult(b1,b1),modulus)
e1=div_(e1,"2")
Loop
Return result
End Function

var base_num="2988348162058574136915891421498819466320163312926952423791023078876139"
var exponent="2351399303373464486466122544523690094744975233415544072992656881240319"
var modulus="10000000000000000000000000000000000000000"
dim as double t=timer
var ans=Modular_Exponentiation(base_num,exponent,modulus)
print "Result:"
Print ans
print "time taken  ";(timer-t)*1000;" milliseconds"
Print "Done"
Sleep```
```Result:
1527229998585248450016808958343740453059
time taken   93.09767815284431 milliseconds
Done

```

## Frink

```a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
println[modPow[a,b,10^40]]```
Output:
```1527229998585248450016808958343740453059
```

## GAP

```# Built-in
a := 2988348162058574136915891421498819466320163312926952423791023078876139;
b := 2351399303373464486466122544523690094744975233415544072992656881240319;
PowerModInt(a, b, 10^40);
1527229998585248450016808958343740453059

# Implementation
PowerModAlt := 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)
# Usage
print 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`

```modPow :: Integer -> Integer -> Integer -> Integer -> Integer
modPow b e 1 r = 0
modPow b 0 m r = r
modPow b e m r
| e `mod` 2 == 1 = modPow b' e' m (r * b `mod` m)
| otherwise = modPow b' e' m r
where
b' = b * b `mod` m
e' = e `div` 2

main = do
print (modPow 2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
(10 ^ 40)
1)
```
Output:
`1527229998585248450016808958343740453059`

or in terms of until:

```powerMod :: Integer -> Integer -> Integer -> Integer
powerMod b e m = x
where
(_, _, x) =
until
(\(_, e, _) -> e <= 0)
(\(b, e, x) ->
( mod (b * b) m
, div e 2
, if 0 /= mod e 2
then mod (b * x) m
else x))
(b, e, 1)

main :: IO ()
main =
print \$
powerMod
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
(10 ^ 40)
```
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
result := 1
while exponent > 0 do {
if exponent % 2 = 1 then
result := (result * base) % modulus
exponent /:= 2
base := base ^ 2 % modulus
}
return result
end
```
Output:
`last 40 digits = 1527229998585248450016808958343740453059`

## J

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

a m&|@^ b
1527229998585248450016808958343740453059
```

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`

## jq

Works with: jq

Also works with gojq, the Go implementation of jq, and with fq.

```# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power(\$b): . as \$in | reduce range(0;\$b) as \$i (1; . * \$in);

# Returns (. ^ \$exp) % \$mod
# where \$exp >= 0, \$mod != 0, and the input are integers.
def modPow(\$exp; \$mod):
if \$mod == 0 then "Cannot take modPow with modulus 0." | error
elif \$exp < 0 then "modPow with exp < 0 is not supported." | error
else . as \$x
| {r: 1, base: (\$x % \$mod), exp: \$exp}
| until( .exp <= 0 or .emit;
if .base == 0 then .emit = 0
else if .exp%2 == 1
then .r = (.r * .base) % \$mod
|    .exp |= (. - 1) / 2
else .exp /= 2
end
| .base |= (. * .) % \$mod
end )
| if .emit then .emit else .r end
end;

2988348162058574136915891421498819466320163312926952423791023078876139 as \$a
| 2351399303373464486466122544523690094744975233415544072992656881240319 as \$b
| (10|power(40)) as \$m
| \$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 = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 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
```

## Lambdatalk

Following scheme

```We will call the lib_BN library for big numbers:

{require lib_BN}

In this library {BN.compare a b} returns 1 if a > b, 0 if a = b and -1 if a < b.
For a better readability we define three small functions

{def BN.= {lambda {:a :b} {= {BN.compare :a :b} 0}}}
-> BN.=
{def BN.even? {lambda {:n} {= {BN.compare {BN.% :n 2} 0} 0}}}
-> BN.even?
{def BN.square {lambda {:n} {BN.* :n :n}}}
-> BN.square

{def mod-exp
{lambda {:a :n :mod}
{if {BN.= :n 0}
then 1
else {if {BN.even? :n}
then {BN.% {BN.square {mod-exp :a {BN./ :n 2} :mod}} :mod}
else {BN.% {BN.* :a   {mod-exp :a {BN.- :n 1} :mod}} :mod}}}}}
-> mod-exp

{mod-exp
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
{BN.pow 10 40}}
-> 1527229998585248450016808958343740453059   // 3300ms
```

## Maple

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

## Mathematica/Wolfram Language

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

## ObjectIcon

Translation of: Icon and Unicon
```# -*- ObjectIcon -*-
#
# This program is close to being an exact copy of the Icon.
#

import io # <-- Object Icon requires this for I/O.

procedure main()
local a, b # <-- Object Icon forces you to declare your variables.

a := 2988348162058574136915891421498819466320163312926952423791023078876139
b := 2351399303373464486466122544523690094744975233415544072992656881240319

# You could leave out the "io." in the call to "write" below,
# because there is some "compatibility with regular Icon" support in
# the io package.
io.write("last 40 digits = ", mod_power(a,b,(10^40)))
end

procedure mod_power(base, exponent, modulus)
local result

result := 1
while exponent > 0 do
{
if exponent % 2 = 1 then
result := (result * base) % modulus
exponent /:= 2
base := base ^ 2 % modulus
}
return result
end```
```\$ oiscript modular-exponentiation-OI.icn
last 40 digits = 1527229998585248450016808958343740453059
```

## OCaml

Using the zarith library:

```let a = Z.of_string "2988348162058574136915891421498819466320163312926952423791023078876139" in
let b = Z.of_string "2351399303373464486466122544523690094744975233415544072992656881240319" in
let m = Z.pow (Z.of_int 10) 40 in
Z.powm a b m
|> Z.to_string
|> print_endline
```
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
```

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

## Phix

Library: Phix/mpfr

Already builtin as mpz_powm, but here is an explicit version

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

## Powershell

Works with: Powershell version 7
```Function Invoke-ModuloExponentiation ([BigInt]\$Base, [BigInt]\$Exponent, \$Modulo) {
\$Result = 1
\$Base = \$Base % \$Modulo
If (\$Base -eq 0) {return 0}

While (\$Exponent -gt 0) {
If ((\$Exponent -band 1) -eq 1) {\$Result = (\$Result * \$Base) % \$Modulo}
\$Exponent = \$Exponent -shr 1
\$Base = (\$Base * \$Base) % \$Modulo
}
return (\$Result % \$Modulo)
}

\$a = [BigInt]::Parse('2988348162058574136915891421498819466320163312926952423791023078876139')
\$b = [BigInt]::Parse('2351399303373464486466122544523690094744975233415544072992656881240319')
\$m = [BigInt]::Pow(10, 40)

Invoke-ModuloExponentiation -Base \$a -Exponent \$b -Modulo \$m
```
Output:
```1527229998585248450016808958343740453059
```

## Prolog

Works with: SWI Prolog

SWI Prolog has a built-in function named powm for this purpose.

```main:-
A = 2988348162058574136915891421498819466320163312926952423791023078876139,
B = 2351399303373464486466122544523690094744975233415544072992656881240319,
M is 10 ** 40,
P is powm(A, B, M),
writeln(P).
```
Output:
```1527229998585248450016808958343740453059
```

## Python

```a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ** 40
print(pow(a, b, m))
```
Output:
`1527229998585248450016808958343740453059`
```def power_mod(b, e, m):
" Without using builtin function "
x = 1
while e > 0:
b, e, x = (
b * b % m,
e // 2,
b * x % m if e % 2 else x
)

return x

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

## Quackery

```  [ temp put 1 unrot
[ dup while
dup 1 & if
[ unrot tuck *
temp share mod
swap rot ]
1 >>
swap dup *
temp share mod
swap again ]
2drop temp release ] is **mod ( n n n --> n )

2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
10 40 ** **mod echo```
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
```

## Raku

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

## REXX

```/* REXX  Modular exponentiation */

say powerMod(,
2988348162058574136915891421498819466320163312926952423791023078876139,,
2351399303373464486466122544523690094744975233415544072992656881240319,,
1e40)
return

powerMod: procedure
parse arg base, exponent, modulus

/* we need a numeric precision of twice the modulus size,    */
/* the exponent size, or the base size, whichever is largest */
numeric digits max(2 * length(format(modulus, , , 0)),,
length(format(exponent, , , 0)), length(format(base, , , 0)))

result = 1
base = base // modulus
do while exponent > 0
if exponent // 2 = 1 then
result = result * base // modulus
base = base * base // modulus
exponent = exponent % 2
end
return result
```
Output:
```1527229998585248450016808958343740453059
```

## Ruby

### Built in since version 2.5.

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

### Using OpenSSL standard library

Library: OpenSSL
```require 'openssl'
a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ** 40
puts a.to_bn.mod_exp(b, m)
```

### Written in Ruby

```def mod_exp(n, e, mod)
fail ArgumentError, 'negative exponent' if e < 0
prod = 1
base = n % mod
until e.zero?
prod = (prod * base) % mod if e.odd?
e >>= 1
base = (base * base) % mod
end
prod
end
```

## 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: [3]

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

## Swift

AttaSwift's BigInt has a built-in modPow method, but here we define a generic modPow.

```import BigInt

func modPow<T: BinaryInteger>(n: T, e: T, m: T) -> T {
guard e != 0 else {
return 1
}

var res = T(1)
var base = n % m
var exp = e

while true {
if exp & 1 == 1 {
res *= base
res %= m
}

if exp == 1 {
return res
}

exp /= 2
base *= base
base %= m
}
}

let a = BigInt("2988348162058574136915891421498819466320163312926952423791023078876139")
let b = BigInt("2351399303373464486466122544523690094744975233415544072992656881240319")

print(modPow(n: a, e: b, m: BigInt(10).power(40)))
```
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 2988348162058574136915891421498819466320163312926952423791023078876139
set b 2351399303373464486466122544523690094744975233415544072992656881240319
set n [expr {10**40}]
puts [expr {modexp(\$a,\$b,\$n)}]
```
Output:
``` 1527229998585248450016808958343740453059
```

### Iterative

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

Demonstrating:

```set a 2988348162058574136915891421498819466320163312926952423791023078876139
set b 2351399303373464486466122544523690094744975233415544072992656881240319
set 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
```

## Wren

Library: Wren-big
```import "./big" for BigInt

var a = BigInt.new("2988348162058574136915891421498819466320163312926952423791023078876139")
var b = BigInt.new("2351399303373464486466122544523690094744975233415544072992656881240319")
var m = BigInt.ten.pow(40)
System.print(a.modPow(b, m))
```
Output:
```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
```

## zsh

```#!/bin/zsh

#
# I use expr from GNU coreutils.
EXPR=`which expr`
#
# One could use GNU bc or any other such program that is capable of
# handling the large integers.
#

unset PATH          # Besides shell commands, the script uses ONLY
# the expression calculator.

mod_power()
{
local base="\${1}"
local exponent="\${2}"
local modulus="\${3}"
local result=1

while [[ `\${EXPR} \${exponent} '>' 0` != 0 ]] ; do
if [[ `\${EXPR} '(' \${exponent} '%' 2 ')' '=' 1` != 0 ]] ; then
result=`\${EXPR} '(' \${result} '*' \${base} ')' '%' \${modulus}`
fi
exponent=`\${EXPR} \${exponent} / 2`
base=`\${EXPR} '(' \${base} '*' \${base} ')' '%' \${modulus}`
done
echo \${result}
}

a=2988348162058574136915891421498819466320163312926952423791023078876139
b=2351399303373464486466122544523690094744975233415544072992656881240319

# One followed by 40 zeros.
modulus=10000000000000000000000000000000000000000

# Or the following will work, as long as the modulus is at least
# 10**40. :)))))
#modulus=\$(mod_power 10 40 1000000000000000000000000000000000000000000000000000000000000)

echo \$(mod_power \${a} \${b} \${modulus})

exit 0
```
Output:
```\$ zsh ./modular-exponentiation.zsh
1527229998585248450016808958343740453059```