Modular inverse

From Rosetta Code
Task
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.


Task

Either by implementing the algorithm, by using a dedicated library or by using a built-in function in your language,   compute the modular inverse of   42 modulo 2017.

8th[edit]

 
\ return "extended gcd" of a and b; The result satisfies the equation:
\ a*x + b*y = gcd(a,b)
: n:xgcd \ a b -- gcd x y
dup 0 n:= if
1 swap \ -- a 1 0
else
tuck n:/mod
-rot recurse
tuck 4 roll
n:* n:neg n:+
then ;
 
\ Return modular inverse of n modulo mod, or null if it doesn't exist (n and mod
\ not coprime):
: n:invmod \ n mod -- invmod
dup >r
n:xgcd rot 1 n:= not if
2drop null
else
drop dup 0 n:< if r@ n:+ then
then
rdrop ;
 
42 2017 n:invmod . cr bye
 
Output:
1969

Ada[edit]

We start with specifying a package "Mod_Inv". The same package is used to solve the Chinese Remainder task [[1]].

package Mod_Inv is
 
procedure X_GCD(A, B: in Natural; D, X, Y: out Integer);
-- the Extended Euclidean Algorithm
-- finds (D, X, Y) with D = GCD(A, B) = A*X + B*Y
 
function Inverse(A, M: Integer) return Integer;
-- computes the multiplicative inverse Inv_A of A mod M, using X_GCD
-- raises Constraint_Error if Inv_A does not exist
 
end Mod_Inv;

The implementation of "Mod_Inv" is as follows.

package body 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;
 
end Mod_Inv;

Now, the main program is a triviality:

with Ada.Text_IO; with Mod_Inv; use Mod_Inv, Ada.text_IO;
 
procedure Mod_Inv_Test is
begin
-- Put_Line(Natural'Image(Inverse(154, 3311)));
-- The above would raise CONSTRAINT_ERROR : GCD ( 154, 3311 ) = 77 /= 1
 
Put_Line(Natural'Image(Inverse(42, 2017)));
end Mod_Inv_Test;

ALGOL 68[edit]

 
BEGIN
PROC modular inverse = (INT a, m) INT :
BEGIN
PROC extended gcd = (INT x, y) []INT :
CO
Algol 68 allows us to return three INTs in several ways. A [3]INT
is used here but it could just as well be a STRUCT.
CO
BEGIN
INT v := 1, a := 1, u := 0, b := 0, g := x, w := y;
WHILE w>0
DO
INT q := g % w, t := a - q * u;
a := u; u := t;
t := b - q * v;
b := v; v := t;
t := g - q * w;
g := w; w := t
OD;
a PLUSAB (a < 0 | u | 0);
(a, b, g)
END;
[] INT egcd = extended gcd (a, m);
(egcd[3] > 1 | 0 | egcd[1] MOD m)
END;
printf (($"42 ^ -1 (mod 2017) = ", g(0)$, modular inverse (42, 2017)))
CO
Note that if ϕ(m) is known, then a^-1 = a^(ϕ(m)-1) mod m which
allows an alternative implementation in terms of modular
exponentiation but, in general, this requires the factorization of
m. If m is prime the factorization is trivial and ϕ(m) = m-1.
2017
is prime which may, or may not, be ironic within the context
of the Rosetta Code conditions.
CO
END
 
Output:
42 ^ -1 (mod 2017) = 1969

AutoHotkey[edit]

Translation of C.

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
}
Output:
1969

AWK[edit]

 
# syntax: GAWK -f MODULAR_INVERSE.AWK
# converted from C
BEGIN {
printf("%s\n",mod_inv(42,2017))
exit(0)
}
function mod_inv(a,b, b0,t,q,x0,x1) {
b0 = b
x0 = 0
x1 = 1
if (b == 1) {
return(1)
}
while (a > 1) {
q = int(a / b)
t = b
b = int(a % b)
a = t
t = x0
x0 = x1 - q * x0
x1 = t
}
if (x1 < 0) {
x1 += b0
}
return(x1)
}
 
Output:
1969

Batch File[edit]

Based from C's second implementation

Translation of: Perl
@echo off
setlocal enabledelayedexpansion
%== Calls the "function" ==%
call :ModInv 42 2017 result
echo !result!
call :ModInv 40 1 result
echo !result!
call :ModInv 52 -217 result
echo !result!
call :ModInv -486 217 result
echo !result!
call :ModInv 40 2018 result
echo !result!
pause>nul
exit /b 0
 
%== The "function" ==%
:ModInv
set a=%1
set b=%2
 
if !b! lss 0 (set /a b=-b)
if !a! lss 0 (set /a a=b - ^(-a %% b^))
 
set t=0&set nt=1&set r=!b!&set /a nr=a%%b
 
:while_loop
if !nr! neq 0 (
set /a q=r/nr
set /a tmp=nt
set /a nt=t - ^(q*nt^)
set /a t=tmp
 
set /a tmp=nr
set /a nr=r - ^(q*nr^)
set /a r=tmp
goto while_loop
)
 
if !r! gtr 1 (set %3=-1&goto :EOF)
if !t! lss 0 set /a t+=b
set %3=!t!
goto :EOF
Output:
1969
0
96
121
-1

Bracmat[edit]

Translation of: Julia
( ( mod-inv
= a b b0 x0 x1 q
.  !arg:(?a.?b)
& ( !b:1
| (!b.0.1):(?b0.?x0.?x1)
& whl
' ( !a:>1
& div$(!a.!b):?q
& (!b.mod$(!a.!b)):(?a.?b)
& (!x1+-1*!q*!x0.!x0):(?x0.?x1)
)
& (!x:>0|!x1+!b0)
)
)
& out$(mod-inv$(42.2017))
};

