Modular inverse: Difference between revisions

Content added Content deleted

Inline

Latest revision as of 10:25, 4 January 2024

In modular arithmetic, the modular multiplicative inverse of an integer a modulo m is an integer x such that

a\,x\equiv 1{\pmod {m}}.

Or in other words, such that:

\exists k\in \mathbb {Z} ,\qquad a\,x=1+k\,m

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.

Action!

INT FUNC ModInverse(INT a,b)
  INT t,nt,r,nr,q,tmp

  IF b<0 THEN b=-b FI
  IF a<0 THEN a=b-(-a MOD b) FI
  t=0 nt=1
  r=b nr=a MOD b
  WHILE nr#0
  DO
    q=r/nr
    tmp=nt nt=t-q*nt t=tmp
    tmp=nr nr=r-q*nr r=tmp
  OD
  IF r>1 THEN
    RETURN (-1)
  FI
  IF t<0 THEN
    t==+b
  FI
RETURN (t)

PROC Test(INT a,b)
  INT res

  res=ModInverse(a,b)
  IF res>=0 THEN
    PrintF("%I MODINV %I=%I%E",a,b,res)
  ELSE
    PrintF("%I MODINV %I has no result%E",a,b)
  FI
RETURN

PROC Main()
  Test(42,2017)
  Test(40,1)
  Test(52,-217)
  Test(-486,217)
  Test(40,2018)
RETURN

Output:

Screenshot from Atari 8-bit computer

42 MODINV 2017=1969
40 MODINV 1=0
52 MODINV -217=96
-486 MODINV 217=121
40 MODINV 2018 has no result

Ada

