Modular exponentiation

From Rosetta Code
Task
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   ,   where

  •  
  •  


A computer is too slow to find the entire value of   .

Instead, the program must use a fast algorithm for modular exponentiation:   .

The algorithm must work for any integers   ,     where     and   .

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

Ada

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

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 "share/atspre_staload.hats"

#include "xprelude/HATS/xprelude.hats"
staload "xprelude/SATS/exrat.sats"
staload _ = "xprelude/DATS/exrat.dats"

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);
  readln;
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
  ! to fulfill the task.
  !

  ! 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_addabs, big_add
  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(+)
     module procedure big_add
  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)
       a = short_addition (a, digit)
    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
          c = big_addabs (a, b)
          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
          c = big_addabs (a, b)
          sign = .false.
       end if
    end if
    c%sign = sign
  end function big_add

  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
          c = big_addabs (a, b)
          sign = .true.
       end if
    else
       if (b%sign) then      ! a >= 0, b <= 0
          c = big_addabs (a, b)
          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
  end function big_addabs

  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
          q = short_addition (q, 1)
          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
  end function short_addition

  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)
       if (carry) call add_back
    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
      logical :: adjust

      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)))
      adjust = .true.
      do while (adjust)
         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
            adjust = .false.
         end if
         if (adjust) then
            qhat = qhat - 1
            rhat = rhat + bk2i (v(n - 1))
            if (rshfti (rhat, 8) == 0) then
               adjust = .false.
            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

    subroutine add_back
      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 add_back

  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)
  end subroutine add_bytes

  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

program modular_exponentiation_task
  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

end program modular_exponentiation_task
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
    Dim answer As String   'THE ANSWER STRING  
    
    '_______ 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
    insert(answer,".",decpos)
    answer=thepoint+zeros+answer
    If dpflag="raw" Then
        answer=Mid(answer,1,decimal_places)
    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
        answer=String(Len(zeros)+runlength+10,"0")
        If dpflag="raw" Then 
            runlength=1
            answer=String(Len(zeros)+runlength+10,"0")
            If decimal_places>Len(zeros) Then
                runlength=runlength+(decimal_places-Len(zeros))
                answer=String(Len(zeros)+runlength+10,"0")
            End If
        End If
        
    Else
        decimal_places=decimal_places+decpos
        runlength=decimal_places
        answer=String(Len(zeros)+runlength+10,"0")
    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
        answer[runcount-1]=count+47
        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
    answer=Rtrim(answer,"0")
    answer=Rtrim(answer,".")
    answer=Ltrim(answer,"0")
    If answer="" Then Return "0"
    Return sign+answer
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

Haskell

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 
   if exponent < 0 then runerr(205,m)          # added for this task
   result := 1
   while exponent > 0 do {
      if exponent % 2 = 1 then 
         result := (result * base) % modulus
      exponent /:= 2   
      base := base ^ 2 % modulus
      }  
   return result
end
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.

Adapted from Wren

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

def task:
    2988348162058574136915891421498819466320163312926952423791023078876139 as $a
  | 2351399303373464486466122544523690094744975233415544072992656881240319 as $b
  | (10|power(40)) as $m
  | $a | modPow($b; $m) ;

task
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

From your system prompt:

$ 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