Output

1969

C[edit]

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

The above method has some problems. Most importantly, when given a pair (a,b) with no solution, it generates an FP exception. When given b=1, it returns 1 which is not a valid result mod 1. When given negative a or b the results are incorrect. The following generates results that should match Pari/GP for numbers in the int range.

Translation of: Perl
#include <stdio.h>
 
int mul_inv(int a, int b)
{
int t, nt, r, nr, q, tmp;
if (b < 0) b = -b;
if (a < 0) a = b - (-a % b);
t = 0; nt = 1; r = b; nr = a % b;
while (nr != 0) {
q = r/nr;
tmp = nt; nt = t - q*nt; t = tmp;
tmp = nr; nr = r - q*nr; r = tmp;
}
if (r > 1) return -1; /* No inverse */
if (t < 0) t += b;
return t;
}
int main(void) {
printf("%d\n", mul_inv(42, 2017));
printf("%d\n", mul_inv(40, 1));
printf("%d\n", mul_inv(52, -217)); /* Pari semantics for negative modulus */
printf("%d\n", mul_inv(-486, 217));
printf("%d\n", mul_inv(40, 2018));
return 0;
}
Output:
1969
0
96
121
-1

C++[edit]

Translation of: C
#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;
}
 

Recursive implementation

#include <iostream>
 
short ObtainMultiplicativeInverse(int a, int b, int s0 = 1, int s1 = 0)
{
return b==0? s0: ObtainMultiplicativeInverse(b, a%b, s1, s0 - s1*(a/b));
}
 
int main(int argc, char* argv[])
{
std::cout << ObtainMultiplicativeInverse(42, 2017) << std::endl;
return 0;
}
 

Common Lisp[edit]

 
;;
;; Calculates the GCD of a and b based on the Extended Euclidean Algorithm. The function also returns
;; the Bézout coefficients s and t, such that gcd(a, b) = as + bt.
;;
;; The algorithm is described on page http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Iterative_method_2
;;
(defun egcd (a b)
(do ((r (cons b a) (cons (- (cdr r) (* (car r) q)) (car r))) ; (r+1 r) i.e. the latest is first.
(s (cons 0 1) (cons (- (cdr s) (* (car s) q)) (car s))) ; (s+1 s)
(u (cons 1 0) (cons (- (cdr u) (* (car u) q)) (car u))) ; (t+1 t)
(q nil))
((zerop (car r)) (values (cdr r) (cdr s) (cdr u))) ; exit when r+1 = 0 and return r s t
(setq q (floor (/ (cdr r) (car r)))))) ; inside loop; calculate the q
 
;;
;; Calculates the inverse module for a = 1 (mod m).
;;
;; Note: The inverse is only defined when a and m are coprimes, i.e. gcd(a, m) = 1.”
;;
(defun invmod (a m)
(multiple-value-bind (r s k) (egcd a m)
(unless (= 1 r) (error "invmod: Values ~a and ~a are not coprimes." a m))
s))
 
