Modular inverse

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

$a^{-1}\equiv x{\pmod {m}}.$

Or in other words, such that:

$\exists k\in \mathbf {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.

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

Ada

<lang Ada>with Ada.Text_IO;

procedure Mod_Inv is

  procedure X_GCD(A, B: in Natural; D, X, Y: out Integer) is
     -- the Extended Euclidean Algorithm
     -- finds (D, X, Y) with D = GCD(A, B) = A*X + B*Y
     R: Natural := A mod B;
  begin
     if R=0 then
        D := B;
        X := 0;
        Y := 1;
     else
        X_GCD(B, R, D, Y, X);
        Y := Y - (A/B)*X;
     end if;
  end X_GCD;

  function Inverse(A, M: Integer) return Integer is
     -- computes the multiplicative inverse of A mod M, using X_GCD
     Result, GCD, Dummy: Integer;
  begin
     X_GCD(A, M, GCD, Result, Dummy);
     if GCD /= 1 then -- inverse does not exist!
        raise Constraint_Error with
          "GCD (" & Integer'Image(A) & "," & Integer'Image(M) & " ) =" &
          Integer'Image(GCD) & " /= 1";
     else -- make sure Result is in {0, ..., M-1}
        if Result < 0 then
           return Result+M;
        else
           return Result;
        end if;
     end if;
  end Inverse;

begin

  Ada.Text_IO.Put_Line(Natural'Image(Inverse(42, 2017)));
  -- Ada.Text_IO.Put_Line(Natural'Image(Inverse(154, 3311)));
  -- The above would raise CONSTRAINT_ERROR : GCD ( 154, 3311 ) = 77 /= 1

end Mod_Inv;</lang>

C

<lang c>#include <stdio.h>

int mul_inv(int a, int b) { int b0 = b, t, q; int x0 = 0, x1 = 1; if (b == 1) return 1; while (a > 1) { q = a / b; t = b, b = a % b, a = t; t = x0, x0 = x1 - q * x0, x1 = t; } if (x1 < 0) x1 += b0; return x1; }

int main(void) { printf("%d\n", mul_inv(42, 2017)); return 0; }</lang>

D

Translation of: C

<lang d>T modInverse(T)(T a, T b) pure nothrow {

   if (b == 1)
       return 1;
   T b0 = b, x0 = 0, x1 = 1;
   while (a > 1) {
       immutable q = a / b;
       auto t = b;
       b = a % b;
       a = t;
       t = x0;
       x0 = x1 - q * x0;
       x1 = t;
   }
   return (x1 < 0) ? (x1 + b0) : x1;

}

void main() {

   import std.stdio;
   writeln(modInverse(42, 2017));

}</lang>

Output:

1969

Go

The standard library function uses the extended Euclidean algorithm internally. <lang go>package main

import ( "fmt" "math/big" )

func main() { a := big.NewInt(42) m := big.NewInt(2017) k := new(big.Int).ModInverse(a, m) fmt.Println(k) }</lang>

Output:

1969

Haskell

<lang haskell>-- Extended Euclidean algorithm. Given non-negative a and b, return x, y and g -- such that ax + by = g, where g = gcd(a,b). Note that x or y may be negative. gcdExt a 0 = (1, 0, a) gcdExt a b = let (q, r) = a `quotRem` b

                (s, t, g) = gcdExt b r
            in (t, s - q * t, g)

-- Given a and m, return Just x such that ax = 1 mod m. If there is no such x -- return Nothing. modInv a m = let (i, _, g) = gcdExt a m

            in if g == 1 then Just (mkPos i) else Nothing
 where mkPos x = if x < 0 then x + m else x

main = do

 print $ 2 `modInv` 4
 print $ 42 `modInv` 2017</lang>

Output

Nothing
Just 1969

J

Solution:<lang j> modInv =: dyad def 'x y&|@^ <: 5 p: y'"0</lang> Example:<lang j> 42 modInv 2017 1969</lang> Notes: Calculates the modular inverse as a^( totient(m) - 1 ) mod m. Note that J has a fast implementation of modular exponentiation (which avoids the exponentiation altogether), invoked with the form m&|@^ (hence, we use explicitly-named arguments for this entry, as opposed to the "variable free" tacit style).

Java

The BigInteger library has a method for this: <lang java>System.out.println(BigInteger.valueOf(42).modInverse(BigInteger.valueOf(2017)));</lang>

Output:

1969

Mathematica

The built-in function FindInstance works well for this <lang Mathematica>modInv[a_, m_] :=

Block[{x,k}, x /. FindInstance[a x == 1 + k m, {x, k}, Integers]]</lang>

modInv[42, 2017]

{1969}

Perl

The modular inverse is not a perl builtin but there is a CPAN module who does the job.

<lang perl>use Math::ModInt qw(mod); print mod(42, 2017)->inverse</lang>

Output:

mod(1969, 2017)

Perl 6

<lang perl6>sub inverse($n, :$modulo) {

   my ($c, $d, $uc, $vc, $ud, $vd) = ($n % $modulo, $modulo, 1, 0, 0, 1);
   my $q;
   while $c != 0 {
       ($q, $c, $d) = ($d div $c, $d % $c, $c);
       ($uc, $vc, $ud, $vd) = ($ud - $q*$uc, $vd - $q*$vc, $uc, $vc);
   }
   return $ud < 0 ?? $ud + $modulo !! $ud;

}

say inverse 42, :modulo(2017)</lang>

Python

Implementation of this pseudocode with this. <lang python>>>> def extended_gcd(aa, bb):

   lastremainder, remainder = abs(aa), abs(bb)
   x, lastx, y, lasty = 0, 1, 1, 0
   while remainder:
       lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
       x, lastx = lastx - quotient*x, x
       y, lasty = lasty - quotient*y, y
   return lastremainder, lastx * (-1 if aa < 0 else 1), lasty * (-1 if bb < 0 else 1)

>>> def modinv(a, m): g, x, y = extended_gcd(a, m) if g != 1: raise ValueError return x % m

>>> modinv(42, 2017) 1969 >>> </lang>

REXX

<lang rexx>/*REXX program calcuates the modular inverse of an integer X modulo Y.*/ parse arg x y . /*get two integers from the C.L. */ say 'modular inverse of ' x " by " y ' ───► ' modInv(x,y) exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────MODINV subroutine───────────────────*/ modInv: parse arg a,b 1 ob; ox=0; $=1

if b \= 1 then do while a>1

                       parse value  a/b a//b b ox   with   q b a t
                       ox=$-q*ox;   $=trunc(t)
                       end    /*while a>1*/

if $<0 then $=$+ob return $</lang> output when using the input of: 42 2017

modular inverse of  42  by  2017  ───►  1969

Run BASIC

<lang runbasic>print multInv(42, 2017) end

function multInv(a,b) b0 = b multInv = 1 if b = 1 then goto [endFun] while a > 1 q = a / b t = b b = a mod b a = t t = x0 x0 = multInv - q * x0 multInv = int(t) wend if multInv < 0 then multInv = multInv + b0 [endFun] end function</lang>1969