Modular inverse
You are encouraged to solve this task according to the task description, using any language you may know.
From Wikipedia:
- In modular arithmetic, the modular multiplicative inverse of an integer a modulo m is an integer x such that
Or in other words, such that:
It can be shown that such an inverse exists if and only if a and m are coprime, but we will ignore this for this task.
Either by implementing the algorithm, by using a dedicated library or by using a builtin function in your language, compute the modular inverse of 42 modulo 2017.
Ada
<lang Ada>with Ada.Text_IO;
procedure Mod_Inv is
procedure X_GCD(A, B: in Natural; D, X, Y: out Integer) is
-- the Extended Euclidean Algorithm
-- finds (D, X, Y) with D = GCD(A, B) = A*X + B*Y
R: Natural := A mod B;
begin
if R=0 then
D := B;
X := 0;
Y := 1;
else
X_GCD(B, R, D, Y, X);
Y := Y - (A/B)*X;
end if;
end X_GCD;
function Inverse(A, M: Integer) return Integer is
-- computes the multiplicative inverse of A mod M, using X_GCD
Result, GCD, Dummy: Integer;
begin
X_GCD(A, M, GCD, Result, Dummy);
if GCD /= 1 then -- inverse does not exist!
raise Constraint_Error with
"GCD (" & Integer'Image(A) & "," & Integer'Image(M) & " ) =" &
Integer'Image(GCD) & " /= 1";
else -- make sure Result is in {0, ..., M-1}
if Result < 0 then
return Result+M;
else
return Result;
end if;
end if;
end Inverse;
begin
Ada.Text_IO.Put_Line(Natural'Image(Inverse(42, 2017))); -- Ada.Text_IO.Put_Line(Natural'Image(Inverse(154, 3311))); -- The above would raise CONSTRAINT_ERROR : GCD ( 154, 3311 ) = 77 /= 1
end Mod_Inv;</lang>
AutoHotkey
Translation of C. <lang AutoHotkey>MsgBox, % ModInv(42, 2017)
ModInv(a, b) { if (b = 1) return 1 b0 := b, x0 := 0, x1 :=1 while (a > 1) { q := a // b , t := b , b := Mod(a, b) , a := t , t := x0 , x0 := x1 - q * x0 , x1 := t } if (x1 < 0) x1 += b0 return x1 }</lang> Output:
1969
C
<lang c>#include <stdio.h>
int mul_inv(int a, int b) { int b0 = b, t, q; int x0 = 0, x1 = 1; if (b == 1) return 1; while (a > 1) { q = a / b; t = b, b = a % b, a = t; t = x0, x0 = x1 - q * x0, x1 = t; } if (x1 < 0) x1 += b0; return x1; }
int main(void) { printf("%d\n", mul_inv(42, 2017)); return 0; }</lang>
C++
<lang cpp>#include <iostream>
using namespace std;
int mul_inv(int a, int b) { int b0 = b, t, q; int x0 = 0, x1 = 1; if (b == 1) return 1; while (a > 1) { q = a / b; t = b, b = a % b, a = t; t = x0, x0 = x1 - q * x0, x1 = t; } if (x1 < 0) x1 += b0; return x1; }
int main(void) { cout<<mul_inv(42, 2017)<<endl; return 0; } </lang>
D
<lang d>T modInverse(T)(T a, T b) pure nothrow {
if (b == 1)
return 1;
T b0 = b,
x0 = 0,
x1 = 1;
while (a > 1) {
immutable q = a / b;
auto t = b;
b = a % b;
a = t;
t = x0;
x0 = x1 - q * x0;
x1 = t;
}
return (x1 < 0) ? (x1 + b0) : x1;
}
void main() {
import std.stdio; writeln(modInverse(42, 2017));
}</lang>
- Output:
1969
Go
The standard library function uses the extended Euclidean algorithm internally. <lang go>package main
import ( "fmt" "math/big" )
func main() { a := big.NewInt(42) m := big.NewInt(2017) k := new(big.Int).ModInverse(a, m) fmt.Println(k) }</lang>
- Output:
1969
Haskell
<lang haskell>-- Extended Euclidean algorithm. Given non-negative a and b, return x, y and g -- such that ax + by = g, where g = gcd(a,b). Note that x or y may be negative. gcdExt a 0 = (1, 0, a) gcdExt a b = let (q, r) = a `quotRem` b
(s, t, g) = gcdExt b r
in (t, s - q * t, g)
-- Given a and m, return Just x such that ax = 1 mod m. If there is no such x -- return Nothing. modInv a m = let (i, _, g) = gcdExt a m
in if g == 1 then Just (mkPos i) else Nothing where mkPos x = if x < 0 then x + m else x
main = do
print $ 2 `modInv` 4 print $ 42 `modInv` 2017</lang>
Output
Nothing Just 1969
Icon and Unicon
<lang unicon>procedure main(args)
a := integer(args[1]) | 42 b := integer(args[2]) | 2017 write(mul_inv(a,b))
end
procedure mul_inv(a,b)
if b == 1 then return 1
(b0 := b, x0 := 0, x1 := 1)
while a > 1 do {
q := a/b
(t := b, b := a%b, a := t)
(t := x0, x0 := x1-q*x0, x1 := t)
}
return if (x1 > 0) then x1 else x1+b0
end</lang>
Sample output:
->mi 1969 ->
Adding a coprime test:
<lang unicon>link numbers
procedure main(args)
a := integer(args[1]) | 42 b := integer(args[2]) | 2017 write(mul_inv(a,b))
end
procedure mul_inv(a,b)
if b == 1 then return 1
if gcd(a,b) ~= 1 then return "not coprime"
(b0 := b, x0 := 0, x1 := 1)
while a > 1 do {
q := a/b
(t := b, b := a%b, a := t)
(t := x0, x0 := x1-q*x0, x1 := t)
}
return if (x1 > 0) then x1 else x1+b0
end</lang>
J
Solution:<lang j> modInv =: dyad def 'x y&|@^ <: 5 p: y'"0</lang> Example:<lang j> 42 modInv 2017 1969</lang> Notes: Calculates the modular inverse as a^( totient(m) - 1 ) mod m. Note that J has a fast implementation of modular exponentiation (which avoids the exponentiation altogether), invoked with the form m&|@^ (hence, we use explicitly-named arguments for this entry, as opposed to the "variable free" tacit style).
Java
The BigInteger library has a method for this:
<lang java>System.out.println(BigInteger.valueOf(42).modInverse(BigInteger.valueOf(2017)));</lang>
- Output:
1969
Julia
Built-in
Julia includes a built-in function for this: <lang julia>invmod(a, b)</lang>
C translation
The following code works on any integer type. To maximize performance, we ensure (via a promotion rule) that the operands are the same type (and use built-ins zero(T) and one(T) for initialization of temporary variables to ensure that they remain of the same type throughout execution).
<lang julia>function modinv{T<:Integer}(a::T, b::T)
b0 = b
x0, x1 = zero(T), one(T)
while a > 1
q = div(a, b)
a, b = b, a % b
x0, x1 = x1 - q * x0, x0
end
x1 < 0 ? x1 + b0 : x1
end modinv(a::Integer, b::Integer) = modinv(promote(a,b)...)</lang>
Output
julia> invmod(42, 2017) 1969 julia> modinv(42, 2017) 1969
Mathematica
The built-in function FindInstance works well for this
<lang Mathematica>modInv[a_, m_] :=
Block[{x,k}, x /. FindInstance[a x == 1 + k m, {x, k}, Integers]]</lang>
modInv[42, 2017]
{1969}
МК-61/52
<lang>П1 П2 <-> П0 0 П5 1 П6 ИП1 1 - x=0 14 С/П ИП0 1 - /-/ x<0 50 ИП0 ИП1 / [x] П4 ИП1 П3 ИП0 ^ ИП1 / [x] ИП1 * - П1 ИП3 П0 ИП5 П3 ИП6 ИП4 ИП5 * - П5 ИП3 П6 БП 14 ИП6 x<0 55 ИП2 + С/П</lang>
OCaml
<lang ocaml>let mul_inv a = function 1 -> 1 | b ->
let rec aux a b x0 x1 = if a <= 1 then x1 else if b = 0 then failwith "mul_inv" else aux b (a mod b) (x1 - (a / b) * x0) x0 in let x = aux a b 0 1 in if x < 0 then x + b else x</lang>
Testing:
# mul_inv 42 2017 ;; - : int = 1969
<lang ocaml>let rec gcd_ext a = function
| 0 -> (1, 0, a)
| b ->
let s, t, g = gcd_ext b (a mod b) in
(t, s - (a / b) * t, g)
let mod_inv a m =
let mk_pos x = if x < 0 then x + m else x in match gcd_ext a m with | i, _, 1 -> mk_pos i | _ -> failwith "mod_inv"</lang>
Testing:
# mod_inv 42 2017 ;; - : int = 1969
PARI/GP
<lang parigp>Mod(1/42,2017)</lang>
Perl
The modular inverse is not a perl builtin but there is a CPAN module who does the job.
<lang perl>use Math::ModInt qw(mod); print mod(42, 2017)->inverse</lang>
- Output:
mod(1969, 2017)
Perl 6
<lang perl6>sub inverse($n, :$modulo) {
my ($c, $d, $uc, $vc, $ud, $vd) = ($n % $modulo, $modulo, 1, 0, 0, 1);
my $q;
while $c != 0 {
($q, $c, $d) = ($d div $c, $d % $c, $c);
($uc, $vc, $ud, $vd) = ($ud - $q*$uc, $vd - $q*$vc, $uc, $vc);
}
return $ud % $modulo;
}
say inverse 42, :modulo(2017)</lang>
Python
Implementation of this pseudocode with this. <lang python>>>> def extended_gcd(aa, bb):
lastremainder, remainder = abs(aa), abs(bb)
x, lastx, y, lasty = 0, 1, 1, 0
while remainder:
lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
x, lastx = lastx - quotient*x, x
y, lasty = lasty - quotient*y, y
return lastremainder, lastx * (-1 if aa < 0 else 1), lasty * (-1 if bb < 0 else 1)
>>> def modinv(a, m): g, x, y = extended_gcd(a, m) if g != 1: raise ValueError return x % m
>>> modinv(42, 2017) 1969 >>> </lang>
Racket
<lang racket> (require math) (modular-inverse 42 2017) </lang> Output: <lang racket> 1969 </lang>
REXX
<lang rexx>/*REXX program calcuates the modular inverse of an integer X modulo Y.*/ parse arg x y . /*get two integers from the C.L. */ say 'modular inverse of ' x " by " y ' ───► ' modInv(x,y) exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────MODINV subroutine───────────────────*/ modInv: parse arg a,b 1 ob; ox=0; $=1
if b \= 1 then do while a>1
parse value a/b a//b b ox with q b a t
ox=$-q*ox; $=trunc(t)
end /*while a>1*/
if $<0 then $=$+ob return $</lang> output when using the input of: 42 2017
modular inverse of 42 by 2017 ───► 1969
Ruby
<lang ruby>#based on pseudo code from http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Iterative_method_2 and from translating the python implementation. def extended_gcd(a, b)
last_remainder, remainder = a.abs, b.abs x, last_x, y, last_y = 0, 1, 1, 0 while remainder != 0 last_remainder, (quotient, remainder) = remainder, last_remainder.divmod(remainder) x, last_x = last_x - quotient*x, x y, last_y = last_y - quotient*y, y end
return last_remainder, last_x * (a < 0 ? -1 : 1)
end
def invmod(e, et)
g, x = extended_gcd(e, et) if g != 1 raise 'Teh maths are broken!' end x % et
end</lang>
?> invmod(42,2017) => 1969
Run BASIC
<lang runbasic>print multInv(42, 2017) end
function multInv(a,b) b0 = b multInv = 1 if b = 1 then goto [endFun] while a > 1 q = a / b t = b b = a mod b a = t t = x0 x0 = multInv - q * x0 multInv = int(t) wend if multInv < 0 then multInv = multInv + b0 [endFun] end function</lang>
- Output:
1969
Seed7
The library bigint.s7i defines the bigInteger function modInverse. It returns the modular multiplicative inverse of a modulo b when a and b are coprime (gcd(a, b) = 1). If a and b are not coprime (gcd(a, b) <> 1) the exception RANGE_ERROR is raised.
<lang seed7>const func bigInteger: modInverse (in var bigInteger: a,
in var bigInteger: b) is func
result
var bigInteger: modularInverse is 0_;
local
var bigInteger: b_bak is 0_;
var bigInteger: x is 0_;
var bigInteger: y is 1_;
var bigInteger: lastx is 1_;
var bigInteger: lasty is 0_;
var bigInteger: temp is 0_;
var bigInteger: quotient is 0_;
begin
if b < 0_ then
raise RANGE_ERROR;
end if;
if a < 0_ and b <> 0_ then
a := a mod b;
end if;
b_bak := b;
while b <> 0_ do
temp := b;
quotient := a div b;
b := a rem b;
a := temp;
temp := x;
x := lastx - quotient * x;
lastx := temp;
temp := y;
y := lasty - quotient * y;
lasty := temp;
end while;
if a = 1_ then
modularInverse := lastx;
if modularInverse < 0_ then
modularInverse +:= b_bak;
end if;
else
raise RANGE_ERROR;
end if;
end func;</lang>
Original source: [1]
Tcl
<lang tcl>proc gcdExt {a b} {
if {$b == 0} {
return [list 1 0 $a]
}
set q [expr {$a / $b}]
set r [expr {$a % $b}]
lassign [gcdExt $b $r] s t g
return [list $t [expr {$s - $q*$t}] $g]
} proc modInv {a m} {
lassign [gcdExt $a $m] i -> g
if {$g != 1} {
return -code error "no inverse exists of $a %! $m"
}
while {$i < 0} {incr i $m}
return $i
}</lang> Demonstrating <lang tcl>puts "42 %! 2017 = [modInv 42 2017]" catch {
puts "2 %! 4 = [modInv 2 4]"
} msg; puts $msg</lang>
- Output:
42 %! 2017 = 1969 no inverse exists of 2 %! 4
XPL0
<lang XPL0>code IntOut=11, Text=12; int X; def A=42, M=2017; [for X:= 2 to M-1 do
if rem(A*X/M) = 1 then [IntOut(0, X); exit];
Text(0, "Does not exist"); ]</lang>
- Output:
1969