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 .
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
AutoHotkey
<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
Bracmat
<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
BBC BASIC
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
C
Given numbers are too big for even 64 bit integers, so might as well take the lazy route and use 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>
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>
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>
<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>
D
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
<lang Dc>2988348162058574136915891421498819466320163312926952423791023078876139 2351399303373464486466122544523690094744975233415544072992656881240319 10 40^|p</lang>
Emacs Lisp
<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>
F#
<lang fsharp>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</lang>
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 = putStrLn . show $
powm 2988348162058574136915891421498819466320163312926952423791023078876139 2351399303373464486466122544523690094744975233415544072992656881240319 (10 ^ 40) 1</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
We can use the built-in powermod
function with the built-in BigInt
type (based on GMP):
<lang julia>julia> a = BigInt("2988348162058574136915891421498819466320163312926952423791023078876139")
b = BigInt("2351399303373464486466122544523690094744975233415544072992656881240319") m = BigInt(10)^40 powermod(a, b, m)
1527229998585248450016808958343740453059</lang>
Maple
<lang Maple>a := 2988348162058574136915891421498819466320163312926952423791023078876139: b := 2351399303373464486466122544523690094744975233415544072992656881240319: a &^ b mod 10^40;</lang>
- Output:
1527229998585248450016808958343740453059
Mathematica
<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
<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
Oforth
Usage : a b mod powmod
<lang Oforth>: powmod(modulus, exponent, base) {
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
<lang parigp>a=2988348162058574136915891421498819466320163312926952423791023078876139; b=2351399303373464486466122544523690094744975233415544072992656881240319; lift(Mod(a,10^40)^b)</lang>
Pascal
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
<lang perl>use bigint;
my $a = 2988348162058574136915891421498819466320163312926952423791023078876139; my $b = 2351399303373464486466122544523690094744975233415544072992656881240319; my $m = 10 ** 40; print $a->bmodpow($b, $m), "\n";</lang>
- Output:
1527229998585248450016808958343740453059
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
PHP
<lang php><?php $a = '2988348162058574136915891421498819466320163312926952423791023078876139'; $b = '2351399303373464486466122544523690094744975233415544072992656881240319'; $m = '1' . str_repeat('0', 40); echo bcpowmod($a, $b, $m), "\n";</lang>
- 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>
Python
<lang python>a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 print(pow(a, b, m))</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
Racket
<lang racket>
- lang racket
(require math) (define a 2988348162058574136915891421498819466320163312926952423791023078876139) (define b 2351399303373464486466122544523690094744975233415544072992656881240319) (define m (expt 10 40)) (modular-expt a b m) </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 REXX, it's around eight million digits for any arithmetic expression or value, which puts limitations on the values of or .
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. <lang rexx>/*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 _</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 LINESIZE.REX REXX program is included here ──► LINESIZE.REX.
- Output:
when using the input of
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
version 2
<lang 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
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.
<lang ruby>require 'openssl'
a = 2988348162058574136915891421498819466320163312926952423791023078876139 b = 2351399303373464486466122544523690094744975233415544072992656881240319 m = 10 ** 40 puts a.to_bn.mod_exp(b, m)</lang>
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?.)
<lang ruby>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)</lang>
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: [2]
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>
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
- Programming Tasks
- Solutions by Programming Task
- Ada
- AutoHotkey
- MPL
- Bracmat
- BBC BASIC
- C
- GMP
- Clojure
- Common Lisp
- D
- Dc
- Emacs Lisp
- Calc
- F Sharp
- GAP
- Gnuplot
- Go
- Groovy
- Haskell
- Icon
- Unicon
- J
- Java
- Julia
- Maple
- Mathematica
- Maxima
- Nim
- Bigints
- Oforth
- OoRexx
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- Python
- OCaml
- Racket
- REXX
- Ruby
- OpenSSL
- Scala
- Scheme
- Seed7
- Tcl
- TXR
- Zkl