Output:
* (invmod 42 2017)

-48
* (mod -48 2017)

1969

D[edit]

Translation of: C
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));
}
Output:
1969

EchoLisp[edit]

 
(lib 'math) ;; for egcd = extended gcd
 
(define (mod-inv x m)
(define-values (g inv q) (egcd x m))
(unless (= 1 g) (error 'not-coprimes (list x m) ))
(if (< inv 0) (+ m inv) inv))
 
(mod-inv 42 2017)1969
(mod-inv 42 666)
🔴 error: not-coprimes (42 666)
 

Elixir[edit]

Translation of: Ruby
defmodule Modular do
def extended_gcd(a, b) do
{last_remainder, last_x} = extended_gcd(abs(a), abs(b), 1, 0, 0, 1)
{last_remainder, last_x * (if a < 0, do: -1, else: 1)}
end
 
defp extended_gcd(last_remainder, 0, last_x, _, _, _), do: {last_remainder, last_x}
defp extended_gcd(last_remainder, remainder, last_x, x, last_y, y) do
quotient = div(last_remainder, remainder)
remainder2 = rem(last_remainder, remainder)
extended_gcd(remainder, remainder2, x, last_x - quotient*x, y, last_y - quotient*y)
end
 
def inverse(e, et) do
{g, x} = extended_gcd(e, et)
if g != 1, do: raise "The maths are broken!"
rem(x+et, et)
end
end
 
IO.puts Modular.inverse(42,2017)
Output:
1969

ERRE[edit]

PROGRAM MOD_INV
 
!$INTEGER
 
PROCEDURE MUL_INV(A,B->T)
LOCAL NT,R,NR,Q,TMP
IF B<0 THEN B=-B
IF A<0 THEN A=B-(-A MOD B)
T=0 NT=1 R=B NR=A MOD B
WHILE NR<>0 DO
Q=R DIV NR
TMP=NT NT=T-Q*NT T=TMP
TMP=NR NR=R-Q*NR R=TMP
END WHILE
IF (R>1) THEN T=-1 EXIT PROCEDURE  ! NO INVERSE
IF (T<0) THEN T+=B
END PROCEDURE
 
 
BEGIN
MUL_INV(42,2017->T) PRINT(T)
MUL_INV(40,1->T) PRINT(T)
MUL_INV(52,-217->T) PRINT(T)  ! pari semantics for negative modulus
MUL_INV(-486,217->T) PRINT(T)
MUL_INV(40,2018->T) PRINT(T)
END PROGRAM
 
Output:
 1969
 0
 96
 121
-1

FunL[edit]

import integers.egcd
 
def modinv( a, m ) =
val (g, x, _) = egcd( a, m )
 
if g != 1 then error( a + ' and ' + m + ' not coprime' )
 
val res = x % m
 
if res < 0 then res + m else res
 
println( modinv(42, 2017) )
Output:
1969

Forth[edit]

ANS Forth with double-number word set

 
: invmod { a m | v b c -- inv }
m to v
1 to c
0 to b
begin a
while v a / >r
c b s>d c s>d r@ 1 m*/ d- d>s to c to b
a v s>d a s>d r> 1 m*/ d- d>s to a to v
repeat b m mod dup to b 0<
if m b + else b then ;
 
42 2017 invmod . 1969 

Go[edit]

The standard library function uses the extended Euclidean algorithm internally.

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)
}
Output:
1969

Haskell[edit]

-- 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
Output:
Nothing
Just 1969

Icon and Unicon[edit]

Translation of: C
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
Output:
->mi
1969
->

Adding a coprime test:

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

J[edit]

Solution:
   modInv =: dyad def 'x y&|@^ <: 5 p: y'"0
Example:
   42 modInv 2017
1969

Notes:

  • Calculates the modular inverse as a^( totient(m) - 1 ) mod m.
  • 5 p: y is Euler's totient function of y.
  • 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: the m&| construct must freeze the value before it can be used but we want to use different values in that expression at different times...).

Java[edit]

The BigInteger library has a method for this:

System.out.println(BigInteger.valueOf(42).modInverse(BigInteger.valueOf(2017)));
Output:
1969

Julia[edit]

Built-in[edit]

Julia includes a built-in function for this:

