Modular exponentiation

From Rosetta Code
Jump to: navigation, search
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 ab, where

  • a = 2988348162058574136915891421498819466320163312926952423791023078876139
  • b = 2351399303373464486466122544523690094744975233415544072992656881240319

A computer is too slow to find the entire value of ab. 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 \ge 0 and m > 0.

Contents

[edit] 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

[edit] 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

[edit] 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

[edit] 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

[edit] 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;
}

[edit] 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))))))

[edit] 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))))

[edit] 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)))

[edit] 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

[edit] Dc

2988348162058574136915891421498819466320163312926952423791023078876139 2351399303373464486466122544523690094744975233415544072992656881240319 10 40^|p

[edit] 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)))

[edit] F#

let expMod a b n =
let mulMod x y n = snd (bigint.DivRem(x * y, n))
let rec loop a b c =
if b = 0I then c else
let (bd, br) = bigint.DivRem(b, 2I)
loop (mulMod a a n) bd (if br = 0I then c else (mulMod 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

[edit] 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);

[edit] 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

[edit] Groovy

println 2988348162058574136915891421498819466320163312926952423791023078876139.modPow(
2351399303373464486466122544523690094744975233415544072992656881240319,
10000000000000000000000000000000000000000)

Ouput:

1527229998585248450016808958343740453059

[edit] 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 = putStrLn . show $
powm
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
(10 ^ 40)
1
Output:
1527229998585248450016808958343740453059

[edit] 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

[edit] J

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

Discussion: The phrase m&|@^ is the natural expression of a^b mod m in J, and is recognized by the interpreter as an opportunity for optimization, by avoiding the exponentiation.

[edit] 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

[edit] Julia

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

julia> a = BigInt("2988348162058574136915891421498819466320163312926952423791023078876139")
b = BigInt("2351399303373464486466122544523690094744975233415544072992656881240319")
m = BigInt(10)^40
powermod(a, b, m)
1527229998585248450016808958343740453059

[edit] Maple

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

[edit] Mathematica

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

[edit] Maxima

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

[edit] Nimrod

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

[edit] ooRexx

/* 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
Output:
1527229998585248450016808958343740453059

[edit] PARI/GP

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

[edit] 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

[edit] Perl

use bigint;
 
my $a = 2988348162058574136915891421498819466320163312926952423791023078876139;
my $b = 2351399303373464486466122544523690094744975233415544072992656881240319;
my $m = 10 ** 40;
print $a->bmodpow($b, $m), "\n";
Output:
1527229998585248450016808958343740453059

[edit] 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

[edit] PHP

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

[edit] 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

[edit] Python

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

[edit] Racket

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

[edit] REXX

[edit] version 1

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

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

/*REXX program displays  modular exponentation:     a**b  mod  M        */
parse arg a b mm /*get optional args from CL.*/
if a=='' | a==',' then a=2988348162058574136915891421498819466320163312926952423791023078876139
if b=='' | b==',' then b=2351399303373464486466122544523690094744975233415544072992656881240319
if mm='' | mm=',' then mm=40 /*MM specified? Use default.*/
say 'a=' a; say ' ('length(a) "digits)" /*show value of A.*/
say 'b=' b; say ' ('length(b) "digits)" /* " " " B.*/
 
do j=1 for words(mm); m=word(mm,j) /*use one of the MM powers.*/
say copies('─',linesize()-1) /*show a nice separator line*/
say 'a**b (mod 10**'m")=" powerModulated(a,b,10**m) /*show ans*/
end /*j*/
exit /*stick a fork in it; done.*/
/*──────────────────────────────────────POWERMODULATED subroutine───────*/
powerModulated: procedure; parse arg x,p,n /*fast modular exponentation*/
if p==0 then return 1 /*special case of P = zero. */
if p==1 then return x /* " " " " = unity.*/
if p<0 then do; say '***error!*** power is negative:' p; exit 13; end
parse value max(x**2,p,n)'E0' with "E" e /*pick biggest of the three.*/
numeric digits max(20,e*2) /*big enough to handle A² */
_=1 /*use this for the 1st value*/
do while p\==0; if p//2==1 then _=_*x//n /*is P odd? */
p=p%2; x=x*x//n /*halve P; calc x² mod n */
end /*while*/ /* [↑] keep moding 'til =0.*/
return _

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

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

Note the REXX program was executing within a window of 600 bytes wide.

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

[edit] version 2

/* 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
Output:
1527229998585248450016808958343740453059

[edit] Ruby

Ruby's core library has no modular exponentiation. OpenSSL, in Ruby's standard library, provides OpenSSL::BN#mod_exp. To reach this method, we call Integer#to_bn to convert a from Integer to OpenSSL::BN. The method implicitly converts b and m.

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

Or we can implement a custom method, Integer#rosetta_mod_exp, to calculate the same result. This method does exponentiation by successive squaring, but replaces each intermediate product with a congruent value. (Program needs Ruby 1.8.7 for Integer#odd?.)

Works with: Ruby version 1.8.7
class Integer
def rosetta_mod_exp(exp, mod)
exp < 0 and raise ArgumentError, "negative exponent"
prod = 1
base = self % mod
until exp.zero?
exp.odd? and prod = (prod * base) % mod
exp >>= 1
base = (base * base) % mod
end
prod
end
end
 
a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ** 40
puts a.rosetta_mod_exp(b, m)

[edit] Scala

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

[edit] 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

[edit] Seed7

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

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

The library bigint.s7i defines modPow with:

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

Original source: [2]

[edit] 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.

[edit] 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

[edit] 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

[edit] TXR

From your system prompt:

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

[edit] 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
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox