Modular exponentiation: Difference between revisions

Content added Content deleted

Inline

Revision as of 13:22, 26 July 2021

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

$a=2988348162058574136915891421498819466320163312926952423791023078876139$
$b=2351399303373464486466122544523690094744975233415544072992656881240319$

A computer is too slow to find the entire value of $a^{b}$ .

Instead, the program must use a fast algorithm for modular exponentiation: $a^{b}\mod m$ .

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

11l

Translation of: D

<lang 11l>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))</lang>

Output:

1527229998585248450016808958343740453059

Ada

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

<lang Ada>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;</lang>

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.

<lang algol68> 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 </lang>

Output:

Last 40 digits = 1527229998585248450016808958343740453059

Arturo

<lang rebol>a: 2988348162058574136915891421498819466320163312926952423791023078876139 b: 2351399303373464486466122544523690094744975233415544072992656881240319

loop [40 80 180 888] 'm ->

   print ["(a ^ b) % 10 ^" m "=" powmod a b 10^m]</lang>

Output:

(a ^ b) % 10 ^ 40 = 1527229998585248450016808958343740453059 
(a ^ b) % 10 ^ 80 = 53259517041910225328867076245698908287781527229998585248450016808958343740453059 
(a ^ b) % 10 ^ 180 = 31857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059 
(a ^ b) % 10 ^ 888 = 261284964380836515397030706363442226571397237057488951313684545241085642329943676248755716124260447188788530017182951051652748425560733974835944416069466176713156182727448301838517000343485327001656948285381173038339073779331230132340669899896448938858785362771190460312412579875349871655999446205426049662261450633448468931573506876255644749155348923523680730999869785472779116009356696816952771965930728940530517799329942590114178284009260298426735086579254282591289756840358811822151307479352856856983393715348870715239020037962938019847992960978849852850613063177471175191444262586321233906926671000476591123695550566585083205841790404069511972417770392822283604206143472509425391114072344402850867571806031857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059

AutoHotkey

Library: MPL

<lang autohotkey>#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</lang>

Output:

1527229998585248450016808958343740453059

BBC BASIC

Works with: BBC BASIC for Windows

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

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

Output:

1527229998585248450016808958343740453059

Bracmat

Translation of: Icon_and_Unicon

<lang bracmat> ( ( mod-power

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

Output:

last 40 digits =  1527229998585248450016808958343740453059

C

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

Library: GMP

<lang c>#include <gmp.h>

int main() { mpz_t a, b, m, r;

mpz_init_set_str(a, "2988348162058574136915891421498819466320" "163312926952423791023078876139", 0); mpz_init_set_str(b, "2351399303373464486466122544523690094744" "975233415544072992656881240319", 0); mpz_init(m); mpz_ui_pow_ui(m, 10, 40);

mpz_init(r); mpz_powm(r, a, b, m);

gmp_printf("%Zd\n", r); /* ...16808958343740453059 */

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

return 0; }</lang> Output:

1527229998585248450016808958343740453059

C#

We can use the built-in function from BigInteger. <lang csharp>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));
   }

}</lang>

Output:

1527229998585248450016808958343740453059

C++

Library: Boost

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

}</lang>

Output:

1527229998585248450016808958343740453059

Clojure

<lang 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))))))</lang>

<lang clojure> (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)))</lang>

Common Lisp

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

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

power (ash power -1)))

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

power (1- power))))

(let ((a 2988348162058574136915891421498819466320163312926952423791023078876139)

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

Works with: CLISP

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

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

Implementation with LOOP

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

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

Crystal

<lang ruby>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)</lang>

Output:

1527229998585248450016808958343740453059

D

Translation of: Icon_and_Unicon

Compile this module with -version=modular_exponentiation to see the output. <lang d>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;
   }</lang>

Output:

1527229998585248450016808958343740453059

Dc

Delphi

Library: System.SysUtils

Library: Velthuis.BigIntegers

Translation of: C#

Thanks for Rudy Velthuis, BigIntegers library [2]. <lang Delphi> 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.</lang>