invmod(a, b)

C translation[edit]

Translation of: C

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

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)...)
Output:
julia> invmod(42, 2017)
1969

julia> modinv(42, 2017)
1969

Maple[edit]

 
1/42 mod 2017;
 
Output:
                                    1969

Mathematica[edit]

The built-in function FindInstance works well for this

modInv[a_, m_] := 
Block[{x,k}, x /. FindInstance[a x == 1 + k m, {x, k}, Integers]]

Another way by using the built-in function PowerMod :

PowerMod[a,-1,m]

For example :

modInv[42, 2017]

{1969}

PowerMod[42, -1, 2017]

1969

МК-61/52[edit]

П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 + С/П

Nim[edit]

Translation of: C
proc mulInv(a0, b0): int =
var (a, b, x0) = (a0, b0, 0)
result = 1
if b == 1: return
while a > 1:
let q = a div b
a = a mod b
swap a, b
result = result - q * x0
swap x0, result
if result < 0: result += b0
 
echo mulInv(42, 2017)

OCaml[edit]

Translation of: C
[edit]

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

Testing:

# mul_inv 42 2017 ;;
- : int = 1969

Translation of: Haskell
[edit]

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"

Testing:

# mod_inv 42 2017 ;;
- : int = 1969

Oforth[edit]

Usage : a modulus invmod

// euclid ( a b -- u v r )
// Return r = gcd(a, b) and (u, v) / r = au + bv
 
: euclid(a, b)
| q u u1 v v1 |
 
b 0 < ifTrue: [ b neg ->b ]
a 0 < ifTrue: [ b a neg b mod - ->a ]
 
1 dup ->u ->v1
0 dup ->v ->u1
 
while(b) [
b a b /mod ->q ->b ->a
u1 u u1 q * - ->u1 ->u
v1 v v1 q * - ->v1 ->v
]
u v a ;
 
: invmod(a, modulus)
a modulus euclid 1 == ifFalse: [ drop drop null return ]
drop dup 0 < ifTrue: [ modulus + ] ;
Output:
42 2017 invmod println
1969

PARI/GP[edit]

Mod(1/42,2017)

Pascal[edit]

 
// increments e step times until bal is greater than t
// repeats until bal = 1 (mod = 1) and returns count
// bal will not be greater than t + e
 
function modInv(e, t : integer) : integer;
var
d : integer;
bal, count, step : integer;
begin
d := 0;
if e < t then
begin
count := 1;
bal := e;
repeat
step := ((t-bal) DIV e)+1;
bal := bal + step * e;
count := count + step;
bal := bal - t;
until bal = 1;
d := count;
end;
modInv := d;
end;

Testing:

    Writeln(modInv(42,2017));
Output:
1969

Perl[edit]

Various CPAN modules can do this, such as:

use bigint; say 42->bmodinv(2017);
# or
use Math::ModInt qw/mod/; say mod(42, 2017)->inverse->residue;
# or
use Math::Pari qw/PARI lift/; say lift PARI "Mod(1/42,2017)";
# or
use Math::GMP qw/:constant/; say 42->bmodinv(2017);
# or
use ntheory qw/invmod/; say invmod(42, 2017);

or we can write our own:

sub invmod {
my($a,$n) = @_;
my($t,$nt,$r,$nr) = (0, 1, $n, $a % $n);
while ($nr != 0) {
# Use this instead of int($r/$nr) to get exact unsigned integer answers
my $quot = int( ($r - ($r % $nr)) / $nr );
($nt,$t) = ($t-$quot*$nt,$nt);
($nr,$r) = ($r-$quot*$nr,$nr);
}
return if $r > 1;
$t += $n if $t < 0;
$t;
}
 
say invmod(42,2017);

Notes: Special cases to watch out for include (1) where the inverse doesn't exist, such as invmod(14,28474), which should return undef or raise an exception, not return a wrong value. (2) the high bit of a or n is set, e.g. invmod(11,2**63), (3) negative first arguments, e.g. invmod(-11,23). The modules and code above handle these cases, but some other language implementations for this task do not.

Perl 6[edit]

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)

PHP[edit]

Algorithm Implementation

<?php
function invmod($a,$n){
if ($n < 0) $n = -$n;
if ($a < 0) $a = $n - (-$a % $n);
$t = 0; $nt = 1; $r = $n; $nr = $a % $n;
while ($nr != 0) {
$quot= intval($r/$nr);
$tmp = $nt; $nt = $t - $quot*$nt; $t = $tmp;
$tmp = $nr; $nr = $r - $quot*$nr; $r = $tmp;
}
if ($r > 1) return -1;
if ($t < 0) $t += $n;
return $t;
}
printf("%d\n", invmod(42, 2017));
?>
Output:
1969

PicoLisp[edit]

Translation of: C
(de modinv (A B)
(let (B0 B X0 0 X1 1 Q 0 T1 0)
(while (< 1 A)
(setq
Q (/ A B)
T1 B
B (% A B)
A T1
T1 X0
X0 (- X1 (* Q X0))
X1 T1 ) )
(if (lt0 X1) (+ X1 B0) X1) ) )
 
(println
(modinv 42 2017) )
 
(bye)

PL/I[edit]

Translation of: REXX
*process source attributes xref or(!);
/*--------------------------------------------------------------------
* 13.07.2015 Walter Pachl
*-------------------------------------------------------------------*/

minv: Proc Options(main);
Dcl (x,y) Bin Fixed(31);
x=42;
y=2017;
Put Edit('modular inverse of',x,' by ',y,' ---> ',modinv(x,y))
(Skip,3(a,f(4)));
modinv: Proc(a,b) Returns(Bin Fixed(31));
Dcl (a,b,ob,ox,d,t) Bin Fixed(31);
ob=b;
ox=0;
d=1;
 
If b=1 Then;
Else Do;
Do While(a>1);
q=a/b;
r=mod(a,b);
a=b;
b=r;
t=ox;
ox=d-q*ox;
d=t;
End;
End;
If d<0 Then
d=d+ob;
Return(d);
End;
End;
Output:
modular inverse of  42 by 2017 ---> 1969

PowerShell[edit]

function invmod($a,$n){
if ([int]$n -lt 0) {$n = -$n}
if ([int]$a -lt 0) {$a = $n - ((-$a) % $n)}
 
$t = 0
$nt = 1
$r = $n
$nr = $a % $n
while ($nr -ne 0) {
$q = [Math]::truncate($r/$nr)
$tmp = $nt
$nt = $t - $q*$nt
$t = $tmp
$tmp = $nr
$nr = $r - $q*$nr
$r = $tmp
}
if ($r -gt 1) {return -1}
if ($t -lt 0) {$t += $n}
return $t
}
 
invmod 42 2017
Output:
PS> .\INVMOD.PS1
1969
PS> 

Python[edit]

Implementation of this pseudocode with this.

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

Racket[edit]

 
(require math)
(modular-inverse 42 2017)
 
Output:
 
1969
 

REXX[edit]

/*REXX program calculates and displays the  modular inverse  of an integer  X  modulo Y.*/
parse arg x y . /*obtain two integers from the C.L. */
if x=='' | x=="," then x= 42 /*Not specified? Then use the default.*/
if y=='' | y=="," then y= 2017 /* " " " " " " */
say 'modular inverse of ' x " by " y ' ───► ' modInv(x,y)
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
modInv: parse arg a,b 1 ob; z=0 /*B & OB are obtained from the 2nd arg.*/
$=1
if b\=1 then do while a>1
parse value a/b a//b b z with q b a t
z=$ - q*z; $=trunc(t)
end /*while*/
if $<0 then $=$+ob
return $

output   when using the default inputs of:   42   2017

modular inverse of  42  by  2017  ───►  1969

Ring[edit]

 
see "42 %! 2017 = " + multInv(42, 2017) + nl
 
func multInv a,b
b0 = b
x0 = 0
multInv = 1
if b = 1 return 0 ok
while a > 1
q = floor(a / b)
t = b
b = a % b
a = t
t = x0
x0 = multInv - q * x0
multInv = t
end
if multInv < 0 multInv = multInv + b0 ok
return multInv
 

Output:

42 %! 2017 = 1969

Ruby[edit]

#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 'The maths are broken!'
end
x % et
end
?> invmod(42,2017)
=> 1969

Run BASIC[edit]

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
Output:
1969

Scala[edit]