with Ada.Text_IO;use Ada.Text_IO;
procedure modular_inverse is
  -- inv_mod calculates the inverse of a mod n. We should have n>0 and, at the end, the contract is a*Result=1 mod n
  -- If this is false then we raise an exception (don't forget the -gnata option when you compile 
  function inv_mod (a : Integer; n : Positive) return Integer with post=> (a * inv_mod'Result) mod n = 1 is 
    -- To calculate the inverse we do as if we would calculate the GCD with the Euclid extended algorithm 
    -- (but we just keep the coefficient on a)
    function inverse (a, b, u, v : Integer) return Integer is
     (if b=0 then u else inverse (b, a mod b, v, u-(v*a)/b));
  begin
    return inverse (a, n, 1, 0);
  end inv_mod;
begin 
  -- This will output -48 (which is correct)
  Put_Line (inv_mod (42,2017)'img);
  -- The further line will raise an exception since the GCD will not be 1
  Put_Line (inv_mod (42,77)'img);
  exception when others => Put_Line ("The inverse doesn't exist."); 
end modular_inverse;

ALGOL 68

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

ATS

Using allocated memory

In addition to solving the task, I demonstrate some aspects of call-by-reference in ATS. In particular, ATS can distinguish at compile time between uninitialized and initialized variables.

The code is written as templates that will expand to code for any of the (non-dependent) signed integer types. In the main program, I use llint (typically exactly the same as long long int in C).

The return value of inverse is a linear optional value. It is allocated in the heap, used once, and freed.

(*

Using the algorithm described at
https://en.wikipedia.org/w/index.php?title=Extended_Euclidean_algorithm&oldid=1135569411#Modular_integers

*)

#include "share/atspre_staload.hats"

fn {tk : tkind}
division_with_nonnegative_remainder
          (n : g0int tk, d : g0int tk,
          (* q and r are called by reference, and start out
             uninitialized. *)
          q : &g0int tk? >> g0int tk,
          r : &g0int tk? >> g0int tk)
    : void =
  let
    (* The C optimizer most likely will reduce these these two
       divisions to just one. They are simply synonyms for C '/' and
       '%', and perform division that rounds the quotient towards
       zero. *)
    val q0 = g0int_div (n, d)
    val r0 = g0int_mod (n, d)
  in
    (* The following calculation results in 'floor division', if the
       divisor is positive, or 'ceiling division', if the divisor is
       negative. This choice of method results in the remainder never
       being negative. *)
    if isgtez n || iseqz r0 then
      (q := q0; r := r0)
    else if isltz d then
      (q := succ q0; r := r0 - d)
    else
      (q := pred q0; r := r0 + d)
  end

fn {tk : tkind}
inverse (a : g0int tk, n : g0int tk) : Option_vt (g0int tk) =
  let
    typedef integer = g0int tk

    fun
    loop (t : integer, newt : integer,
          r : integer, newr : integer) : Option_vt integer =
      if iseqz newr then
        begin
          if r > g0i2i 1 then
            None_vt ()
          else if t < g0i2i 0 then
            Some_vt (t + n)
          else
            Some_vt t
        end
      else
        let
          (* These become C variables. *)
          var quotient : g0int tk?
          var remainder : g0int tk?

          (* Show the type AT COMPILE TIME. *)
          prval _ = $showtype quotient
          prval _ = $showtype remainder

          val () =
            division_with_nonnegative_remainder
              (r, newr, quotient, remainder)

          (* THE TYPES WILL HAVE CHANGED, because the storage is
             initialized by the call to
             division_with_nonnegative_remainder. *)
          prval _ = $showtype quotient
          prval _ = $showtype remainder

          val t = newt
          and newt = t - (quotient * newt)
          and r = newr
          and newr = remainder
        in
          loop (t, newt, r, newr)
        end
  in
    loop (g0i2i 0, g0i2i 1, n, a)
  end

implement
main0 () =
  case+ inverse (42LL, 2017LL) of
  | ~ None_vt () => println! "There is no inverse."
  | ~ Some_vt value => println! value

Output:

First compile the program:

$ patscc -DATS_MEMALLOC_LIBC -g -O2 modular-inverse.dats
**SHOWTYPE[UP]**(/home/trashman/src/chemoelectric/rosettacode-contributions/modular-inverse.dats: 1786(line=62, offs=31) -- 1794(line=62, offs=39)): S2Etop(knd=0; S2Eapp(S2Ecst(g0int_t0ype); S2Evar(tk(8481)))): S2RTbas(S2RTBASimp(1; t@ype))
**SHOWTYPE[UP]**(/home/trashman/src/chemoelectric/rosettacode-contributions/modular-inverse.dats: 1825(line=63, offs=31) -- 1834(line=63, offs=40)): S2Etop(knd=0; S2Eapp(S2Ecst(g0int_t0ype); S2Evar(tk(8481)))): S2RTbas(S2RTBASimp(1; t@ype))
**SHOWTYPE[UP]**(/home/trashman/src/chemoelectric/rosettacode-contributions/modular-inverse.dats: 2137(line=72, offs=31) -- 2145(line=72, offs=39)): S2Eapp(S2Ecst(g0int_t0ype); S2Evar(tk(8481))): S2RTbas(S2RTBASimp(1; t@ype))
**SHOWTYPE[UP]**(/home/trashman/src/chemoelectric/rosettacode-contributions/modular-inverse.dats: 2176(line=73, offs=31) -- 2185(line=73, offs=40)): S2Eapp(S2Ecst(g0int_t0ype); S2Evar(tk(8481))): S2RTbas(S2RTBASimp(1; t@ype))

You may notice there is a subtle change in the type of quotient and remainder, once they have been initialized. ATS is making it safe not to initialize variables (with phony values) when you declare them.

Now run the program:

$ ./a.out
1969

Safely avoiding the need for an allocator

Here I demonstrate an optional value that requires no runtime overhead at all, but which is safe. If there is no inverse, then the compiler knows that inverse_value is still uninitialized, and will not let you use its value. (Try it and see.)

(*
Using the algorithm described at
https://en.wikipedia.org/w/index.php?title=Extended_Euclidean_algorithm&oldid=1135569411#Modular_integers
*)

#include "share/atspre_staload.hats"

fn {tk : tkind}
division_with_nonnegative_remainder
          (n : g0int tk, d : g0int tk,
          (* q and r are called by reference, and start out
             uninitialized. *)
          q : &g0int tk? >> g0int tk,
          r : &g0int tk? >> g0int tk)
    : void =
  let
    (* The C optimizer most likely will reduce these these two
       divisions to just one. They are simply synonyms for C '/' and
       '%', and perform division that rounds the quotient towards
       zero. *)
    val q0 = g0int_div (n, d)
    val r0 = g0int_mod (n, d)
  in
    (* The following calculation results in 'floor division', if the
       divisor is positive, or 'ceiling division', if the divisor is
       negative. This choice of method results in the remainder never
       being negative. *)
    if isgtez n || iseqz r0 then
      (q := q0; r := r0)
    else if isltz d then
      (q := succ q0; r := r0 - d)
    else
      (q := pred q0; r := r0 + d)
  end

fn {tk : tkind}
inverse (a : g0int tk, n : g0int tk,
         inverse_exists : &bool? >> bool exists,
         inverse_value  : &g0int tk? >> opt (g0int tk, exists))
    : #[exists: bool] void =
  let
    typedef integer = g0int tk

    fun
    loop (t : integer, newt : integer,
          r : integer, newr : integer,
          inverse_exists : &bool? >> bool exists,
          inverse_value  : &g0int tk? >> opt (g0int tk, exists))
        : #[exists: bool] void =
      if iseqz newr then
        begin
          if r > g0i2i 1 then
            let
              val () = inverse_exists := false
              prval () = opt_none inverse_value
            in
            end
          else if t < g0i2i 0 then
            let
              val () = inverse_exists := true
              val () = inverse_value := t + n
              prval () = opt_some inverse_value
            in
            end
          else
            let
              val () = inverse_exists := true
              val () = inverse_value := t
              prval () = opt_some inverse_value
            in
            end
        end
      else
        let
          (* These become C variables. *)
          var quotient : g0int tk?
          var remainder : g0int tk?

          val () =
            division_with_nonnegative_remainder
              (r, newr, quotient, remainder)

          val t = newt
          and newt = t - (quotient * newt)
          and r = newr
          and newr = remainder
        in
          loop (t, newt, r, newr, inverse_exists, inverse_value)
        end
  in
    loop (g0i2i 0, g0i2i 1, n, a, inverse_exists, inverse_value)
  end

implement
main0 () =
  let
    var inverse_exists : bool?
    var inverse_value  : llint?
  in
    inverse (42LL, 2017LL, inverse_exists, inverse_value);
    if inverse_exists then
      let
        prval () = opt_unsome inverse_value
      in
        println! inverse_value
      end
    else
      let
        prval () = opt_unnone inverse_value
      in
        println! "There is no inverse."
      end
  end

Output:

There is no need to tell patscc what allocator to use, because none is used.

$ patscc -g -O2 modular-inverse-noheap.dats && ./a.out
1969

AutoHotkey

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:

AWK

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

BASIC

ASIC

Translation of: Pascal

REM Modular inverse
E = 42
T = 2017
GOSUB CalcModInv:
PRINT ModInv
END

CalcModInv:
REM Increments E Step times until Bal is greater than T 
REM Repeats until Bal = 1 (MOD = 1) and returns Count
REM Bal will not be greater than T + E
D = 0
IF E < T THEN
  Bal = E
  Count = 1
  Loop:
    Step = T - Bal
    Step = Step / E
    Step = Step + 1
    REM So ... Step = (T - Bal) / E + 1
    StepTimesE = Step * E
    Bal = Bal + StepTimesE
    Count = Count + Step
    Bal = Bal - T
    IF Bal <> 1 THEN Loop:
  D = Count
ENDIF
ModInv = D
RETURN

Output:

BASIC256

print multInv(42, 2017)
end

function multInv(a,b)
	x0 = 0
	b0	= b
	multInv	= 1
	if b = 1 then return
	while a > 1
		q	= a / b
		t	= b 
		b	= a mod b
		a	= t
		t	= x0 
		x0	= multInv - q * x0
		multInv	= int(t)
	end while
	if multInv < 0 then return multInv + b0
end function

Output:

Chipmunk Basic

Works with: Chipmunk Basic version 3.6.4

Works with: QBasic

Translation of: BASIC256

10 CLS
20 CALL modularinverse(42, 2017)
30 CALL modularinverse(40, 1)
40 END
50 SUB modularinverse(e,t)
60 d = 0
70 IF e < t THEN
80  b = e
90  c = 1
100  WHILE b > 1
110   s = INT(((t-b)/e)+1)
120   b = b+s*e
130   c = c+s
140   b = b-t
150  WEND
160 d = c
170 ENDIF
180 m = d
190 PRINT m
200 END SUB

Minimal BASIC

Translation of: Pascal

Works with: Applesoft BASIC

Works with: Commodore BASIC version 3.5

Works with: Nascom ROM BASIC version 4.7

10 REM Modular inverse
20 LET E = 42
30 LET T = 2017
40 GOSUB 500
50 PRINT M
60 END
490 REM Calculate modular inverse
500 LET D = 0
510 IF E >= T THEN 600
520 LET B = E
530 LET C = 1
540 LET S1 = INT((T-B)/E)+1
550 LET B = B+S1*E
560 LET C = C+S1
570 LET B = B-T
580 IF B <> 1 THEN 540
590 LET D = C
600 LET M = D
610 RETURN

QBasic

The Chipmunk Basic solution works without any changes.

True BASIC

SUB modularinverse(e,t)
    LET d = 0
    IF e < t then
       LET b = e
       LET c = 1
       DO WHILE b > 1
          LET s = int(((t-b)/e)+1)
          LET b = b+s*e
          LET c = c+s
          LET b = b-t
       LOOP
       LET d = c
    END IF
    LET m = d
    PRINT m
END SUB

CALL modularinverse(42,2017)
CALL modularinverse(40,1)
END

Batch File

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:

BCPL

get "libhdr"

let mulinv(a, b) =
    b<0 -> mulinv(a, -b),
    a<0 -> mulinv(b - (-a rem b), b),
    valof
    $(  let t, nt, r, nr = 0, 1, b, a rem b
        until nr = 0
        $(  let tmp, q = ?, r / nr
            tmp := nt ; nt := t - q*nt ; t := tmp
            tmp := nr ; nr := r - q*nr ; r := tmp
        $)
        resultis r>1 -> -1, 
                 t<0 -> t + b, 
                 t
    $)
    
let show(a, b) be
$(  let mi = mulinv(a, b)
    test mi>=0 
        do writef("%N, %N -> %N*N", a, b, mi)
        or writef("%N, %N -> no inverse*N", a, b)
$)

let start() be
$(  show(42, 2017)
    show(40, 1)
    show(52, -217)
    show(-486, 217)
    show(40, 2018)
$)

Output:

42, 2017 -> 1969
40, 1 -> 0
52, -217 -> 96
-486, 217 -> 121
40, 2018 -> no inverse

Bracmat

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

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

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:

C#

public class Program
{
    static void Main()
    {
        System.Console.WriteLine(42.ModInverse(2017));
    }
}

public static class IntExtensions
{
    public static int ModInverse(this int a, int m)
    {
        if (m == 1) return 0;
        int m0 = m;
        (int x, int y) = (1, 0);

        while (a > 1) {
            int q = a / m;
            (a, m) = (m, a % m);
            (x, y) = (y, x - q * y);
        }
        return x < 0 ? x + m0 : x;
    }
}

C++

Iterative implementation

Translation of: C

#include <iostream>
 
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) {
	std::cout << mul_inv(42, 2017) << std::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;
}

Clojure

(ns test-p.core
  (:require [clojure.math.numeric-tower :as math]))

(defn extended-gcd
  "The extended Euclidean algorithm--using Clojure code from RosettaCode for Extended Eucliean
  (see http://en.wikipedia.orwiki/Extended_Euclidean_algorithm)
  Returns a list containing the GCD and the Bézout coefficients
  corresponding to the inputs with the result: gcd followed by bezout coefficients "
  [a b]
  (cond (zero? a) [(math/abs b) 0 1]
        (zero? b) [(math/abs a) 1 0]
        :else (loop [s 0
                     s0 1
                     t 1
                     t0 0
                     r (math/abs b)
                     r0 (math/abs a)]
                (if (zero? r)
                  [r0 s0 t0]
                  (let [q (quot r0 r)]
                    (recur (- s0 (* q s)) s
                           (- t0 (* q t)) t
                           (- r0 (* q r)) r))))))

(defn mul_inv
  " Get inverse using extended gcd.  Extended GCD returns
    gcd followed by bezout coefficients. We want the 1st coefficients
   (i.e. second of extend-gcd result).  We compute mod base so result
    is between 0..(base-1) "
  [a b]
  (let [b (if (neg? b) (- b) b)
        a (if (neg? a) (- b (mod (- a) b)) a)
        egcd (extended-gcd a b)]
      (if (= (first egcd) 1)
        (mod (second egcd) b)
        (str "No inverse since gcd is: " (first egcd)))))


(println (mul_inv 42 2017))
(println (mul_inv 40 1))
(println (mul_inv 52 -217))
(println (mul_inv -486 217))
(println (mul_inv 40 2018))

Output:

1969
0
96
121
No inverse since gcd is: 2

CLU

Translation of: Perl

mul_inv = proc (a, b: int) returns (int) signals (no_inverse)
    if b<0 then b := -b end
    if a<0 then a := b - (-a // b) end
    t: int := 0
    nt: int := 1
    r: int := b
    nr: int := a // b
    while nr ~= 0 do
        q: int := r / nr
        t, nt := nt, t - q*nt
        r, nr := nr, r - q*nr
    end
    if r>1 then signal no_inverse end
    if t<0 then t := t+b end
    return(t)
end mul_inv

start_up = proc ()
    pair = struct[a, b: int]
    tests: sequence[pair] := sequence[pair]$
       [pair${a: 42, b: 2017},
        pair${a: 40, b: 1},
        pair${a: 52, b: -217},
        pair${a: -486, b: 217},
        pair${a: 40, b: 2018}]
        
    po: stream := stream$primary_output()
    for test: pair in sequence[pair]$elements(tests) do
        stream$puts(po, int$unparse(test.a) || ", "
                     || int$unparse(test.b) || " -> ")
        stream$putl(po, int$unparse(mul_inv(test.a, test.b)))
        except when no_inverse:
            stream$putl(po, "no modular inverse")
        end 
    end
end start_up

Output:

42, 2017 -> 1969
40, 1 -> 0
52, -217 -> 96
-486, 217 -> 121
40, 2018 -> no modular inverse

Comal

0010 FUNC mulinv#(a#,b#) CLOSED
0020   IF b#<0 THEN b#:=-b#
0030   IF a#<0 THEN a#:=b#-(-a# MOD b#)
0040   t#:=0;nt#:=1;r#:=b#;nr#:=a# MOD b#
0050   WHILE nr#<>0 DO
0060     q#:=r# DIV nr#
0070     tmp#:=nt#;nt#:=t#-q#*nt#;t#:=tmp#
0080     tmp#:=nr#;nr#:=r#-q#*nr#;r#:=tmp#
0090   ENDWHILE
0100   IF r#>1 THEN RETURN -1
0110   IF t#<0 THEN t#:+b#
0120   RETURN t#
0130 ENDFUNC mulinv#
0140 //
0150 WHILE NOT EOD DO
0160   READ a#,b#
0170   PRINT a#,", ",b#," -> ",mulinv#(a#,b#)
0180 ENDWHILE
0190 END
0200 //
0210 DATA 42,2017,40,1,52,-217,-486,217,40,2018

Output:

42, 2017 -> 1969
40, 1 -> 0
52, -217 -> 96
-486, 217 -> 121
40, 2018 -> -1

Common Lisp

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

Cowgol

include "cowgol.coh";

sub mulinv(a: int32, b: int32): (t: int32) is
    if b<0 then b := -b; end if;
    if a<0 then a := b - (-a % b); end if;
    t := 0;
    var nt: int32 := 1;
    var r := b;
    var nr := a % b;
    while nr != 0 loop
        var q := r / nr;
        var tmp := nt; nt := t - q*nt; t := tmp;
        tmp := nr; nr := r - q*nr; r := tmp;
    end loop;
    if r>1 then t := -1;
    elseif t<0 then t := t + b;
    end if;
end sub;

record Pair is
    a: int32;
    b: int32;
end record;

var data: Pair[] := {
    {42, 2017}, 
    {40, 1}, 
    {52, -217}, 
    {-486, 217}, 
    {40, 2018}
};

var i: @indexof data := 0;
while i < @sizeof data loop
    print_i32(data[i].a as uint32);
    print(", ");
    print_i32(data[i].b as uint32);
    print(" -> ");
    var mi := mulinv(data[i].a, data[i].b);
    if mi<0
        then print("no inverse");
        else print_i32(mi as uint32);
    end if;
    print_nl();
    i := i + 1;
end loop;

Output:

42, 2017 -> 1969
40, 1 -> 0
52, 4294967079 -> 96
4294966810, 217 -> 121
40, 2018 -> no inverse

Craft Basic

let e = 42
let t = 2017

gosub modularinverse

end

sub modularinverse

	let d = 0

	if e < t then

		let b = e
		let c = 1

		do

			let s = int(((t - b) / e) + 1)
			let b = b + s * e
			let c = c + s
			let b = b - t

		loop b <> 1

		let d = c

	endif

	let m = d

	print m

return

Output:

Crystal

Translation of: Ruby

def modinv(a0, m0)
  return 1 if m0 == 1
  a, m = a0, m0
  x0, inv = 0, 1
  while a > 1
    inv -= (a // m) * x0
    a, m = m, a % m
    x0, inv = inv, x0
  end
  inv += m0 if inv < 0
  inv
end

Output:

> modinv(42,2017)
=> 1969

D

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:

dc

Translation of: C

This solution prints the inverse u only if it exists (a*u = 1 mod m).

dc -e "[m=]P?dsm[a=]P?dsa1sv[dsb~rsqlbrldlqlv*-lvsdsvd0<x]dsxxldd[dlmr+]sx0>xdla*lm%[p]sx1=x"

If ~ is not implemented, it can be replaced by SdSnlnld/LnLd%.

Replace [p]sx1=x at the end by [pq]sx1=x16i6E6F7420636F7072696D65P if an error message "not coprime" is desired.

Output:

m=2 800^1+
a=37
342411551958695219479776173037037562556082184118925013641969995739234\
344644689214483533004909620355470582887300743869103978073598454778206\
829469635119691272637318902731800747596752668736012071540136041369140\
1228044652005748974399041408477572

m=2017
a=42
1969

m=42
a=7

Delphi

See #Pascal.

Draco

proc mulinv(int a, b) int:
    int t, nt, r, nr, q, tmp;
    if b<0 then b := -b fi;
    if a<0 then a := b - (-a % b) fi;
    t := 0; nt := 1; r := b; nr := a % b;
    while nr /= 0 do
        q := r / nr;
        tmp := nt; nt := t - q*nt; t := tmp;
        tmp := nr; nr := r - q*nr; r := tmp
    od;
    if r>1 then -1
    elif t<0 then t+b
    else t
    fi
corp

proc show(int a, b) void:
    int mi;
    mi := mulinv(a, b);
    if mi>=0 
        then writeln(a:5, ", ", b:5, " -> ", mi:5)
        else writeln(a:5, ", ", b:5, " -> no inverse")
    fi
corp

proc main() void:
    show(42, 2017);
    show(40, 1);
    show(52, -217);
    show(-486, 217);
    show(40, 2018)
corp

Output:

   42,  2017 ->  1969
   40,     1 ->     0
   52,  -217 ->    96
 -486,   217 ->   121
   40,  2018 -> no inverse

EasyLang

Translation of: AWK

func mod_inv a b .
   b0 = b
   x1 = 1
   if b = 1
      return 1
   .
   while a > 1
      q = a div b
      t = b
      b = a mod b
      a = t
      t = x0
      x0 = x1 - q * x0
      x1 = t
   .
   if x1 < 0
      x1 += b0
   .
   return x1
.
print mod_inv 42 2017

EchoLisp

(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

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:

ERRE

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:

F#

// Calculate the inverse of a (mod m)
// See here for eea specs: 
// https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
let modInv m a = 
  let rec eea t t' r r' =
    match r' with
    | 0 -> t
    | _ -> 
      let div = r/r'
      eea t' (t - div * t') r' (r - div * r')
  (m + eea 0 1 m a) % m

Input:

// Inverse of 347 (mod 29)
modInv 29 347

Output:

Factor

USE: math.functions
42 2017 mod-inv

Output:

Forth

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 ;

ANS Forth version without locals

: modinv ( a m - inv)
  dup 1-              \ a m (m != 1)?
  if                  \ a m
    tuck 1 0          \ m0 a m 1 0
    begin             \ m0 a m inv x0
      2>r over 1 >    \ m0 a m (a > 1)?       R: inv x0
    while             \ m0 a m                R: inv x0
      tuck /mod       \ m0 m (a mod m) (a/m)  R: inv x0
      r> tuck *       \ m0 a' m' x0 (a/m)*x0  R: inv
      r> swap -       \ m0 a' m' x0 (inv-q)   R:
    repeat            \ m0 a' m' inv' x0'
    2drop             \ m0                    R: inv x0
    2r> drop          \ m0 inv                R:
    dup 0<            \ m0 inv (inv < 0)?
    if over + then    \ m0 (inv + m0)
  then                \ x inv'
  nip                 \ inv
;

42 2017 invmod . 1969
42 2017 modinv . 1969

Fortran

program modular_inverse_task

  implicit none

  write (*,*) inverse (42, 2017)

contains

  ! Returns -1 if there is no inverse. I assume n > 0. The algorithm
  ! is described at
  ! https://en.wikipedia.org/w/index.php?title=Extended_Euclidean_algorithm&oldid=1135569411#Modular_integers
  function inverse (a, n) result (inverse_value)
    integer, intent(in) :: a, n
    integer :: inverse_value

    integer :: t, newt
    integer :: r, newr
    integer :: quotient, remainder, tmp

    if (n <= 0) error stop
    t = 0;  newt = 1
    r = n;  newr = a
    do while (newr /= 0)
       remainder = modulo (r, newr) ! Floor division.
       quotient = (r - remainder) / newr
       tmp = newt;  newt = t - (quotient * newt);  t = tmp
       r   = newr;  newr = remainder
    end do
    if (r > 1) then
       inverse_value = -1
    else if (t < 0) then
       inverse_value = t + n
    else
       inverse_value = t
    end if
  end function inverse

end program modular_inverse_task

Output:

$ gfortran -Wall -Wextra modular_inverse_task.f90 && ./a.out
        1969

FreeBASIC

' version 10-07-2018
' compile with: fbc -s console

Type ext_euclid
    Dim As Integer a, b
End Type

' "Table method" aka "The Magic Box"
Function magic_box(x As Integer, y As Integer) As ext_euclid

    Dim As Integer a(1 To 128), b(1 To 128), d(1 To 128), k(1 To 128)

    a(1) = 1 : b(1) = 0 : d(1) = x
    a(2) = 0 : b(2) = 1 : d(2) = y : k(2) = x \ y

    Dim As Integer i = 2

    While Abs(d(i)) <> 1
        i += 1
        a(i) = a(i -2) - k(i -1) * a(i -1)
        b(i) = b(i -2) - k(i -1) * b(i -1)
        d(i) = d(i -2) Mod d(i -1)
        k(i) = d(i -1) \ d(i)
        'Print a(i),b(i),d(i),k(i)
        If d(i -1) Mod d(i) = 0 Then Exit While
    Wend
    
    If d(i) = -1 Then '  -1 * (ab + by) = -1 * -1 ==> -ab -by = 1 
        a(i) = -a(i)
        b(i) = -b(i)
    End If

    Function = Type( a(i), b(i) )

End Function
' ------=< MAIN >=------

Dim As Integer x, y, gcd
Dim As ext_euclid result

Do
    Read x, y
    If x = 0 AndAlso y = 0 Then Exit Do
    result = magic_box(x, y)
    With result
        gcd = .a * x + .b * y
        Print "a * "; Str(x); " + b * "; Str(y);  
        Print " = GCD("; Str(x); ", "; Str(y); ") ="; gcd
        If gcd > 1 Then 
            Print "No solution, numbers are not coprime"
        Else 
            Print "a = "; .a; ", b = ";.b
            Print "The Modular inverse of "; x; " modulo "; y; " = ";
            While .a < 0 : .a += IIf(y > 0, y, -y) : Wend 
            Print .a
            'Print "The Modular inverse of "; y; " modulo "; x; " = ";
            'While .b < 0 : .b += IIf(x > 0, x, -x) : Wend 
            'Print .b
        End if
    End With
    Print
Loop

Data 42, 2017
Data 40, 1
Data 52, -217
Data -486, 217
Data 40, 2018
Data 0, 0

' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End

Output:

a * 42 + b * 2017 = GCD(42, 2017) = 1
a = -48, b =  1
The Modular inverse of  42 modulo  2017 =  1969

a * 40 + b * 1 = GCD(40, 1) = 1
a =  0, b =  1
The Modular inverse of  40 modulo  1 =  0

a * 52 + b * -217 = GCD(52, -217) = 1
a =  96, b =  23
The Modular inverse of  52 modulo -217 =  96

a * -486 + b * 217 = GCD(-486, 217) = 1
a = -96, b = -215
The Modular inverse of -486 modulo  217 =  121

a * 40 + b * 2018 = GCD(40, 2018) = 2
No solution, numbers are not coprime

Frink

println[modInverse[42, 2017]]

Output:

FunL

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:

Go

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:

GW-BASIC

Translation of: Pascal

Works with: PC-BASIC version any

10 ' Modular inverse
20 LET E% = 42
30 LET T% = 2017
40 GOSUB 1000
50 PRINT MODINV%
60 END

990  ' increments e stp (step) times until bal is greater than t
992  ' repeats until bal = 1 (mod = 1) and returns count
994  ' bal will not be greater than t + e 
1000 LET D% = 0
1010 IF E% >= T% THEN GOTO 1140
1020  LET BAL% = E%
1025  ' At least one iteration is necessary
1030  LET STP% = ((T% - BAL%) \ E%) + 1
1040  LET BAL% = BAL% + STP% * E%
1050  LET COUNT% = 1 + STP%
1060  LET BAL% = BAL% - T%
1070  WHILE BAL% <> 1
1080   LET STP% = ((T% - BAL%) \ E%) + 1
1090   LET BAL% = BAL% + STP% * E%
1100   LET COUNT% = COUNT% + STP%
1110   LET BAL% = BAL% - T%
1120  WEND
1130  LET D% = COUNT%
1140 LET MODINV% = D%
1150 RETURN

Output:

Haskell

-- Given a and m, return Just x such that ax = 1 mod m.
-- If there is no such x return Nothing.
modInv :: Int -> Int -> Maybe Int
modInv a m
  | 1 == g = Just (mkPos i)
  | otherwise = Nothing
  where
    (i, _, g) = gcdExt a m
    mkPos x
      | x < 0 = x + m
      | otherwise = x

-- 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 :: Int -> Int -> (Int, Int, Int)
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)

main :: IO ()
main = mapM_ print [2 `modInv` 4, 42 `modInv` 2017]

Output:

Nothing
Just 1969

Icon and Unicon

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

IS-BASIC

100 PRINT MODINV(42,2017)
120 DEF MODINV(A,B)
130   LET B=ABS(B)
140   IF A<0 THEN LET A=B-MOD(-A,B)
150   LET T=0:LET NT=1:LET R=B:LET NR=MOD(A,B)
160   DO WHILE NR<>0
170     LET Q=INT(R/NR)
180     LET TMP=NT:LET NT=T-Q*NT:LET T=TMP
190     LET TMP=NR:LET NR=R-Q*NR:LET R=TMP
200   LOOP
210   IF R>1 THEN
220     LET MODINV=-1
230   ELSE IF T<0 THEN
240     LET MODINV=T+B
250   ELSE
260     LET MODINV=T
270   END IF
280 END DEF

J

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

The BigInteger library has a method for this:

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

Output:

Alternatively, working from first principles.

public final class ModularInverse {

	public static void main(String[] aArgs) {
		System.out.println(inverseModulus(42, 2017));
	}
	
	private static Egcd extendedGCD(int aOne, int aTwo) {
		if ( aOne == 0 ) {
	        return new Egcd(aTwo, 0, 1);
		}
	    Egcd value = extendedGCD(aTwo % aOne, aOne);	    
	    return new Egcd(value.aG, value.aY - ( aTwo / aOne ) * value.aX, value.aX);	    
	}
	    		
	private static int inverseModulus(int aNumber, int aModulus) {
	    Egcd value = extendedGCD(aNumber, aModulus); 
	    return ( value.aG == 1 ) ? ( value.aX + aModulus ) % aModulus : 0;
	}
	
	private static record Egcd(int aG, int aX, int aY) {}

}

Output:

JavaScript

Using brute force.

var modInverse = function(a, b) {
    a %= b;
    for (var x = 1; x < b; x++) {
        if ((a*x)%b == 1) {
            return x;
        }
    }
}

jq

Works with: jq

Works with gojq, the Go implementation of jq

# Integer division:
# If $j is 0, then an error condition is raised;
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
def idivide($j):
  . as $i
  | ($i % $j) as $mod
  | ($i - $mod) / $j ;

# the multiplicative inverse of . modulo $n
def modInv($n):
  if $n == 1 then 1
  else . as $this
  | { r   : $n,
      t   : 0,
      newR: length, # abs
      newT: 1}
  | until(.newR == 0;
        .newR as $newR
        | (.r | idivide($newR)) as $q
        | {r   : $newR,
           t   : .newT,
           newT: (.t - $q * .newT),
           newR: (.r - $q * $newR) } )
	     
  | if (.r|length) != 1 then "\($this) and \($n) are not co-prime." | error
    else .t
    | if .     < 0 then . + $n
      elif $this < 0 then - .
      else .
      end
    end
  end ;

# Example:
42 | modInv(2017)

Output:

Julia

Works with: Julia version 1.2

Built-in

Julia includes a built-in function for this:

invmod(a, b)

C translation

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(a::T, b::T) where T <: Integer
    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

Kotlin

// version 1.0.6

import java.math.BigInteger

fun main(args: Array<String>) {
    val a = BigInteger.valueOf(42)
    val m = BigInteger.valueOf(2017)
    println(a.modInverse(m))
}

Output:

Lambdatalk

Translation of: Phix

 

{def mulinv
 {def mulinv.loop
  {lambda {:t :nt :r :nr}
   {if {not {= :nr 0}}
    then {mulinv.loop :nt 
                      {- :t {* {floor {/ :r :nr}} :nt}}
                      :nr 
                      {- :r {* {floor {/ :r :nr}} :nr}} }
    else {cons :t :r} }}}
 {lambda {:a :n}
  {let { {:a :a} {:n :n}
         {:cons {mulinv.loop 0 
                             1 
                             {if {< :n 0} then {- :n} else :n} 
                             {if {< :a 0} then {- :n {% {- :a} :n}} else :a}}}
       } {if {> {cdr :cons} 1}
          then not invertible
          else {if {< {car :cons} 0}
          then {+ {car :cons} :n}
          else {car :cons} }}}}}
-> mulinv

{mulinv 42 2017}
-> 1969
{mulinv 40 1}
-> 0
{mulinv 52 -217}
-> 96
{mulinv -486 217}
-> 121
{mulinv 40 218}
-> not invertible

m4

Translation of: Mercury

Note that $0 is the name of the macro being evaluated. Therefore, in the following, _$0 is the name of another macro, the same as the name of the first macro, except for an underscore prepended. This is a common idiom.

The core of the program is __inverse recursively calling itself.

m4 is a macro-preprocessor, and so some of the following is there simply to keep things from being echoed to the output. :)

divert(-1)

# I assume non-negative arguments. The algorithm is described at
# https://en.wikipedia.org/w/index.php?title=Extended_Euclidean_algorithm&oldid=1135569411#Modular_integers

define(`inverse',`_$0(eval(`$1'), eval(`$2'))')
define(`_inverse',`_$0($2, 0, 1, $2, $1)')
define(`__inverse',
`dnl  n = $1,  t = $2,  newt = $3,  r = $4,  newr = $5
ifelse(eval($5 != 0), 1, `$0($1, $3,
                             eval($2 - (($4 / $5) * $3)),
                             $5,eval($4 % $5))',
       eval($4 > 1), 1, `no inverse',
       eval($2 < 0), 1, eval($2 + $1),
       $2)')

divert`'dnl
inverse(42, 2017)

Output:

$ m4 modular-inverse-task.m4
1969

MAD

            NORMAL MODE IS INTEGER
            INTERNAL FUNCTION(AA, BB)
            ENTRY TO MULINV.
            A = AA
            B = BB 
            WHENEVER B.L.0, B = -B
            WHENEVER A.L.0, A = B - (-(A-A/B*B))
            T = 0
            NT = 1
            R = B
            NR = A-A/B*B
LOOP        WHENEVER NR.NE.0
                Q = R/NR
                TMP = NT
                NT = T - Q*NT
                T = TMP
                TMP = NR
                NR = R - Q*NR
                R = TMP
                TRANSFER TO LOOP
            END OF CONDITIONAL
            WHENEVER R.G.1, FUNCTION RETURN -1
            WHENEVER T.L.0, T = T+B
            FUNCTION RETURN T
            END OF FUNCTION
            
            INTERNAL FUNCTION(AA, BB)
            VECTOR VALUES FMT = $I5,2H, ,I5,2H: ,I5*$
            ENTRY TO SHOW.
            PRINT FORMAT FMT, AA, BB, MULINV.(AA, BB)
            END OF FUNCTION
            
            SHOW.(42,2017)
            SHOW.(40,1)
            SHOW.(52,-217)
            SHOW.(-486,217)
            SHOW.(40,2018)
            END OF PROGRAM

Output:

   42,  2017:  1969
   40,     1:     0
   52,  -217:    96
 -486,   217:   121
   40,  2018:    -1

Maple

1/42 mod 2017;

Output:

Mathematica/Wolfram Language

Use the built-in function ModularInverse:

ModularInverse[a, m]

For example:

ModularInverse[42, 2017]
1969

Maxima

Using built-in function inv_mod

inv_mod(42,2017);

Output:

Mercury

Works with: Mercury version 22.01.1

%%% -*- mode: mercury; prolog-indent-width: 2; -*-
%%%
%%% Compile with:
%%%    mmc --make --use-subdirs modular_inverse_task
%%%

:- module modular_inverse_task.
:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.

:- implementation.
:- import_module exception.
:- import_module int.

%% inverse(A, N, Inverse). I assume N > 0, and throw an exception if
%% it is not. The predicate fails if there is no inverse (and thus is
%% "semidet"). The algorithm is described at
%% https://en.wikipedia.org/w/index.php?title=Extended_Euclidean_algorithm&oldid=1135569411#Modular_integers
:- pred inverse(int::in, int::in, int::out) is semidet.
inverse(A, N, Inverse) :-
  if (N =< 0) then throw(domain_error("inverse"))
  else inverse_(N, 0, 1, N, A, Inverse).

:- pred inverse_(int::in, int::in, int::in, int::in, int::in,
                 int::out) is semidet.
inverse_(N, T, NewT, R, NewR, Inverse) :-
  if (NewR \= 0)
  then (Quotient = div(R, NewR), % Floor division.
        inverse_(N,
                 NewT, T - (Quotient * NewT),
                 NewR, R - (Quotient * NewR),
                 Inverse))      % Tail recursion.
  else (R =< 1,                 % R =< 1 FAILS if R > 1.
        (if (T < 0) then Inverse = T + N else Inverse = T)).

main(!IO) :-
  if inverse(42, 2017, Inverse)
  then (print(Inverse, !IO), nl(!IO))
  else (print("There is no inverse.", !IO), nl(!IO)).

:- end_module modular_inverse_task.

Output:

$ mmc --make --use-subdirs modular_inverse_task && ./modular_inverse_task
Making Mercury/int3s/modular_inverse_task.int3
Making Mercury/ints/modular_inverse_task.int
Making Mercury/cs/modular_inverse_task.c
Making Mercury/os/modular_inverse_task.o
Making modular_inverse_task
1969

МК-61/52

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

Modula-2

Translation of: C

MODULE ModularInverse;
  FROM InOut IMPORT WriteString, WriteInt, WriteLn;

  TYPE Data = RECORD x : INTEGER;
                     y : INTEGER
              END;

  VAR c  : INTEGER;
      ab : ARRAY [1..5] OF Data;

PROCEDURE mi(VAR a, b : INTEGER): INTEGER;
  VAR t, nt, r, nr, q, tmp : INTEGER;

BEGIN
  b := ABS(b);
  IF a < 0 THEN a := b - (-a MOD b) END;
  t := 0; nt := 1; r := b; nr := a MOD b;
  WHILE (nr # 0) DO
    q := r / nr;
    tmp := nt; nt := t - q * nt; t := tmp;
    tmp := nr; nr := r - q * nr; r := tmp;
  END;
  IF (r > 1) THEN RETURN -1 END;
  IF (t < 0) THEN RETURN t + b END;
  RETURN t;
END mi;

BEGIN
  ab[1].x := 42;   ab[1].y := 2017;
  ab[2].x := 40;   ab[2].y := 1;
  ab[3].x := 52;   ab[3].y := -217;
  ab[4].x := -486; ab[4].y := 217;
  ab[5].x := 40;   ab[5].y := 2018;
  WriteLn;
  WriteString("Modular inverse");
  WriteLn;
  FOR c := 1 TO 5 DO
    WriteInt(ab[c].x, 6); WriteString(", ");
    WriteInt(ab[c].y, 6); WriteString(" = ");
    WriteInt(mi(ab[c].x, ab[c].y),6);
    WriteLn;
  END;
END ModularInverse.

Output:

Modular inverse
    42,   2017 =   1969
    40,      1 =      0
    52,   -217 =     96
  -486,    217 =    121
    40,   2018 =     -1

newLISP

(define (modular-multiplicative-inverse a n)
    (if (< n 0)
        (setf n (abs n)))
        
    (if (< a 0)
        (setf a (- n (% (- 0 a) n))))
    
    (setf t 0)
    (setf nt 1)
    (setf r n)
    (setf nr (mod a n))
    
    (while (not (zero? nr))
        (setf q (int (div r nr)))
        (setf tmp nt)
        (setf nt (sub t (mul q nt)))
        (setf t tmp)
        (setf tmp nr)
        (setf nr (sub r (mul q nr)))
        (setf r tmp))
    
    (if (> r 1)
        (setf retvalue nil))
    
    (if (< t 0)
        (setf retvalue (add t n))
        (setf retvalue t))
    retvalue)  

(println (modular-multiplicative-inverse 42 2017))

Output:

Nim

Translation of: C

proc modInv(a0, b0: int): int =
  var (a, b, x0) = (a0, b0, 0)
  result = 1
  if b == 1: return
  while a > 1:
    result = result - (a div b) * x0
    a = a mod b
    swap a, b
    swap x0, result
  if result < 0: result += b0

echo modInv(42, 2017)

Output:

Oberon-2

Works with: Oxford Oberon-2 Compiler version 3.3.0

(*-*- mode: indented-text; tab-width: 2; -*-*)

MODULE modularInverseInOberon2;

IMPORT Out;

(* Division with a non-negative remainder. This will work no matter
   how your compiler handles DIV (and mine seems not to do what the
   Oberon-2 specification says). *)
PROCEDURE euclidDiv (x, y : INTEGER) : INTEGER;
  VAR q : INTEGER;
BEGIN
  IF 0 <= y THEN      (* Do floor division. *)
    IF 0 <= x THEN
      q := x DIV y
    ELSE
      q := -((-x) DIV y);
      IF (-x) MOD y # 0 THEN q := q - 1 END
    END;
  ELSE                (* Do ceiling division. *)
    IF 0 <= x THEN
      q := -(x DIV (-y))
    ELSE
      q := ((-x) DIV (-y));
      IF (-x) MOD (-y) # 0 THEN q := q + 1 END
    END
  END;
  RETURN q
END euclidDiv;

(* I have added this unit test because, earlier, I posted a buggy
   version of euclidDiv. *)
PROCEDURE testEuclidDiv;
  VAR x, y, q, r : INTEGER;
BEGIN
  FOR x := -100 TO 100 DO
    FOR y := -100 TO 100 DO
      IF y # 0 THEN
        q := euclidDiv (x, y);
        r := x - (q * y);
        IF (r < 0) OR (ABS (y) <= r) THEN
          (* A remainder was outside the expected range. *)
          Out.String ("euclidDiv fails its test")
        END
      END
    END
  END
END testEuclidDiv;

PROCEDURE inverse (a, n : INTEGER) : INTEGER;
  VAR t, newt   : INTEGER;
  VAR r, newr   : INTEGER;
  VAR quotient  : INTEGER;
  VAR tmp       : INTEGER;
BEGIN
  t := 0;  newt := 1;
  r := n;  newr := a;
  WHILE newr # 0 DO
    quotient := euclidDiv (r, newr);
    tmp := newt;  newt := t - (quotient * newt);  t := tmp;
    tmp := newr;  newr := r - (quotient * newr);  r := tmp
  END;
  IF r > 1 THEN
    t := -1
  ELSIF t < 0 THEN
    t := t + n
  END;
  RETURN t
END inverse;

BEGIN
  testEuclidDiv;
  Out.Int (inverse (42, 2017), 0);
  Out.Ln
END modularInverseInOberon2.

Output:

$ obc modularInverseInOberon2.Mod && ./a.out
1969

ObjectIcon

Translation of: Oberon-2

Translation of: Mercury

# -*- ObjectIcon -*-

import exception
import io

procedure main ()
  test_euclid_div ()
  io.write (inverse (42, 2017))
end

procedure inverse (a, n)        # FAILS if there is no inverse.
  local t, newt, r, newr, quotient, tmp

  if n <= 0 then throw ("non-positive modulus")
  t := 0;  newt := 1
  r := n;  newr := a
  while newr ~= 0 do
  {
    quotient := euclid_div (r, newr)
    tmp := newt;  newt := t - (quotient * newt);  t := tmp
    tmp := newr;  newr := r - (quotient * newr);  r := tmp
  }
  r <= 1 | fail
  return (if t < 0 then t + n else t)
end

procedure euclid_div (x, y)
  # This kind of integer division always gives a remainder between 0
  # and abs(y)-1, inclusive. Thus the remainder is always a LEAST
  # RESIDUE modulo abs(y). (If y is a positive modulus, then only the
  # floor division branch is used.)
  return \
    if 0 <= y then              # Do floor division.
      (if 0 <= x then x / y
      else if (-x) % y = 0 then -((-x) / y)
      else -((-x) / y) - 1)
    else                        # Do ceiling division.
      (if 0 <= x then -(x / (-y))
      else if (-x) % (-y) = 0 then ((-x) / (-y))
      else ((-x) / (-y)) + 1)
end

procedure test_euclid_div ()
  local x, y, q, r

  every x := -100 to 100 do
    every y := -100 to 100 & y ~= 0 do
    {
      q := euclid_div (x, y)
      r := x - (q * y)
      if r < 0 | abs (y) <= r then
        # A remainder was outside the expected range.
        throw ("Test of euclid_div failed.")
    }
end

Output:

$ oiscript modular-inverse-task-OI.icn
1969

OCaml

Translation of: C

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

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

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

Owl Lisp

Works with: Owl Lisp version 0.2.1

(import (owl math)
        (owl math extra))

(define (euclid-quotient x y)
  (if (<= 0 y)
      (cond ((<= 0 x) (quotient x y))
            ((zero? (remainder (negate x) y))
             (negate (quotient (negate x) y)))
            (else (- (negate (quotient (negate x) y)) 1)))
      (cond ((<= 0 x) (negate (quotient x (negate y))))
            ((zero? (remainder (negate x) (negate y)))
             (quotient (negate x) (negate y)))
            (else (+ (quotient (negate x) (negate y)) 1)))))

;; A unit test of euclid-quotient.
(let repeat ((x -100)
             (y -100))
  (cond ((= x 101) #t)
        ((= y 0) (repeat x (+ y 1)))
        ((= y 101) (repeat (+ x 1) -100))
        (else (let* ((q (euclid-quotient x y))
                     (r (- x (* q y))))
                (cond ((< r 0) (display "negative remainder\n"))
                      ((<= (abs y) r) (display "remainder too large\n"))
                      (else (repeat x (+ y 1))))))))

(define (inverse a n)
  (let repeat ((t 0) (newt 1)
               (r n) (newr a))
    (cond ((not (zero? newr))
           (let ((quotient (euclid-quotient r newr)))
             (repeat newt (- t (* quotient newt))
                     newr (- r (* quotient newr)))))
          ((< 1 r) #f)                  ; The inverse does not exist.
          ((negative? t) (+ t n))
          (else t))))

(display (inverse 42 2017))
(newline)

Output:

$ ol modular-inverse-task-Owl.scm
1969

PARI/GP

Mod(1/42,2017)

Pascal

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

Perl

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.

Phix

Translation of: C

function mul_inv(integer a, n)
    if n<0 then n = -n end if
    if a<0 then a = n - mod(-a,n) end if
    integer t = 0,  nt = 1,
            r = n,  nr = a;   
    while nr!=0 do
        integer q = floor(r/nr)
        {t, nt} = {nt, t-q*nt}
        {r, nr} = {nr, r-q*nr}
    end while
    if r>1 then return "a is not invertible" end if
    if t<0 then t += n end if
    return t
end function
 
?mul_inv(42,2017)
?mul_inv(40, 1)
?mul_inv(52, -217)  /* Pari semantics for negative modulus */
?mul_inv(-486, 217)
?mul_inv(40, 2018)

Output:

1969
0
96
121
"a is not invertible"

PHP

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:

PicoLisp

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

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

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>

Prolog

egcd(_, 0, 1, 0) :- !.
egcd(A, B, X, Y) :-
    divmod(A, B, Q, R),
    egcd(B, R, S, X),
    Y is S - Q*X.

modinv(A, B, N) :-
    egcd(A, B, X, Y),
    A*X + B*Y =:= 1,
    N is X mod B.

Output:

?- modinv(42, 2017, N).
N = 1969.

?- modinv(42, 64, X).
false.

PureBasic

Using brute force.

EnableExplicit
Declare main()
Declare.i mi(a.i, b.i)

If OpenConsole("MODULAR-INVERSE")
  main() : Input() : End
EndIf

Macro ModularInverse(a, b)  
  PrintN(~"\tMODULAR-INVERSE(" + RSet(Str(a),5) + "," + 
         RSet(Str(b),5)+") = " + 
         RSet(Str(mi(a, b)),5))   
EndMacro

Procedure main()
  ModularInverse(42, 2017)  ; = 1969
  ModularInverse(40, 1)     ; = 0
  ModularInverse(52, -217)  ; = 96
  ModularInverse(-486, 217) ; = 121
  ModularInverse(40, 2018)  ; = -1
EndProcedure

Procedure.i mi(a.i, b.i)
  Define x.i = 1,
         y.i = Int(Abs(b)),
         r.i = 0  
  If y = 1 : ProcedureReturn 0 : EndIf  
  While x < y
    r = (a * x) % b   
    If r = 1 Or (y + r) = 1
      Break
    EndIf
    x + 1
  Wend
  If x > y - 1 : x = -1 : EndIf  
  ProcedureReturn x  
EndProcedure

Output:

        MODULAR-INVERSE(   42, 2017) =  1969
        MODULAR-INVERSE(   40,    1) =     0
        MODULAR-INVERSE(   52, -217) =    96
        MODULAR-INVERSE( -486,  217) =   121
        MODULAR-INVERSE(   40, 2018) =    -1

Python

Builtin function

Since python3.8, builtin function "pow" can be used directly to compute modular inverses by specifying an exponent of -1:

>>> pow(42, -1, 2017)
1969

Iteration and error-handling

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

Recursion and an option type

Or, using functional composition as an alternative to iterative mutation, and wrapping the resulting value in an option type, to allow for the expression of computations which establish the absence of a modular inverse:

from functools import (reduce)
from itertools import (chain)


# modInv :: Int -> Int -> Maybe Int
def modInv(a):
    return lambda m: (
        lambda ig=gcdExt(a)(m): (
            lambda i=ig[0]: (
                Just(i + m if 0 > i else i) if 1 == ig[2] else (
                    Nothing()
                )
            )
        )()
    )()


# gcdExt :: Int -> Int -> (Int, Int, Int)
def gcdExt(x):
    def go(a, b):
        if 0 == b:
            return (1, 0, a)
        else:
            (q, r) = divmod(a, b)
            (s, t, g) = go(b, r)
        return (t, s - q * t, g)
    return lambda y: go(x, y)


#  TEST ---------------------------------------------------

# Numbers between 2010 and 2015 which do yield modular inverses for 42:

# main :: IO ()
def main():
    print (
        mapMaybe(
            lambda y: bindMay(modInv(42)(y))(
                lambda mInv: Just((y, mInv))
            )
        )(
            enumFromTo(2010)(2025)
        )
    )

# -> [(2011, 814), (2015, 48), (2017, 1969), (2021, 1203)]


# GENERIC ABSTRACTIONS ------------------------------------


# enumFromTo :: Int -> Int -> [Int]
def enumFromTo(m):
    return lambda n: list(range(m, 1 + n))


# bindMay (>>=) :: Maybe  a -> (a -> Maybe b) -> Maybe b
def bindMay(m):
    return lambda mf: (
        m if m.get('Nothing') else mf(m.get('Just'))
    )


# Just :: a -> Maybe a
def Just(x):
    return {'type': 'Maybe', 'Nothing': False, 'Just': x}


# mapMaybe :: (a -> Maybe b) -> [a] -> [b]
def mapMaybe(mf):
    return lambda xs: reduce(
        lambda a, x: maybe(a)(lambda j: a + [j])(mf(x)),
        xs,
        []
    )


# maybe :: b -> (a -> b) -> Maybe a -> b
def maybe(v):
    return lambda f: lambda m: v if m.get('Nothing') else (
        f(m.get('Just'))
    )


# Nothing :: Maybe a
def Nothing():
    return {'type': 'Maybe', 'Nothing': True}


# MAIN ---
main()

Output:

[(2011, 814), (2015, 48), (2017, 1969), (2021, 1203)]

Quackery

Translation of: Forth

  [ dup 1 != if
      [ tuck 1 0
        [ swap temp put
          temp put
          over 1 > while
          tuck /mod swap
          temp take tuck *
          temp take swap -
          again ]
        2drop
        temp release
        temp take
        dup 0 < if
          [ over + ] ]
    nip ]                 is modinv ( n n --> n )

  42 2017 modinv echo

Output:

Using Extended Euclidean Algorithm

Handles negative args. Returns -1 for non-coprime args.

  [ dup 0 = iff
      [ 2drop 1 0 ]
      done
    dup unrot /mod
    dip swap recurse
    tuck 2swap *
    dip swap - ]        is egcd    ( n n --> n n )

  [ dup 0 < if negate
    over 0 < if
      [ swap negate
        over tuck mod
        - swap ]
    dup rot 2dup egcd
    2swap 2over rot *
    unrot * + 1 != iff
      [ drop 2drop -1 ]
      done
    nip swap mod ]      is modinv ( n n --> n   )

say "  42 2017 modinv --> "   42 2017 modinv echo cr ( --> 1969 )
say "  40    1 modinv --> "   40    1 modinv echo cr ( --> 0    )
say "  52 -217 modinv --> "   52 -217 modinv echo cr ( --> 96   )
say "-486  217 modinv --> " -486  217 modinv echo cr ( --> 121  )
say "  40 2018 modinv --> "   40 2018 modinv echo cr ( --> -1   )

Output:

  42 2017 modinv --> 1969
  40    1 modinv --> 0
  52 -217 modinv --> 96
-486  217 modinv --> 121
  40 2018 modinv --> -1

Racket

(require math)
(modular-inverse 42 2017)

Output:

Raku

(formerly Perl 6)

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)

REXX

/*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                                     /*initialize modular inverse to unity. */
        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                   /*Negative?  Then add the original  B. */
        return $

output when using the default inputs of: 42 2017

modular inverse of  42  by  2017  ───►  1969

Ring

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

RPL

Using complex numbers allows to ‘parallelize’ calculations and keeps the stack depth low: never more than 4 levels despite the simultaneous use of 6 variables: r, r’, u, u’, q - and b for the final touch.

Works with: Halcyon Calc version 4.2.7

RPL code

Comment

 ≪ 
   DUP ROT 1 R→C ROT 0 R→C
   WHILE DUP RE REPEAT
      OVER RE OVER RE / FLOOR
      OVER * NEG ROT +
   END
   DROP C→R ROT MOD
   SWAP 1 == SWAP 0 IFTE
≫ ‘MODINV’ STO

MODINV ( a b -- x )  with ax = 1 mod b
3: b   2: (r,u)←(a,1)   1:(r',u')←(b,0)
While r' ≠ 0
    q ← r // r'
   (r - q*r', u - q*u')

Forget (r',u') and calculate u mod b
Test r and return zero if a and b are not co-prime

Input:

123 456 MODINV
42 2017 MODINV

Output:

2: 1969
1: 0

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 'The maths are broken!'
  end
  x % et
end

 
> invmod(42,2017)
=> 1969

Simplified equivalent implementation

def modinv(a, m) # compute a^-1 mod m if possible
  raise "NO INVERSE - #{a} and #{m} not coprime" unless a.gcd(m) == 1
  return m if m == 1
  m0, inv, x0 = m, 1, 0
  while a > 1
    inv -= (a / m) * x0
    a, m = m, a % m
    inv, x0 = x0, inv
  end
  inv += m0 if inv < 0
  inv
end

> modinv(42,2017)
=> 1969

The OpenSSL module has modular inverse built in:

require 'openssl'
p OpenSSL::BN.new(42).mod_inverse(2017).to_i

Run BASIC

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:

Rust

fn mod_inv(a: isize, module: isize) -> isize {
  let mut mn = (module, a);
  let mut xy = (0, 1);
  
  while mn.1 != 0 {
    xy = (xy.1, xy.0 - (mn.0 / mn.1) * xy.1);
    mn = (mn.1, mn.0 % mn.1);
  }
  
  while xy.0 < 0 {
    xy.0 += module;
  }
  xy.0
}

fn main() {
  println!("{}", mod_inv(42, 2017))
}

Output:

Alternative implementation

fn modinv(a0: isize, m0: isize) -> isize {
    if m0 == 1 { return 1 }

    let (mut a, mut m, mut x0, mut inv) = (a0, m0, 0, 1);

    while a > 1 {
        inv -= (a / m) * x0;
        a = a % m;
        std::mem::swap(&mut a, &mut m);
        std::mem::swap(&mut x0, &mut inv);
    }
    
    if inv < 0 { inv += m0 }
    inv
}

fn main() {
  println!("{}", modinv(42, 2017))
}

Output:

Scala

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

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: [1]

Sidef

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:

Swift

extension BinaryInteger {
  @inlinable
  public func modInv(_ mod: Self) -> Self {
    var (m, n) = (mod, self)
    var (x, y) = (Self(0), Self(1))

    while n != 0 {
      (x, y) = (y, x - (m / n) * y)
      (m, n) = (n, m % n)
    }

    while x < 0 {
      x += mod
    }

    return x
  }
}

print(42.modInv(2017))

Output:

Tcl

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

Tiny BASIC

2017 causes integer overflow, so I'll do the inverse of 42 modulo 331 instead.

    PRINT "Modular inverse."
    PRINT "What is the modulus?"
    INPUT M
    PRINT "What number is to be inverted?"
    INPUT X
    PRINT "Solution is:"
10  LET A = A + 1
    GOTO 20
15  IF B = 1 THEN PRINT A
    IF B = 1 THEN END
    IF A = M-1 THEN PRINT "nonexistent"
    IF A = M-1 THEN END
    GOTO 10
20  LET B = A*X
30  IF B < M THEN GOTO 15
    LET B = B - M
    GOTO 30

Output:

Modular inverse.
What is the modulus?
331
What number is to be inverted?
42
Solution is:
134

Another version:

Translation of: GW-BASIC

    REM Modular inverse
    LET E=42
    LET T=2017
    GOSUB 100
    PRINT M
    END

    REM Increments E S (step) times until B is greater than T
    REM Repeats until B = 1 (MOD = 1) and C (count)
    REM B will not be greater than T + E
100 LET D=0
    IF E>=T THEN GOTO 130
    LET B=E
    REM At least one iteration is necessary
    LET S=((T-B)/E)+1
    LET B=B+S*E
    LET C=1+S
    LET B=B-T
110 IF B=1 THEN GOTO 120
    LET S=((T-B)/E)+1
    LET B=B+S*E
    LET C=C+S
    LET B=B-T
    GOTO 110
120 LET D=C
130 LET M=D
    RETURN

Output:

tsql

;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

TypeScript

Translation of: Pascal

// Modular inverse

function modInv(e: number, t: number): number {
  var d = 0;
  if (e < t) {
    var count = 1;
    var bal = e;
    do {
      var step = Math.floor((t - bal) / e) + 1;
      bal += step * e;
      count += step;
      bal -= t;
    } while (bal != 1);
    d = count;
  }
  return d;
}

console.log(`${modInv(42, 2017)}`); // 1969

Output:

uBasic/4tH

Translation of: C

Print FUNC(_MulInv(42, 2017))
End

_MulInv Param(2)
  Local(5)

  c@ = b@
  f@ = 0
  g@ = 1

  If b@ = 1 Then Return

  Do While a@ > 1
    e@ = a@ / b@
    d@ = b@
    b@ = a@ % b@
    a@ = d@

    d@ = f@
    f@ = g@ - e@ * f@
    g@ = d@
  Loop

  If g@ < 0 Then g@ = g@ + c@
Return (g@)

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

UNIX Shell

Works with: Bourne Again SHell

Works with: Korn Shell

Works with: Zsh

Translation of: PowerShell

function invmod {
  typeset -i a=$1 n=$2
  if (( n < 0 )); then (( n = -n )); fi
  if (( a < 0 )); then (( a = n - (-a) % n )); fi

  typeset -i t=0 nt=1 r=n nr q tmp
  (( nr = a % n ))
  while ((  nr )); do
    (( q = r/nr ))
    (( tmp = nt ))
    (( nt = t - q*nt ))
    (( t = tmp ))
    (( tmp = nr ))
    (( nr = r - q*nr ))
    (( r = tmp ))
  done
  if (( r > 1 )); then
    return 1
  fi
  while (( t < 0 )); do (( t += n )); done
  printf '%s\n' "$t"
}

invmod 42 2017

Output:

VBA

Translation of: Phix

Private Function mul_inv(a As Long, n As Long) As Variant
    If n < 0 Then n = -n
    If a < 0 Then a = n - ((-a) Mod n)
    Dim t As Long: t = 0
    Dim nt As Long: nt = 1
    Dim r As Long: r = n
    Dim nr As Long: nr = a
    Dim q As Long
    Do While nr <> 0
        q = r \ nr
        tmp = t
        t = nt
        nt = tmp - q * nt
        tmp = r
        r = nr
        nr = tmp - q * nr
    Loop
    If r > 1 Then
        mul_inv = "a is not invertible"
    Else
        If t < 0 Then t = t + n
        mul_inv = t
    End If
End Function
Public Sub mi()
    Debug.Print mul_inv(42, 2017)
    Debug.Print mul_inv(40, 1)
    Debug.Print mul_inv(52, -217) '/* Pari semantics for negative modulus */
    Debug.Print mul_inv(-486, 217)
    Debug.Print mul_inv(40, 2018)
End Sub

Output:

 1969 
 0 
 96 
 121 
a is not invertible

V (Vlang)

fn main() {
	println("42 %! 2017 = ${mult_inv(42, 2017)}")
}

fn mult_inv(aa int, bb int) int {
	mut a, mut b := aa, bb
	mut x0, mut t := 0, 0
	mut b0 := b  
	mut x1 := 1
	if b == 1 {return 1}
	for 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
}

Output:

42 %! 2017 = 1969

Wren

Library: Wren-big

import "./big" for BigInt

var a = BigInt.new(42)
var b = BigInt.new(2017)
System.print(a.modInv(b))

Output:

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

Output:

zkl

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