Output:

1527229998585248450016808958343740453059

EchoLisp

<lang scheme> (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

</lang>

Emacs Lisp

Library: Calc

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

     (b "2351399303373464486466122544523690094744975233415544072992656881240319"))
 ;; "$ ^ $$ mod (10 ^ 40)" performs modular exponentiation.
 ;; "unpack(-5, x)_1" unpacks the integer from the modulo form.
 (message "%s" (calc-eval "unpack(-5, $ ^ $$ mod (10 ^ 40))_1" nil a b)))</lang>

Erlang

<lang 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)).

</lang>

34> modexp_rosetta:test().
modexp_rosetta:test().
1527229998585248450016808958343740453059
35>

F#

<lang fsharp>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</lang>

Factor

<lang factor>! Built-in 2988348162058574136915891421498819466320163312926952423791023078876139 2351399303373464486466122544523690094744975233415544072992656881240319 10 40 ^ ^mod .</lang>

Output:

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

</lang>

Result:
1527229998585248450016808958343740453059
time taken   93.09767815284431 milliseconds
Done

GAP

<lang 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);</lang>

gnuplot

<lang 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</lang>

Go

<lang go>package main

import (

   "fmt"
   "math/big"

)

func main() {

   a, _ := new(big.Int).SetString(
       "2988348162058574136915891421498819466320163312926952423791023078876139", 10)
   b, _ := new(big.Int).SetString(
       "2351399303373464486466122544523690094744975233415544072992656881240319", 10)
   m := big.NewInt(10)
   r := big.NewInt(40)
   m.Exp(m, r, nil)

   r.Exp(a, b, m)
   fmt.Println(r)

}</lang>

Output:

1527229998585248450016808958343740453059

Groovy

<lang groovy>println 2988348162058574136915891421498819466320163312926952423791023078876139.modPow(

           2351399303373464486466122544523690094744975233415544072992656881240319,
           10000000000000000000000000000000000000000)</lang>

Ouput:

1527229998585248450016808958343740453059

Haskell

<lang haskell>powm :: Integer -> Integer -> Integer -> Integer -> Integer powm b 0 m r = r powm b e m r

 | e `mod` 2 == 1 = powm (b * b `mod` m) (e `div` 2) m (r * b `mod` m)

powm b e m r = powm (b * b `mod` m) (e `div` 2) m r

main :: IO () main =

 print $
 powm
   2988348162058574136915891421498819466320163312926952423791023078876139
   2351399303373464486466122544523690094744975233415544072992656881240319
   (10 ^ 40)
   1</lang>

Output:

1527229998585248450016808958343740453059

or in terms of until: <lang haskell>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)</lang>

Output:

1527229998585248450016808958343740453059

Icon and Unicon

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

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

end

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

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

end</lang>

Output:

last 40 digits = 1527229998585248450016808958343740453059

J

Solution:<lang j> m&|@^</lang> Example:<lang j> a =: 2988348162058574136915891421498819466320163312926952423791023078876139x

  b =: 2351399303373464486466122544523690094744975233415544072992656881240319x
  m =: 10^40x

  a m&|@^ b

1527229998585248450016808958343740453059</lang> Discussion: The phrase 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.

<lang java>import java.math.BigInteger;

public class PowMod {

   public static void main(String[] args){
       BigInteger a = new BigInteger(
     "2988348162058574136915891421498819466320163312926952423791023078876139");
       BigInteger b = new BigInteger(
     "2351399303373464486466122544523690094744975233415544072992656881240319");
       BigInteger m = new BigInteger("10000000000000000000000000000000000000000");
       
       System.out.println(a.modPow(b, m));
   }

}</lang>

Output:

1527229998585248450016808958343740453059

Julia

Works with: Julia version 1.0

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

<lang julia>a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = big(10) ^ 40 @show powermod(a, b, m)</lang>

Output:

powermod(a, b, m) = 1527229998585248450016808958343740453059

Kotlin

<lang scala>// 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))

}</lang>

Output:

1527229998585248450016808958343740453059

Lambdatalk

Following scheme

<lang 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 </lang>

Maple

Output:

1527229998585248450016808958343740453059

Mathematica/Wolfram Language

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

Maxima

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

Nim

Library: bigints

<lang nim>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)</lang>

Output:

1527229998585248450016808958343740453059

OCaml

Using the zarith library:

<lang ocaml> 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</lang>

Output:

1527229998585248450016808958343740453059

Oforth

Usage : a b mod powmod

<lang Oforth>: 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 ;</lang>

Output:

>2988348162058574136915891421498819466320163312926952423791023078876139
ok
>2351399303373464486466122544523690094744975233415544072992656881240319
ok
>10 40 pow
ok
>powmod println
1527229998585248450016808958343740453059
ok

ooRexx

<lang rexx>/* Modular exponentiation */

numeric digits 100 say powerMod(,

2988348162058574136915891421498819466320163312926952423791023078876139,,
2351399303373464486466122544523690094744975233415544072992656881240319,,
1e40)

exit

powerMod: procedure

use strict arg base, exponent, modulus

exponent=exponent~d2x~x2b~strip('L','0') result=1 base = base // modulus do exponentPos=exponent~length to 1 by -1

 if (exponent~subChar(exponentPos) == '1')
   then result = (result * base) // modulus
 base = (base * base) // modulus

end return result</lang>

Output:

1527229998585248450016808958343740453059

PARI/GP

Pascal

Works with: Free_Pascal

Library: GMP

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

uses

 gmp;

var

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

begin

 mpz_init_set_str(a, '2988348162058574136915891421498819466320163312926952423791023078876139', 10);
 mpz_init_set_str(b, '2351399303373464486466122544523690094744975233415544072992656881240319', 10);
 mpz_init(m);
 mpz_ui_pow_ui(m, 10, 40);

 mpz_init(r);
 mpz_powm(r, a, b, m);

 fmt := '%Zd' + chr(13) + chr(10);
 mp_printf(fmt, @r); (* ...16808958343740453059 *)
 
 mpz_clear(a);
 mpz_clear(b);
 mpz_clear(m);
 mpz_clear(r);

end.</lang>

Output:

% ./ModularExponentiation
1527229998585248450016808958343740453059

Perl

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. <lang perl>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";</lang>

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

Output:

1527229998585248450016808958343740453059

PicoLisp

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

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

Test: <lang PicoLisp>: (**Mod

  2988348162058574136915891421498819466320163312926952423791023078876139
  2351399303373464486466122544523690094744975233415544072992656881240319
  10000000000000000000000000000000000000000 )

-> 1527229998585248450016808958343740453059</lang>

Prolog

Works with: SWI Prolog

SWI Prolog has a built-in function named powm for this purpose. <lang prolog>main:-

   A = 2988348162058574136915891421498819466320163312926952423791023078876139,
   B = 2351399303373464486466122544523690094744975233415544072992656881240319,
   M is 10 ** 40,
   P is powm(A, B, M),
   writeln(P).</lang>

Output:

1527229998585248450016808958343740453059

Python

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

Output:

1527229998585248450016808958343740453059

<lang python>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))</lang>

Output:

1527229998585248450016808958343740453059

Quackery

<lang 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</lang>

Output:

1527229998585248450016808958343740453059

Racket

lang racket

(require math) (define a 2988348162058574136915891421498819466320163312926952423791023078876139) (define b 2351399303373464486466122544523690094744975233415544072992656881240319) (define m (expt 10 40)) (modular-expt a b m) </lang>

Output:

1527229998585248450016808958343740453059

Raku

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

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

}

say expmod

   2988348162058574136915891421498819466320163312926952423791023078876139,
   2351399303373464486466122544523690094744975233415544072992656881240319,
   10**40;</lang>

Output:

1527229998585248450016808958343740453059

REXX

version 1

This REXX program attempts to handle any a, b, or m, but there are limits for any computer language.

For some REXXes, it's around eight million digits for any arithmetic expression or value, which puts limitations on
the values of a² or 10^m.

This REXX program code has code to automatically adjust the number of decimal digits to accommodate huge
numbers which are computed when raising large numbers to some arbitrary power. <lang rexx>/*REXX program displays the modular exponentiation of: a**b mod m */ parse arg a b m /*obtain optional args from the CL*/ if a== | a=="," then a= 2988348162058574136915891421498819466320163312926952423791023078876139 if b== | b=="," then b= 2351399303373464486466122544523690094744975233415544072992656881240319 if m = | m ="," then m= 40 /*Not specified? Then use default.*/ say 'a=' a " ("length(a) 'digits)' /*display the value of A (&length)*/ say 'b=' b " ("length(b) 'digits)' /* " " " " B " */

      do j=1  for words(m);           y= word(m, j)  /*use one of the MM powers (list).*/
      say copies('═', linesize() - 1)                /*show a nice separator fence line*/
      say 'a**b (mod 10**'y")="  powerMod(a,b,10**y) /*display the answer ───► console.*/
      end   /*j*/

exit 0 /*stick a fork in it, we're done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ powerMod: procedure; parse arg x,p,mm /*fast modular exponentiation code*/

         parse value max(x*x, p, mm)'E0'  with "E" e /*obtain the biggest of the three.*/
         numeric digits max(40, e+e)                 /*ensure big enough to handle  A².*/
         $= 1                                        /*use this for the first value.   */
                 do  until p==0                      /*perform until P is equal to zero*/
                 if p // 2  then $= $ * x  //  mm    /*is P odd?  (is ÷ remainder ≡ 1?)*/
                 p= p % 2;       x= x * x  //  mm    /*halve  P;   calculate  x² mod n */
                 end   /*until*/;         return $   /* [↑]  keep mod'ing until P=zero.*/</lang>

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

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

a= 2988348162058574136915891421498819466320163312926952423791023078876139      (70 digits)
b= 2351399303373464486466122544523690094744975233415544072992656881240319      (70 digits)
═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
a**b (mod 10**40)= 1527229998585248450016808958343740453059
═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
a**b (mod 10**80)= 53259517041910225328867076245698908287781527229998585248450016808958343740453059
═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
a**b (mod 10**180)= 31857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059
═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
a**b (mod 10**888)= 261284964380836515397030706363442226571397237057488951313684545241085642329943676248755716124260447188788530017182951051652748425560733974835944416069466176713156182727448301838517000343485327001656948285381173038339073779331230132340669899896448938858785362771190460312412579875349871655999446205426049662261450633448468931573506876255644749155348923523680730999869785472779116009356696816952771965930728940530517799329942590114178284009260298426735086579254282591289756840358811
822151307479352856856983393715348870715239020037962938019847992960978849852850613063177471175191444262586321233906926671000476591123695550566585083205841790404069511972417770392822283604206143472509425391114072344402850867571806031857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059

version 2

This REXX version handles only up to 100 decimal digits. <lang rexx> /* REXX Modular exponentiation */

numeric digits 100 say powerMod(,

2988348162058574136915891421498819466320163312926952423791023078876139,, 
2351399303373464486466122544523690094744975233415544072992656881240319,, 
1e40)

exit

powerMod: procedure

parse arg base, exponent, modulus

exponent = strip(x2b(d2x(exponent)), 'L', '0') result = 1 base = base // modulus do exponentPos=length(exponent) to 1 by -1

 if substr(exponent, exponentPos, 1) = '1'    
   then result = (result * base) // modulus  
 base = (base * base) // modulus

end return result</lang>

Output:

1527229998585248450016808958343740453059

Ruby

Built in since version 2.5.

Using OpenSSL standard library

Library: OpenSSL

<lang ruby>require 'openssl' a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 puts a.to_bn.mod_exp(b, m)</lang>

Written in Ruby

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

</lang>

Rust

<lang 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;
   }

}</lang>

Test code:

<lang rust>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);

}</lang>

Output:

The last 40 digits of
2988348162058574136915891421498819466320163312926952423791023078876139
to the power of
2351399303373464486466122544523690094744975233415544072992656881240319
are:
1527229998585248450016808958343740453059