Based on the Handbook of Applied Cryptography, Chapter 2. See http://cacr.uwaterloo.ca/hac/ .

 
def gcdExt(u: Int, v: Int): (Int, Int, Int) = {
@tailrec
def aux(a: Int, b: Int, x: Int, y: Int, x1: Int, x2: Int, y1: Int, y2: Int): (Int, Int, Int) = {
if(b == 0) (x, y, a) else {
val (q, r) = (a / b, a % b)
aux(b, r, x2 - q * x1, y2 - q * y1, x, x1, y, y1)
}
}
aux(u, v, 1, 0, 0, 1, 1, 0)
}
 
def modInv(a: Int, m: Int): Option[Int] = {
val (i, j, g) = gcdExt(a, m)
if (g == 1) Option(if (i < 0) i + m else i) else Option.empty
}

Translated from C++ (on this page)

 
def modInv(a: Int, m: Int, x:Int = 1, y:Int = 0) : Int = if (m == 0) x else modInv(m, a%m, y, x - y*(a/m))
 
Output:
scala> modInv(2,4)
res1: Option[Int] = None

scala> modInv(42, 2017)
res2: Option[Int] = Some(1976)

Seed7[edit]

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.

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;

Original source: [2]

Sidef[edit]

Built-in:

say 42.modinv(2017)

Algorithm implementation:

func invmod(a, n) {
var (t, nt, r, nr) = (0, 1, n, a % n)
while (nr != 0) {
var quot = int((r - (r % nr)) / nr);
(nt, t) = (t - quot*nt, nt);
(nr, r) = (r - quot*nr, nr);
}
r > 1 && return()
t < 0 && (t += n)
t
}
 
say invmod(42, 2017)
Output:
1969

Tcl[edit]

Translation of: Haskell
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
}

Demonstrating

puts "42 %! 2017 = [modInv 42 2017]"
catch {
puts "2 %! 4 = [modInv 2 4]"
} msg; puts $msg
Output:
42 %! 2017 = 1969
no inverse exists of 2 %! 4

tsql[edit]

;WITH Iterate(N,A,B,X0,X1)
AS (
SELECT
1
,CASE WHEN @a < 0 THEN @b-(-@a % @b) ELSE @a END
,CASE WHEN @b < 0 THEN -@b ELSE @b END
,0
,1
UNION ALL
SELECT
N+1
,B
,A%B
,X1-((A/B)*X0)
,X0
FROM Iterate
WHERE A != 1 AND B != 0
),
ModularInverse(Result)
AS (
SELECT
-1
FROM Iterate
WHERE A != 1 AND B = 0
UNION ALL
SELECT
TOP(1)
CASE WHEN X1 < 0 THEN X1+@b ELSE X1 END AS Result
FROM Iterate
WHERE (SELECT COUNT(*) FROM Iterate WHERE A != 1 AND B = 0) = 0
ORDER BY N DESC
)
SELECT *
FROM ModularInverse

uBasic/4tH[edit]

Translation of: Perl
Print FUNC(_mul_inv(42, 2017))
Print FUNC(_mul_inv(40, 1))
Print FUNC(_mul_inv(52, -217))
Print FUNC(_mul_inv(-486, 217))
Print FUNC(_mul_inv(40, 2018))
 
End
 
_mul_inv Param(2)
Local(6)
 
If (b@ < 0) b@ = -b@
If (a@ < 0) a@ = b@ - (-a@ % b@)
c@ = 0 : d@ = 1 : e@ = b@ : f@ = a@ % b@
 
Do Until (f@ = 0)
g@ = e@/f@
h@ = d@ : d@ = c@ - g@*d@ : c@ = h@
h@ = f@ : f@ = e@ - g@*f@ : e@ = h@
Loop
 
If (e@ > 1) Return (-1) ' No inverse'
If (c@ < 0) c@ = c@ + b@
Return (c@)
Output:
1969
0
96
121
-1

0 OK, 0:156

XPL0[edit]

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");
]
Output:
1969

zkl[edit]

Translation of: Haskell
fcn gcdExt(a,b){
if(b==0) return(1,0,a);
q,r:=a.divr(b); s,t,g:=gcdExt(b,r); return(t,s-q*t,g);
}
fcn modInv(a,m){i,_,g:=gcdExt(a,m); if(g==1) {if(i<0)i+m} else Void}

divr(a,b) is [integer] (a/b,remainder)

Output:
modInv(2,4)  //-->Void
modInv(42,2017)  //-->1969