Modular exponentiation: Difference between revisions
m
→{{header|Wren}}: Minor tidy
(Include Erlang solution) |
m (→{{header|Wren}}: Minor tidy) |
||
| (23 intermediate revisions by 11 users not shown) | |||
Line 1:
{{task}}
Find the last '''40''' decimal digits of <math>a^b</math>, where
::* <math>a = 2988348162058574136915891421498819466320163312926952423791023078876139</math>
::* <math>b = 2351399303373464486466122544523690094744975233415544072992656881240319</math>
The algorithm must work for any integers <math>a, b, m</math>, where <math>b \ge 0</math> and <math>m > 0</math>.
<br><br>
Line 16:
{{trans|D}}
<
BigInt result = 1
Line 29:
print(pow_mod(BigInt(‘2988348162058574136915891421498819466320163312926952423791023078876139’),
BigInt(‘2351399303373464486466122544523690094744975233415544072992656881240319’),
BigInt(10) ^ 40))</
{{out}}
Line 40:
Using the big integer implementation from a cryptographic library [https://github.com/cforler/Ada-Crypto-Library/].
<
procedure Mod_Exp is
Line 65:
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;</
{{out}}
Line 73:
The code below uses Algol 68 Genie which provides arbitrary precision arithmetic for LONG LONG modes.
<
BEGIN
PR precision=1000 PR
Line 97:
printf (($"Last 40 digits = ", 40dl$, mod power (a, b, m)))
END
</syntaxhighlight>
{{out}}
Line 105:
=={{header|Arturo}}==
<
b: 2351399303373464486466122544523690094744975233415544072992656881240319
loop [40 80 180 888] 'm ->
print ["(a ^ b) % 10 ^" m "=" powmod a b 10^m]</
{{out}}
Line 117:
(a ^ b) % 10 ^ 180 = 31857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059
(a ^ b) % 10 ^ 888 = 261284964380836515397030706363442226571397237057488951313684545241085642329943676248755716124260447188788530017182951051652748425560733974835944416069466176713156182727448301838517000343485327001656948285381173038339073779331230132340669899896448938858785362771190460312412579875349871655999446205426049662261450633448468931573506876255644749155348923523680730999869785472779116009356696816952771965930728940530517799329942590114178284009260298426735086579254282591289756840358811822151307479352856856983393715348870715239020037962938019847992960978849852850613063177471175191444262586321233906926671000476591123695550566585083205841790404069511972417770392822283604206143472509425391114072344402850867571806031857295076204937005344367438778481743660325586328069392203762862423884839076695547212682454523811053259517041910225328867076245698908287781527229998585248450016808958343740453059</pre>
=={{header|ATS}}==
{{libheader|ats2-xprelude}}
{{libheader|GMP}}
For its multiple precision support, you will need [https://sourceforge.net/p/chemoelectric/ats2-xprelude/ ats2-xprelude] with GNU MP support enabled.
There is GNU MP support that comes with the ATS2 compiler, however, so perhaps someone will write a demonstration using that. Unlike ats2-xprelude (which assumes there is garbage collection), that support represents multi-precision numbers as linear types.
<syntaxhighlight lang="ats">
(* You will need
https://sourceforge.net/p/chemoelectric/ats2-xprelude/ *)
#include "share/atspre_staload.hats"
#include "xprelude/HATS/xprelude.hats"
staload "xprelude/SATS/exrat.sats"
staload _ = "xprelude/DATS/exrat.dats"
val a = exrat_make_string_exn "2988348162058574136915891421498819466320163312926952423791023078876139"
val b = exrat_make_string_exn "2351399303373464486466122544523690094744975233415544072992656881240319"
val modulus = exrat_make (10, 1) ** 40
(* xprelude/SATS/exrat.sats includes the "exrat_numerator_modular_pow"
function, based on GMP's mpz_powm. *)
val result1 = exrat_numerator_modular_pow (a, b, modulus)
(* But that was too easy. Here is the right-to-left binary method,
https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
*)
val result2 =
(lam (base : exrat,
exponent : exrat,
modulus : exrat) : exrat =>
let
val zero = exrat_make (0, 1)
and one = exrat_make (1, 1)
and two = exrat_make (2, 1)
macdef divrem = exrat_numerator_euclid_division
macdef rem = exrat_numerator_euclid_remainder
in
if modulus = one then
zero
else
let
fun
loop (result : exrat,
base : exrat,
exponent : exrat) : exrat =
if iseqz exponent then
result
else
let
val @(exponent, remainder) = exponent \divrem two
val result =
if remainder = one then
(result * base) \rem modulus
else
result
val base = (base * base) \rem modulus
in
loop (result, base, exponent)
end
in
loop (one, base \rem modulus, exponent)
end
end) (a, b, modulus)
implement
main0 () =
begin
println! result1;
println! result2
end
</syntaxhighlight>
{{out}}
<pre>$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW $(pkg-config --cflags ats2-xprelude) $(pkg-config --variable=PATSCCFLAGS ats2-xprelude) modular-exponentiation.dats $(pkg-config --libs ats2-xprelude) -lgc && ./a.out
1527229998585248450016808958343740453059
1527229998585248450016808958343740453059
</pre>
=={{header|AutoHotkey}}==
{{libheader|MPL}}
<
#SingleInstance, Force
SetBatchLines, -1
Line 145 ⟶ 227:
msgbox % MP_DEC(result)
Return</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
Line 152 ⟶ 234:
{{works with|BBC BASIC for Windows}}
Uses the Huge Integer Math & Encryption library.
<
PROC_himeinit("")
Line 160 ⟶ 242:
h1% = 1 : h2% = 2 : h3% = 3 : h4% = 4
SYS `hi_PowMod`, ^h1%, ^h2%, ^h3%, ^h4%
PRINT FN_higetdec(4)</
{{out}}
<pre>
1527229998585248450016808958343740453059
</pre>
=={{header|bc}}==
<syntaxhighlight lang="bc">define p(n, e, m) {
auto r
for (r = 1; e > 0; e /= 2) {
if (e % 2 == 1) r = n * r % m
n = n * n % m
}
return(r)
}
a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
p(a, b, 10 ^ 40)</syntaxhighlight>
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|Bracmat}}==
{{trans|Icon_and_Unicon}}
<
= base exponent modulus result
. !arg:(?base,?exponent,?modulus)
Line 195 ⟶ 293:
)
& out$("last 40 digits = " mod-power$(!a,!b,10^40))
)</
{{out}}
<pre>last 40 digits = 1527229998585248450016808958343740453059</pre>
Line 202 ⟶ 300:
Given numbers are too big for even 64 bit integers, so might as well take the lazy route and use GMP:
{{libheader|GMP}}
<
int main()
Line 226 ⟶ 324:
return 0;
}</
Output:
<pre>
Line 234 ⟶ 332:
=={{header|C sharp}}==
We can use the built-in function from BigInteger.
<
using System.Numerics;
Line 245 ⟶ 343:
Console.WriteLine(BigInteger.ModPow(a, b, m));
}
}</
{{out}}
<pre>
Line 253 ⟶ 351:
=={{header|C++}}==
{{libheader|Boost}}
<
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/integer.hpp>
Line 265 ⟶ 363:
std::cout << powm(a, b, pow(cpp_int(10), 40)) << '\n';
return 0;
}</
{{out}}
Line 273 ⟶ 371:
=={{header|Clojure}}==
<
(defn m* [p q] (mod (* p q) m))
(loop [b b, e e, x 1]
(if (zero? e) x
(if (even? e) (recur (m* b b) (/ e 2) x)
(recur (m* b b) (quot e 2) (m* b x))))))</
<
(defn modpow
" b^e mod m (using Java which solves some cases the pure clojure method has to be modified to tackle--i.e. with large b & e and
calculation simplications when gcd(b, m) == 1 and gcd(e, m) == 1) "
[b e m]
(.modPow (biginteger b) (biginteger e) (biginteger m)))</
=={{header|Common Lisp}}==
<
"Return BASE raised to the POWER, modulo DIVISOR.
This function is faster than (MOD (EXPT BASE POWER) DIVISOR), but
Line 305 ⟶ 403:
(let ((a 2988348162058574136915891421498819466320163312926952423791023078876139)
(b 2351399303373464486466122544523690094744975233415544072992656881240319))
(format t "~A~%" (rosetta-mod-expt a b (expt 10 40))))</
{{works with|CLISP}}
<
(let ((a 2988348162058574136915891421498819466320163312926952423791023078876139)
(b 2351399303373464486466122544523690094744975233415544072992656881240319))
(format t "~A~%" (mod-expt a b (expt 10 40))))</
===Implementation with LOOP===
<
(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)))</
=={{header|Crystal}}==
<
module Integer
Line 346 ⟶ 444:
m = 10.to_big_i ** 40
puts a.powmod(b, m)</
{{out}}
Line 354 ⟶ 452:
{{trans|Icon_and_Unicon}}
Compile this module with <code>-version=modular_exponentiation</code> to see the output.
<
private import std.bigint;
Line 383 ⟶ 481:
"4744975233415544072992656881240319"),
BigInt(10) ^^ 40).writeln;
}</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|Dc}}==
<syntaxhighlight lang
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
Line 394 ⟶ 492:
{{Trans|C#}}
Thanks for Rudy Velthuis, BigIntegers library [https://github.com/rvelthuis/DelphiBigNumbers].
<syntaxhighlight lang="delphi">
program Modular_exponentiation;
Line 412 ⟶ 510:
Writeln(BigInteger.ModPow(a, b, m).ToString);
readln;
end.</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|EchoLisp}}==
<
(lib 'bigint)
Line 440 ⟶ 538:
(xpowmod a b m)
→ 1527229998585248450016808958343740453059
</syntaxhighlight>
=={{header|Emacs Lisp}}==
{{libheader|Calc}}
<
(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)))</
=={{header|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].
Line 469 ⟶ 570:
Result;
binary_exp(Base,Exponent,Acc) ->
binary_exp(Base*Base,Exponent bsr 1,Acc
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
Exponent bsr 1,Mod,(Acc * exp_factor(Base,Exponent))rem Mod).
exp_factor(_,0) ->
1;
Base;
exp_factor(Base,Exponent) ->
exp_factor(Base,Exponent band 1).
test() ->
Line 496 ⟶ 595:
binary_exp(10,40)).
</syntaxhighlight>
<pre>
34> modexp_rosetta:test().
Line 503 ⟶ 602:
35>
</pre>
=={{header|F_Sharp|F#}}==
<
let rec loop a b c =
if b = 0I then c else
Line 515 ⟶ 615:
let b = 2351399303373464486466122544523690094744975233415544072992656881240319I
printfn "%A" (expMod a b (10I**40)) // -> 1527229998585248450016808958343740453059
0</
=={{header|Factor}}==
<
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
10 40 ^
^mod .</
{{out}}
<pre>
1527229998585248450016808958343740453059
</pre>
=={{header|Fortran}}==
<syntaxhighlight lang="Fortran">
module big_integers ! Big (but not very big) integers.
!
! A very primitive multiple precision module, using the classical
! algorithms (long multiplication, long division) and a mere 8-bit
! "digit" size. Some procedures might assume integers are in a
! two's-complement representation. This module is good enough for us
! to fulfill the task.
!
! NOTE: I assume that iachar and achar do not alter the most
! significant bit.
use, intrinsic :: iso_fortran_env, only: int16
implicit none
private
public :: big_integer
public :: integer2big
public :: string2big
public :: big2string
public :: big_sgn
public :: big_cmp, big_cmpabs
public :: big_neg, big_abs
public :: big_addabs, big_add
public :: big_subabs, big_sub
public :: big_mul ! One might also include a big_muladd.
public :: big_divrem ! One could also include big_div and big_rem.
public :: big_pow
public :: operator(+)
public :: operator(-)
public :: operator(*)
public :: operator(**)
type :: big_integer
! The representation is sign-magnitude. The radix is 256, which
! is not speed-efficient, but which seemed relatively easy to
! work with if one were writing in standard Fortran (and assuming
! iachar and achar were "8-bit clean").
logical :: sign = .false. ! .false. => +sign, .true. => -sign.
character, allocatable :: bytes(:)
end type big_integer
character, parameter :: zero = achar (0)
character, parameter :: one = achar (1)
! An integer type capable of holding an unsigned 8-bit value.
integer, parameter :: bytekind = int16
interface operator(+)
module procedure big_add
end interface
interface operator(-)
module procedure big_neg
module procedure big_sub
end interface
interface operator(*)
module procedure big_mul
end interface
interface operator(**)
module procedure big_pow
end interface
contains
elemental function logical2byte (bool) result (byte)
logical, intent(in) :: bool
character :: byte
if (bool) then
byte = one
else
byte = zero
end if
end function logical2byte
elemental function logical2i (bool) result (i)
logical, intent(in) :: bool
integer :: i
if (bool) then
i = 1
else
i = 0
end if
end function logical2i
elemental function byte2i (c) result (i)
character, intent(in) :: c
integer :: i
i = iachar (c)
end function byte2i
elemental function i2byte (i) result (c)
integer, intent(in) :: i
character :: c
c = achar (i)
end function i2byte
elemental function byte2bk (c) result (i)
character, intent(in) :: c
integer(bytekind) :: i
i = iachar (c, kind = bytekind)
end function byte2bk
elemental function bk2byte (i) result (c)
integer(bytekind), intent(in) :: i
character :: c
c = achar (i)
end function bk2byte
elemental function bk2i (i) result (j)
integer(bytekind), intent(in) :: i
integer :: j
j = int (i)
end function bk2i
elemental function i2bk (i) result (j)
integer, intent(in) :: i
integer(bytekind) :: j
j = int (iand (i, 255), kind = bytekind)
end function i2bk
! Left shift of the least significant 8 bits of a bytekind integer.
elemental function lshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ishft (ibits (a, 0, 8 - i), i)
end function lshftbk
! Right shift of the least significant 8 bits of a bytekind integer.
elemental function rshftbk (a, i) result (c)
integer(bytekind), intent(in) :: a
integer, intent(in) :: i
integer(bytekind) :: c
c = ibits (a, i, 8 - i)
end function rshftbk
! Left shift an integer.
elemental function lshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, i)
end function lshfti
! Right shift an integer.
elemental function rshfti (a, i) result (c)
integer, intent(in) :: a
integer, intent(in) :: i
integer :: c
c = ishft (a, -i)
end function rshfti
function integer2big (i) result (a)
integer, intent(in) :: i
type(big_integer), allocatable :: a
!
! To write a more efficient implementation of this procedure is
! left as an exercise for the reader.
!
character(len = 100) :: buffer
write (buffer, '(I0)') i
a = string2big (trim (buffer))
end function integer2big
function string2big (s) result (a)
character(len = *), intent(in) :: s
type(big_integer), allocatable :: a
integer :: n, i, istart, iend
integer :: digit
if ((s(1:1) == '-') .or. s(1:1) == '+') then
istart = 2
else
istart = 1
end if
iend = len (s)
n = (iend - istart + 2) / 2
allocate (a)
allocate (a%bytes(n))
a%bytes = zero
do i = istart, iend
digit = ichar (s(i:i)) - ichar ('0')
if (digit < 0 .or. 9 < digit) error stop
a = short_multiplication (a, 10)
a = short_addition (a, digit)
end do
a%sign = (s(1:1) == '-')
call normalize (a)
end function string2big
function big2string (a) result (s)
type(big_integer), intent(in) :: a
character(len = :), allocatable :: s
type(big_integer), allocatable :: q
integer :: r
integer :: sgn
sgn = big_sgn (a)
if (sgn == 0) then
s = '0'
else
q = a
s = ''
do while (big_sgn (q) /= 0)
call short_division (q, 10, q, r)
s = achar (r + ichar ('0')) // s
end do
if (sgn < 0) s = '-' // s
end if
end function big2string
function big_sgn (a) result (sgn)
type(big_integer), intent(in) :: a
integer :: sgn
integer :: n, i
n = size (a%bytes)
i = 1
sgn = 1234
do while (sgn == 1234)
if (i == n + 1) then
sgn = 0
else if (a%bytes(i) /= zero) then
if (a%sign) then
sgn = -1
else
sgn = 1
end if
else
i = i + 1
end if
end do
end function big_sgn
function big_cmp (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp
if (a%sign) then
if (b%sign) then
cmp = -big_cmpabs (a, b)
else
cmp = -1
end if
else
if (b%sign) then
cmp = 1
else
cmp = big_cmpabs (a, b)
end if
end if
end function big_cmp
function big_cmpabs (a, b) result (cmp)
type(big_integer(*)), intent(in) :: a, b
integer :: cmp
integer :: n, i
integer :: ia, ib
cmp = 1234
n = max (size (a%bytes), size (b%bytes))
i = n
do while (cmp == 1234)
if (i == 0) then
cmp = 0
else
ia = byteval (a, i)
ib = byteval (b, i)
if (ia < ib) then
cmp = -1
else if (ia > ib) then
cmp = 1
else
i = i - 1
end if
end if
end do
end function big_cmpabs
function big_neg (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .not. c%sign
end function big_neg
function big_abs (a) result (c)
type(big_integer), intent(in) :: a
type(big_integer), allocatable :: c
c = a
c%sign = .false.
end function big_abs
function big_add (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c
logical :: sign
if (a%sign) then
if (b%sign) then ! a <= 0, b <= 0
c = big_addabs (a, b)
sign = .true.
else ! a <= 0, b >= 0
c = big_subabs (a, b)
sign = .not. c%sign
end if
else
if (b%sign) then ! a >= 0, b <= 0
c = big_subabs (a, b)
sign = c%sign
else ! a >= 0, b >= 0
c = big_addabs (a, b)
sign = .false.
end if
end if
c%sign = sign
end function big_add
function big_sub (a, b) result (c)
type(big_integer), intent(in) :: a
type(big_integer), intent(in) :: b
type(big_integer), allocatable :: c
logical :: sign
if (a%sign) then
if (b%sign) then ! a <= 0, b <= 0
c = big_subabs (a, b)
sign = .not. c%sign
else ! a <= 0, b >= 0
c = big_addabs (a, b)
sign = .true.
end if
else
if (b%sign) then ! a >= 0, b <= 0
c = big_addabs (a, b)
sign = .false.
else ! a >= 0, b >= 0
c = big_subabs (a, b)
sign = c%sign
end if
end if
c%sign = sign
end function big_sub
function big_addabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c
! Compute abs(a) + abs(b).
integer :: n, nc, i
logical :: carry
type(big_integer), allocatable :: tmp
n = max (size (a%bytes), size (b%bytes))
nc = n + 1
allocate(tmp)
allocate(tmp%bytes(nc))
call add_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call add_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp
end function big_addabs
function big_subabs (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c
! Compute abs(a) - abs(b). The result is signed.
integer :: n, i
logical :: carry
type(big_integer), allocatable :: tmp
n = max (size (a%bytes), size (b%bytes))
allocate(tmp)
allocate(tmp%bytes(n))
if (big_cmpabs (a, b) >= 0) then
tmp%sign = .false.
call sub_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
end do
else
tmp%sign = .true.
call sub_bytes (get_byte (b, 1), get_byte (a, 1), .false., tmp%bytes(1), carry)
do i = 2, n
call sub_bytes (get_byte (b, i), get_byte (a, i), carry, tmp%bytes(i), carry)
end do
end if
call normalize (tmp)
c = tmp
end function big_subabs
function big_mul (a, b) result (c)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable :: c
!
! This is Knuth, Volume 2, Algorithm 4.3.1M.
!
integer :: na, nb, nc
integer :: i, j
integer :: ia, ib, ic
integer :: carry
type(big_integer), allocatable :: tmp
na = size (a%bytes)
nb = size (b%bytes)
nc = na + nb + 1
allocate (tmp)
allocate (tmp%bytes(nc))
tmp%bytes = zero
j = 1
do j = 1, nb
ib = byte2i (b%bytes(j))
if (ib /= 0) then
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = byte2i (tmp%bytes(i + j - 1))
ic = (ia * ib) + ic + carry
tmp%bytes(i + j - 1) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(na + j) = i2byte (carry)
end if
end do
tmp%sign = (a%sign .neqv. b%sign)
call normalize (tmp)
c = tmp
end function big_mul
subroutine big_divrem (a, b, q, r)
type(big_integer), intent(in) :: a, b
type(big_integer), allocatable, intent(inout) :: q, r
!
! Division with a remainder that is never negative. Equivalently,
! this is floor division if the divisor is positive, and ceiling
! division if the divisor is negative.
!
! See Raymond T. Boute, "The Euclidean definition of the functions
! div and mod", ACM Transactions on Programming Languages and
! Systems, Volume 14, Issue 2, pp. 127-144.
! https://doi.org/10.1145/128861.128862
!
call nonnegative_division (a, b, .true., .true., q, r)
if (a%sign) then
if (big_sgn (r) /= 0) then
q = short_addition (q, 1)
r = big_sub (big_abs (b), r)
end if
q%sign = .not. b%sign
else
q%sign = b%sign
end if
end subroutine big_divrem
function short_addition (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c
! Compute abs(a) + b.
integer :: na, nc, i
logical :: carry
type(big_integer), allocatable :: tmp
na = size (a%bytes)
nc = na + 1
allocate(tmp)
allocate(tmp%bytes(nc))
call add_bytes (a%bytes(1), i2byte (b), .false., tmp%bytes(1), carry)
do i = 2, na
call add_bytes (a%bytes(i), zero, carry, tmp%bytes(i), carry)
end do
tmp%bytes(nc) = logical2byte (carry)
call normalize (tmp)
c = tmp
end function short_addition
function short_multiplication (a, b) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), allocatable :: c
integer :: i, na, nc
integer :: ia, ic
integer :: carry
type(big_integer), allocatable :: tmp
na = size (a%bytes)
nc = na + 1
allocate (tmp)
allocate (tmp%bytes(nc))
tmp%sign = a%sign
carry = 0
do i = 1, na
ia = byte2i (a%bytes(i))
ic = (ia * b) + carry
tmp%bytes(i) = i2byte (iand (ic, 255))
carry = ishft (ic, -8)
end do
tmp%bytes(nc) = i2byte (carry)
call normalize (tmp)
c = tmp
end function short_multiplication
! Division without regard to signs.
subroutine nonnegative_division (a, b, want_q, want_r, q, r)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_q, want_r
type(big_integer), intent(inout), allocatable :: q, r
integer :: na, nb
integer :: remainder
na = size (a%bytes)
nb = size (b%bytes)
! It is an error if b has "significant" zero-bytes or is equal to
! zero.
if (b%bytes(nb) == zero) error stop
if (nb == 1) then
if (want_q) then
call short_division (a, byte2i (b%bytes(1)), q, remainder)
else
block
type(big_integer), allocatable :: bit_bucket
call short_division (a, byte2i (b%bytes(1)), bit_bucket, remainder)
end block
end if
if (want_r) then
if (allocated (r)) deallocate (r)
allocate (r)
allocate (r%bytes(1))
r%bytes(1) = i2byte (remainder)
end if
else
if (na >= nb) then
call long_division (a, b, want_q, want_r, q, r)
else
if (want_q) q = string2big ("0")
if (want_r) r = a
end if
end if
end subroutine nonnegative_division
subroutine short_division (a, b, q, r)
type(big_integer), intent(in) :: a
integer, intent(in) :: b
type(big_integer), intent(inout), allocatable :: q
integer, intent(inout) :: r
!
! This is Knuth, Volume 2, Exercise 4.3.1.16.
!
! The divisor is assumed to be positive.
!
integer :: n, i
integer :: ia, ib, iq
type(big_integer), allocatable :: tmp
ib = b
n = size (a%bytes)
allocate (tmp)
allocate (tmp%bytes(n))
r = 0
do i = n, 1, -1
ia = (256 * r) + byte2i (a%bytes(i))
iq = ia / ib
r = mod (ia, ib)
tmp%bytes(i) = i2byte (iq)
end do
tmp%sign = a%sign
call normalize (tmp)
q = tmp
end subroutine short_division
subroutine long_division (a, b, want_quotient, want_remainder, quotient, remainder)
type(big_integer), intent(in) :: a, b
logical, intent(in) :: want_quotient, want_remainder
type(big_integer), intent(inout), allocatable :: quotient
type(big_integer), intent(inout), allocatable :: remainder
!
! This is Knuth, Volume 2, Algorithm 4.3.1D.
!
! We do not deal here with the signs of the inputs and outputs.
!
! It is assumed size(a%bytes) >= size(b%bytes), and that b has no
! leading zero-bytes and is at least two bytes long. If b is one
! byte long and nonzero, use short division.
!
integer :: na, nb, m, n
integer :: num_lz, num_nonlz
integer :: j
integer :: qhat
logical :: carry
!
! We will NOT be working with VERY large numbers, and so it will
! be safe to put temporary storage on the stack. (Note: your
! Fortran might put this storage in a heap instead of the stack.)
!
! v = b, normalized to put its most significant 1-bit all the
! way left.
!
! u = a, shifted left by the same amount as b.
!
! q = the quotient.
!
! The remainder, although shifted left, will end up in u.
!
integer(bytekind) :: u(0:size (a%bytes) + size (b%bytes))
integer(bytekind) :: v(0:size (b%bytes) - 1)
integer(bytekind) :: q(0:size (a%bytes) - size (b%bytes))
na = size (a%bytes)
nb = size (b%bytes)
n = nb
m = na - nb
! In the most significant byte of the divisor, find the number of
! leading zero bits, and the number of bits after that.
block
integer(bytekind) :: tmp
tmp = byte2bk (b%bytes(n))
num_nonlz = bit_size (tmp) - leadz (tmp)
num_lz = 8 - num_nonlz
end block
call normalize_v (b%bytes) ! Make the most significant bit of v be one.
call normalize_u (a%bytes) ! Shifted by the same amount as v.
! Assure ourselves that the most significant bit of v is a one.
if (.not. btest (v(n - 1), 7)) error stop
do j = m, 0, -1
call calculate_qhat (qhat)
call multiply_and_subtract (carry)
q(j) = i2bk (qhat)
if (carry) call add_back
end do
if (want_quotient) then
if (allocated (quotient)) deallocate (quotient)
allocate (quotient)
allocate (quotient%bytes(m + 1))
quotient%bytes = bk2byte (q)
call normalize (quotient)
end if
if (want_remainder) then
if (allocated (remainder)) deallocate (remainder)
allocate (remainder)
allocate (remainder%bytes(n))
call unnormalize_u (remainder%bytes)
call normalize (remainder)
end if
contains
subroutine normalize_v (b_bytes)
character, intent(in) :: b_bytes(n)
!
! Normalize v so its most significant bit is a one. Any
! normalization factor that achieves this goal will suffice; we
! choose 2**num_lz. (Knuth uses (2**32) div (y[n-1] + 1).)
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!
integer :: i
integer(bytekind) :: btmp(0:n - 1)
btmp = byte2bk (b_bytes)
v(0) = lshftbk (btmp(0), num_lz)
do i = 1, n - 1
v(i) = ior (lshftbk (btmp(i), num_lz), &
& rshftbk (btmp(i - 1), num_nonlz))
end do
end subroutine normalize_v
subroutine normalize_u (a_bytes)
character, intent(in) :: a_bytes(m + n)
!
! Shift a leftwards to get u. Shift by as much as b was shifted
! to get v.
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!
integer :: i
integer(bytekind) :: atmp(0:m + n - 1)
atmp = byte2bk (a_bytes)
u(0) = lshftbk (atmp(0), num_lz)
do i = 1, m + n - 1
u(i) = ior (lshftbk (atmp(i), num_lz), &
& rshftbk (atmp(i - 1), num_nonlz))
end do
u(m + n) = rshftbk (atmp(m + n - 1), num_nonlz)
end subroutine normalize_u
subroutine unnormalize_u (r_bytes)
character, intent(out) :: r_bytes(n)
!
! Strictly for readability, we use linear stack space for an
! intermediate result.
!
integer :: i
integer(bytekind) :: rtmp(0:n - 1)
do i = 0, n - 1
rtmp(i) = ior (rshftbk (u(i), num_lz), &
& lshftbk (u(i + 1), num_nonlz))
end do
rtmp(n - 1) = rshftbk (u(n - 1), num_lz)
r_bytes = bk2byte (rtmp)
end subroutine unnormalize_u
subroutine calculate_qhat (qhat)
integer, intent(out) :: qhat
integer :: itmp, rhat
logical :: adjust
itmp = ior (lshfti (bk2i (u(j + n)), 8), &
& bk2i (u(j + n - 1)))
qhat = itmp / bk2i (v(n - 1))
rhat = mod (itmp, bk2i (v(n - 1)))
adjust = .true.
do while (adjust)
if (rshfti (qhat, 8) /= 0) then
continue
else if (qhat * bk2i (v(n - 2)) &
& > ior (lshfti (rhat, 8), &
& bk2i (u(j + n - 2)))) then
continue
else
adjust = .false.
end if
if (adjust) then
qhat = qhat - 1
rhat = rhat + bk2i (v(n - 1))
if (rshfti (rhat, 8) == 0) then
adjust = .false.
end if
end if
end do
end subroutine calculate_qhat
subroutine multiply_and_subtract (carry)
logical, intent(out) :: carry
integer :: i
integer :: qhat_v
integer :: mul_carry, sub_carry
integer :: diff
mul_carry = 0
sub_carry = 0
do i = 0, n
! Multiplication.
qhat_v = mul_carry
if (i /= n) qhat_v = qhat_v + (qhat * bk2i (v(i)))
mul_carry = rshfti (qhat_v, 8)
qhat_v = iand (qhat_v, 255)
! Subtraction.
diff = bk2i (u(j + i)) - qhat_v + sub_carry
sub_carry = -(logical2i (diff < 0)) ! Carry 0 or -1.
u(j + i) = i2bk (diff)
end do
carry = (sub_carry /= 0)
end subroutine multiply_and_subtract
subroutine add_back
integer :: i, carry, sum
q(j) = q(j) - 1_bytekind
carry = 0
do i = 0, n - 1
sum = bk2i (u(j + i)) + bk2i (v(i)) + carry
carry = ishft (sum, -8)
u(j + i) = i2bk (sum)
end do
end subroutine add_back
end subroutine long_division
function big_pow (a, i) result (c)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
type(big_integer), allocatable :: c
type(big_integer), allocatable :: base
integer :: exponent, exponent_halved
integer :: j, last_set
if (i < 0) error stop
if (i == 0) then
c = integer2big (1)
else
last_set = bit_size (i) - leadz (i)
base = a
exponent = i
c = integer2big (1)
do j = 0, last_set - 1
exponent_halved = exponent / 2
if (2 * exponent_halved /= exponent) then
c = c * base
end if
exponent = exponent_halved
base = base * base
end do
end if
end function big_pow
subroutine add_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out
integer :: ia, ib, ic
ia = byte2i (a)
if (carry_in) ia = ia + 1
ib = byte2i (b)
ic = ia + ib
c = i2byte (iand (ic, 255))
carry_out = (ic >= 256)
end subroutine add_bytes
subroutine sub_bytes (a, b, carry_in, c, carry_out)
character, intent(in) :: a, b
logical, value :: carry_in
character, intent(inout) :: c
logical, intent(inout) :: carry_out
integer :: ia, ib, ic
ia = byte2i (a)
ib = byte2i (b)
if (carry_in) ib = ib + 1
ic = ia - ib
carry_out = (ic < 0)
if (carry_out) ic = ic + 256
c = i2byte (iand (ic, 255))
end subroutine sub_bytes
function get_byte (a, i) result (byte)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
character :: byte
if (size (a%bytes) < i) then
byte = zero
else
byte = a%bytes(i)
end if
end function get_byte
function byteval (a, i) result (v)
type(big_integer), intent(in) :: a
integer, intent(in) :: i
integer :: v
if (size (a%bytes) < i) then
v = 0
else
v = byte2i (a%bytes(i))
end if
end function byteval
subroutine normalize (a)
type(big_integer), intent(inout) :: a
logical :: done
integer :: i
character, allocatable :: fewer_bytes(:)
! Shorten to the minimum number of bytes.
i = size (a%bytes)
done = .false.
do while (.not. done)
if (i == 1) then
done = .true.
else if (a%bytes(i) /= zero) then
done = .true.
else
i = i - 1
end if
end do
if (i /= size (a%bytes)) then
allocate (fewer_bytes (i))
fewer_bytes = a%bytes(1:i)
call move_alloc (fewer_bytes, a%bytes)
end if
! If the magnitude is zero, then clear the sign bit.
if (size (a%bytes) == 1) then
if (a%bytes(1) == zero) then
a%sign = .false.
end if
end if
end subroutine normalize
end module big_integers
program modular_exponentiation_task
use, non_intrinsic :: big_integers
implicit none
type(big_integer), allocatable :: zero, one, two
type(big_integer), allocatable :: a, b, modulus
zero = integer2big (0)
one = integer2big (1)
two = integer2big (2)
a = string2big ("2988348162058574136915891421498819466320163312926952423791023078876139")
b = string2big ("2351399303373464486466122544523690094744975233415544072992656881240319")
modulus = string2big ("10") ** 40
write (*,*) big2string (modular_pow (a, b, modulus))
contains
! The right-to-left binary method,
! https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
function modular_pow (base, exponent, modulus) result (retval)
type(big_integer), intent(in) :: base, exponent, modulus
type(big_integer), allocatable :: retval
type(big_integer), allocatable :: bas, expnt, remainder, bit_bucket
if (big_sgn (modulus - one) == 0) then
retval = zero
else
retval = one
bas = base
expnt = exponent
do while (big_sgn (expnt) /= 0)
call big_divrem (expnt, two, expnt, remainder)
if (big_sgn (remainder) /= 0) then
call big_divrem (retval * bas, modulus, bit_bucket, retval)
end if
call big_divrem (bas * bas, modulus, bit_bucket, bas)
end do
end if
end function modular_pow
end program modular_exponentiation_task
</syntaxhighlight>
{{out}}
<pre>$ gfortran -g -fbounds-check -Wall -Wextra modular_exponentiation.f90 && ./a.out
1527229998585248450016808958343740453059
</pre>
=={{header|FreeBASIC}}==
<syntaxhighlight lang="freebasic">
'From first principles (No external library)
Line 792 ⟶ 1,907:
</syntaxhighlight>
<pre>
Result:
Line 799 ⟶ 1,914:
Done
</pre>
=={{header|Frink}}==
<syntaxhighlight lang="frink">a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
println[modPow[a,b,10^40]]</syntaxhighlight>
{{out}}
<pre>1527229998585248450016808958343740453059
</pre>
=={{header|GAP}}==
<
a := 2988348162058574136915891421498819466320163312926952423791023078876139;
b := 2351399303373464486466122544523690094744975233415544072992656881240319;
Line 822 ⟶ 1,946:
end;
PowerModAlt(a, b, 10^40);</
=={{header|gnuplot}}==
<
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</
=={{header|Go}}==
<
import (
Line 851 ⟶ 1,975:
r.Exp(a, b, m)
fmt.Println(r)
}</
{{out}}
<pre>
Line 858 ⟶ 1,982:
=={{header|Groovy}}==
<
2351399303373464486466122544523690094744975233415544072992656881240319,
10000000000000000000000000000000000000000)</
Ouput:
<pre>1527229998585248450016808958343740453059</pre>
=={{header|Haskell}}==
<syntaxhighlight lang="haskell">
modPow :: Integer -> Integer -> Integer -> Integer -> Integer
modPow b 0 m r = r
modPow b e m r
| e `mod` 2 == 1 = modPow b' e' m (r * b `mod` m)
| otherwise = modPow b' e' m r
where
b' = b * b `mod` m
e' = e `div` 2
main
print (modPow 2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
(10 ^ 40)
1
</syntaxhighlight>
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
or in terms of ''until'':
<
powerMod b e m = x
where
Line 903 ⟶ 2,030:
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
(10 ^ 40)</
{{Out}}
<pre>1527229998585248450016808958343740453059</pre>
Line 909 ⟶ 2,036:
=={{header|Icon}} and {{header|Unicon}}==
This uses the exponentiation procedure from [[RSA_code#Icon_and_Unicon|RSA Code]] an example of the right to left binary method.
<
a := 2988348162058574136915891421498819466320163312926952423791023078876139
b := 2351399303373464486466122544523690094744975233415544072992656881240319
write("last 40 digits = ",mod_power(a,b,(10^40)))
end
Line 925 ⟶ 2,052:
}
return result
end</
{{out}}
Line 931 ⟶ 2,058:
=={{header|J}}==
'''Solution''':<syntaxhighlight lang
'''Example''':<
b =: 2351399303373464486466122544523690094744975233415544072992656881240319x
m =: 10^40x
a m&|@^ b
1527229998585248450016808958343740453059</
'''Discussion''': The phrase <tt>a m&|@^ b</tt> is the natural expression of <tt>a^b mod m</tt> in J, and is recognized by the interpreter as an opportunity for optimization, by [http://www.jsoftware.com/help/dictionary/special.htm#recognized%20phrase avoiding the exponentiation].
Line 943 ⟶ 2,070:
<code>java.math.BigInteger.modPow</code> solves this task. Inside [[OpenJDK]], [http://hg.openjdk.java.net/jdk7/jdk7/jdk/file/f097ca2434b1/src/share/classes/java/math/BigInteger.java BigInteger.java] implements <code>BigInteger.modPow</code> with a fast algorithm from [http://philzimmermann.com/EN/bnlib/index.html 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]].
<
public class PowMod {
Line 955 ⟶ 2,082:
System.out.println(a.modPow(b, m));
}
}</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|jq}}==
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq, and with fq.'''
'''Adapted from [[#Wren|Wren]]'''
<syntaxhighlight lang=jq>
# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
# Returns (. ^ $exp) % $mod
# where $exp >= 0, $mod != 0, and the input are integers.
def modPow($exp; $mod):
if $mod == 0 then "Cannot take modPow with modulus 0." | error
elif $exp < 0 then "modPow with exp < 0 is not supported." | error
else . as $x
| {r: 1, base: ($x % $mod), exp: $exp}
| until( .exp <= 0 or .emit;
if .base == 0 then .emit = 0
else if .exp%2 == 1
then .r = (.r * .base) % $mod
| .exp |= (. - 1) / 2
else .exp /= 2
end
| .base |= (. * .) % $mod
end )
| if .emit then .emit else .r end
end;
def task:
2988348162058574136915891421498819466320163312926952423791023078876139 as $a
| 2351399303373464486466122544523690094744975233415544072992656881240319 as $b
| (10|power(40)) as $m
| $a | modPow($b; $m) ;
task
</syntaxhighlight>
{{output}}
<pre>
1527229998585248450016808958343740453059
</pre>
=={{header|Julia}}==
Line 963 ⟶ 2,131:
We can use the built-in <code>powermod</code> function with the built-in <code>BigInt</code> type (based on GMP):
<
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = big(10) ^ 40
@show powermod(a, b, m)</
{{out}}
Line 972 ⟶ 2,140:
=={{header|Kotlin}}==
<
import java.math.BigInteger
Line 981 ⟶ 2,149:
val m = BigInteger.TEN.pow(40)
println(a.modPow(b, m))
}</
{{out}}
Line 992 ⟶ 2,160:
Following scheme
<
We will call the lib_BN library for big numbers:
Line 1,021 ⟶ 2,189:
{BN.pow 10 40}}
-> 1527229998585248450016808958343740453059 // 3300ms
</syntaxhighlight>
=={{header|Maple}}==
<
b := 2351399303373464486466122544523690094744975233415544072992656881240319:
a &^ b mod 10^40;</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
b = 2351399303373464486466122544523690094744975233415544072992656881240319;
m = 10^40;
PowerMod[a, b, m]
-> 1527229998585248450016808958343740453059</
=={{header|Maxima}}==
<
b: 2351399303373464486466122544523690094744975233415544072992656881240319$
power_mod(a, b, 10^40);
/* 1527229998585248450016808958343740453059 */</
=={{header|Nim}}==
{{libheader|bigints}}
<
proc powmod(b, e, m: BigInt): BigInt =
Line 1,062 ⟶ 2,230:
b = initBigInt("2351399303373464486466122544523690094744975233415544072992656881240319")
echo powmod(a, b, 10.pow 40)</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|ObjectIcon}}==
{{trans|Icon and Unicon}}
<syntaxhighlight lang="objecticon">
# -*- ObjectIcon -*-
#
# This program is close to being an exact copy of the Icon.
#
import io # <-- Object Icon requires this for I/O.
procedure main()
local a, b # <-- Object Icon forces you to declare your variables.
a := 2988348162058574136915891421498819466320163312926952423791023078876139
b := 2351399303373464486466122544523690094744975233415544072992656881240319
# You could leave out the "io." in the call to "write" below,
# because there is some "compatibility with regular Icon" support in
# the io package.
io.write("last 40 digits = ", mod_power(a,b,(10^40)))
end
procedure mod_power(base, exponent, modulus)
local result
result := 1
while exponent > 0 do
{
if exponent % 2 = 1 then
result := (result * base) % modulus
exponent /:= 2
base := base ^ 2 % modulus
}
return result
end
</syntaxhighlight>
<pre>$ oiscript modular-exponentiation-OI.icn
last 40 digits = 1527229998585248450016808958343740453059
</pre>
=={{header|OCaml}}==
Line 1,070 ⟶ 2,280:
Using the zarith library:
<
let a = Z.of_string "2988348162058574136915891421498819466320163312926952423791023078876139" in
let b = Z.of_string "2351399303373464486466122544523690094744975233415544072992656881240319" in
Line 1,076 ⟶ 2,286:
Z.powm a b m
|> Z.to_string
|> print_endline</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
Line 1,084 ⟶ 2,294:
Usage : a b mod powmod
<
1 exponent dup ifZero: [ return ]
while ( dup 0 > ) [
dup isEven ifFalse: [ swap base * modulus mod swap ]
2 / base sq modulus mod ->base
] drop ;</
{{out}}
Line 1,102 ⟶ 2,312:
1527229998585248450016808958343740453059
ok
</pre>
=={{header|PARI/GP}}==
<
b=2351399303373464486466122544523690094744975233415544072992656881240319;
lift(Mod(a,10^40)^b)</
=={{header|Pascal}}==
Line 1,141 ⟶ 2,323:
{{libheader|GMP}}
A port of the C example using gmp.
<
uses
Line 1,166 ⟶ 2,348:
mpz_clear(m);
mpz_clear(r);
end.</
{{out}}
<pre>% ./ModularExponentiation
Line 1,174 ⟶ 2,356:
=={{header|Perl}}==
Calculating the result both with an explicit algorithm (borrowed from Raku example) and with a built-in available when the <code>use bigint</code> pragma is invoked. Note that <code>bmodpow</code> modifies the base value (here <code>$a</code>) while <code>expmod</code> does not.
<
sub expmod {
Line 1,191 ⟶ 2,373:
print expmod($a, $b, $m), "\n";
print $a->bmodpow($b, $m), "\n";</
{{out}}
<pre>1527229998585248450016808958343740453059
Line 1,200 ⟶ 2,382:
Already builtin as mpz_powm, but here is an explicit version
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
Line 1,234 ⟶ 2,416:
<span style="color: #7060A8;">mpz_powm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">result</span><span style="color: #0000FF;">,</span><span style="color: #000000;">base</span><span style="color: #0000FF;">,</span><span style="color: #000000;">exponent</span><span style="color: #0000FF;">,</span><span style="color: #000000;">modulus</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">result</span><span style="color: #0000FF;">)</span>
<!--</
{{out}}
Line 1,243 ⟶ 2,425:
=={{header|PHP}}==
<
$a = '2988348162058574136915891421498819466320163312926952423791023078876139';
$b = '2351399303373464486466122544523690094744975233415544072992656881240319';
$m = '1' . str_repeat('0', 40);
echo bcpowmod($a, $b, $m), "\n";</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
Line 1,253 ⟶ 2,435:
=={{header|PicoLisp}}==
The following function is taken from "lib/rsa.l":
<
(let M 1
(loop
Line 1,260 ⟶ 2,442:
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )</
Test:
<
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
10000000000000000000000000000000000000000 )
-> 1527229998585248450016808958343740453059</
=={{header|Powershell}}==
{{works with|Powershell|7}}
<syntaxhighlight lang="powershell">Function Invoke-ModuloExponentiation ([BigInt]$Base, [BigInt]$Exponent, $Modulo) {
$Result = 1
$Base = $Base % $Modulo
If ($Base -eq 0) {return 0}
While ($Exponent -gt 0) {
If (($Exponent -band 1) -eq 1) {$Result = ($Result * $Base) % $Modulo}
$Exponent = $Exponent -shr 1
$Base = ($Base * $Base) % $Modulo
}
return ($Result % $Modulo)
}
$a = [BigInt]::Parse('2988348162058574136915891421498819466320163312926952423791023078876139')
$b = [BigInt]::Parse('2351399303373464486466122544523690094744975233415544072992656881240319')
$m = [BigInt]::Pow(10, 40)
Invoke-ModuloExponentiation -Base $a -Exponent $b -Modulo $m</syntaxhighlight>
{{out}}
<pre>
1527229998585248450016808958343740453059
</pre>
=={{header|Prolog}}==
{{works with|SWI Prolog}}
SWI Prolog has a built-in function named powm for this purpose.
<
A = 2988348162058574136915891421498819466320163312926952423791023078876139,
B = 2351399303373464486466122544523690094744975233415544072992656881240319,
M is 10 ** 40,
P is powm(A, B, M),
writeln(P).</
{{out}}
Line 1,284 ⟶ 2,492:
=={{header|Python}}==
<
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ** 40
print(pow(a, b, m))</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
<
" Without using builtin function "
x = 1
Line 1,307 ⟶ 2,515:
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ** 40
print(power_mod(a, b, m))</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
Line 1,313 ⟶ 2,521:
=={{header|Quackery}}==
<
[ dup while
dup 1 & if
Line 1,327 ⟶ 2,535:
2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
10 40 ** **mod echo</
{{out}}
Line 1,334 ⟶ 2,542:
=={{header|Racket}}==
<
#lang racket
(require math)
Line 1,341 ⟶ 2,549:
(define m (expt 10 40))
(modular-expt a b m)
</syntaxhighlight>
{{out}}
<pre>
Line 1,350 ⟶ 2,558:
(formerly Perl 6)
This is specced as a built-in, but here's an explicit version:
<syntaxhighlight lang="raku"
my $c = 1;
repeat while $b div= 2 {
Line 1,362 ⟶ 2,570:
2988348162058574136915891421498819466320163312926952423791023078876139,
2351399303373464486466122544523690094744975233415544072992656881240319,
10**40;</
{{out}}
<pre>1527229998585248450016808958343740453059</pre>
=={{header|REXX}}==
<syntaxhighlight lang="rexx">
/* REXX Modular exponentiation */
say powerMod(,
2988348162058574136915891421498819466320163312926952423791023078876139,,
2351399303373464486466122544523690094744975233415544072992656881240319,,
1e40)
return
powerMod: procedure
parse arg base, exponent, modulus
/* we need a numeric precision of twice the modulus size, */
/* the exponent size, or the base size, whichever is largest */
numeric digits max(2 * length(format(modulus, , , 0)),,
length(format(exponent, , , 0)), length(format(base, , , 0)))
result = 1
base = base // modulus
do while exponent > 0
result = result * base // modulus
end
return result
</syntaxhighlight>
{{out}}
<pre>
Line 1,447 ⟶ 2,609:
=={{header|Ruby}}==
===Built in since version 2.5.===
<
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10**40
puts a.pow(b, m)</
===Using OpenSSL standard library===
{{libheader|OpenSSL}}
<
a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = 10 ** 40
puts a.to_bn.mod_exp(b, m)</
===Written in Ruby===
<
def mod_exp(n, e, mod)
fail ArgumentError, 'negative exponent' if e < 0
Line 1,471 ⟶ 2,633:
prod
end
</syntaxhighlight>
=={{header|Rust}}==
<
num = "0.2.0"
*/
Line 1,521 ⟶ 2,683:
base %= &m;
}
}</
'''Test code:'''
<
let (a, b, num_digits) = (
"2988348162058574136915891421498819466320163312926952423791023078876139",
Line 1,544 ⟶ 2,706:
println!("The last {} digits of\n{}\nto the power of\n{}\nare:\n{}",
num_digits, a, b, result);
}</
{{out}}
<pre>The last 40 digits of
Line 1,554 ⟶ 2,716:
=={{header|Scala}}==
<
val a = BigInt(
Line 1,561 ⟶ 2,723:
"2351399303373464486466122544523690094744975233415544072992656881240319")
println(a.modPow(b, BigInt(10).pow(40)))</
=={{header|Scheme}}==
<
(define (square n)
(* n n))
Line 1,580 ⟶ 2,742:
(mod-exp 2988348162058574136915891421498819466320163312926952423791023078876139
2351399303373464486466122544523690094744975233415544072992656881240319
(expt 10 40)))</
{{out}}
Line 1,592 ⟶ 2,754:
[http://seed7.sourceforge.net/libraries/bigint.htm#modPow%28in_var_bigInteger,in_var_bigInteger,in_bigInteger%29 modPow],
which does modular exponentiation.
<
include "bigint.s7i";
Line 1,600 ⟶ 2,762:
2351399303373464486466122544523690094744975233415544072992656881240319_,
10_ ** 40));
end func;</
{{out}}
Line 1,608 ⟶ 2,770:
The library bigint.s7i defines modPow with:
<
in var bigInteger: exponent, in bigInteger: modulus) is func
result
Line 1,624 ⟶ 2,786:
end while;
end if;
end func;</
Original source: [http://seed7.sourceforge.net/algorith/math.htm#modPow]
Line 1,630 ⟶ 2,792:
=={{header|Sidef}}==
Built-in:
<
2988348162058574136915891421498819466320163312926952423791023078876139,
2351399303373464486466122544523690094744975233415544072992656881240319,
10**40)</
User-defined:
<
var c = 1
do {
Line 1,643 ⟶ 2,805:
} while (b //= 2)
c
}</
{{out}}
<pre>
Line 1,656 ⟶ 2,818:
AttaSwift's BigInt has a built-in modPow method, but here we define a generic modPow.
<
func modPow<T: BinaryInteger>(n: T, e: T, m: T) -> T {
Line 1,686 ⟶ 2,848:
let b = BigInt("2351399303373464486466122544523690094744975233415544072992656881240319")
print(modPow(n: a, e: b, m: BigInt(10).power(40)))</
{{out}}
Line 1,698 ⟶ 2,860:
===Recursive===
<
# Algorithm from http://introcs.cs.princeton.edu/java/78crypto/ModExp.java.html
Line 1,710 ⟶ 2,872:
}
return $c
}</
Demonstrating:
<
set b 2351399303373464486466122544523690094744975233415544072992656881240319
set n [expr {10**40}]
puts [expr {modexp($a,$b,$n)}]</
{{out}}
<pre>
Line 1,722 ⟶ 2,884:
===Iterative===
<
proc modexp {a b n} {
for {set c 1} {$b} {set a [expr {$a*$a % $n}]} {
Line 1,731 ⟶ 2,893:
}
return $c
}</
Demonstrating:
<
set b 2351399303373464486466122544523690094744975233415544072992656881240319
set n [expr {10**40}]
puts [modexp $a $b $n]</
{{out}}
<pre>
Line 1,746 ⟶ 2,908:
From your system prompt:
<
2351399303373464486466122544523690094744975233415544072992656881240319
(expt 10 40)))'
1527229998585248450016808958343740453059</
=={{header|Visual Basic .NET}}==
{{works with|Visual Basic .NET|2011}}
<
Imports System.Numerics
Line 1,776 ⟶ 2,938:
Loop
Return result
End Function 'ModPowBig</
{{out}}
<pre>
Line 1,784 ⟶ 2,946:
=={{header|Wren}}==
{{libheader|Wren-big}}
<
var a = BigInt.new("2988348162058574136915891421498819466320163312926952423791023078876139")
var b = BigInt.new("2351399303373464486466122544523690094744975233415544072992656881240319")
var m = BigInt.ten.pow(40)
System.print(a.modPow(b, m))</
{{out}}
Line 1,798 ⟶ 2,960:
=={{header|zkl}}==
Using the GMP big num library:
<
a:=BN("2988348162058574136915891421498819466320163312926952423791023078876139");
b:=BN("2351399303373464486466122544523690094744975233415544072992656881240319");
m:=BN(10).pow(40);
a.powm(b,m).println();
a.powm(b,m) : "%,d".fmt(_).println();</
{{out}}
<pre>
Line 1,809 ⟶ 2,971:
1,527,229,998,585,248,450,016,808,958,343,740,453,059
</pre>
=={{header|zsh}}==
<syntaxhighlight lang="zsh">
#!/bin/zsh
#
# I use expr from GNU coreutils.
EXPR=`which expr`
#
# One could use GNU bc or any other such program that is capable of
# handling the large integers.
#
unset PATH # Besides shell commands, the script uses ONLY
# the expression calculator.
mod_power()
{
local base="${1}"
local exponent="${2}"
local modulus="${3}"
local result=1
while [[ `${EXPR} ${exponent} '>' 0` != 0 ]] ; do
if [[ `${EXPR} '(' ${exponent} '%' 2 ')' '=' 1` != 0 ]] ; then
result=`${EXPR} '(' ${result} '*' ${base} ')' '%' ${modulus}`
fi
exponent=`${EXPR} ${exponent} / 2`
base=`${EXPR} '(' ${base} '*' ${base} ')' '%' ${modulus}`
done
echo ${result}
}
a=2988348162058574136915891421498819466320163312926952423791023078876139
b=2351399303373464486466122544523690094744975233415544072992656881240319
# One followed by 40 zeros.
modulus=10000000000000000000000000000000000000000
# Or the following will work, as long as the modulus is at least
# 10**40. :)))))
#modulus=$(mod_power 10 40 1000000000000000000000000000000000000000000000000000000000000)
echo $(mod_power ${a} ${b} ${modulus})
exit 0
</syntaxhighlight>
{{out}}
<pre>$ zsh ./modular-exponentiation.zsh
1527229998585248450016808958343740453059</pre>
| |||