Scala

<lang scala>import scala.math.BigInt

val a = BigInt(

 "2988348162058574136915891421498819466320163312926952423791023078876139")

val b = BigInt(

 "2351399303373464486466122544523690094744975233415544072992656881240319")

println(a.modPow(b, BigInt(10).pow(40)))</lang>

Scheme

<lang 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)))</lang>

Output:

> result
1527229998585248450016808958343740453059

Seed7

The library bigint.s7i defines the function modPow, which does modular exponentiation. <lang seed7>$ include "seed7_05.s7i";

 include "bigint.s7i";

const proc: main is func

 begin
   writeln(modPow(2988348162058574136915891421498819466320163312926952423791023078876139_,
                  2351399303373464486466122544523690094744975233415544072992656881240319_,
                  10_ ** 40));
 end func;</lang>

Output:

1527229998585248450016808958343740453059

The library bigint.s7i defines modPow with: <lang seed7>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;</lang>

Original source: [3]

Sidef

Built-in: <lang ruby>say expmod(

   2988348162058574136915891421498819466320163312926952423791023078876139,
   2351399303373464486466122544523690094744975233415544072992656881240319,
   10**40)</lang>

User-defined: <lang ruby>func expmod(a, b, n) {

   var c = 1
   do {
       (c *= a) %= n if b.is_odd
       (a *= a) %= n
   } while (b //= 2)
   c

}</lang>

Output:

1527229998585248450016808958343740453059

Swift

Library: AttaSwift BigInt

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

<lang swift>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)))</lang>

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

<lang tcl>package require Tcl 8.5

Algorithm from http://introcs.cs.princeton.edu/java/78crypto/ModExp.java.html
but Tcl has arbitrary-width integers and an exponentiation operator, which
helps simplify the code.

proc tcl::mathfunc::modexp {a b n} {

   if {$b == 0} {return 1}
   set c [expr {modexp($a, $b / 2, $n)**2 % $n}]
   if {$b & 1} {

set c [expr {($c * $a) % $n}]

   }
   return $c

}</lang> Demonstrating: <lang tcl>set a 2988348162058574136915891421498819466320163312926952423791023078876139 set b 2351399303373464486466122544523690094744975233415544072992656881240319 set n [expr {10**40}] puts [expr {modexp($a,$b,$n)}]</lang>

Output:

 1527229998585248450016808958343740453059

Iterative

<lang tcl>package require Tcl 8.5 proc modexp {a b n} {

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

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

   }
   return $c

}</lang> Demonstrating: <lang tcl>set a 2988348162058574136915891421498819466320163312926952423791023078876139 set b 2351399303373464486466122544523690094744975233415544072992656881240319 set n [expr {10**40}] puts [modexp $a $b $n]</lang>

Output:

 1527229998585248450016808958343740453059

TXR

From your system prompt:

<lang sh>$ txr -p '(exptmod 2988348162058574136915891421498819466320163312926952423791023078876139

                  2351399303373464486466122544523690094744975233415544072992656881240319
                  (expt 10 40)))'

1527229998585248450016808958343740453059</lang>

Visual Basic .NET

Works with: Visual Basic .NET version 2011

<lang vbnet>' 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</lang>

Output:

x=1527229998585248450016808958343740453059

Wren

Library: Wren-big

<lang ecmascript>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))</lang>

Output:

1527229998585248450016808958343740453059

zkl

Using the GMP big num library: <lang zkl>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();</lang>

Output:

1527229998585248450016808958343740453059
1,527,229,998,585,248,450,016,808,958,343,740,453,059

@@ Line 1,032: / Line 1,032: @@
 <pre>1527229998585248450016808958343740453059</pre>
-=={{header|Mathematica}}==
+=={{header|Mathematica}}/{{header|Wolfram Language}}==
 <lang Mathematica>a = 2988348162058574136915891421498819466320163312926952423791023078876139;
 b = 2351399303373464486466122544523690094744975233415544072992656881240319;