# Compute modular inverse

Given two positive numbers x and n with x<2^n, write the shortest possible function to compute x^-1 mod 2^n. In other words, find y such that x*y=1 mod 2^n.

Your function must complete in a reasonable time for at least n=64, so exhaustive search will not work.

If the inverse does not exist, you must indicate that to the caller somehow (throw an exception, return a sentinel value, etc).

If you're wondering where to start, try the Extended Euclidean Algorithm.

• this is going to be a single statement in some math softwares – st0le Jan 29 '11 at 23:45
• @st0le: Right, and you wouldn't be allowed to use such a function in such systems. :-D – Chris Jester-Young Jan 30 '11 at 0:20

## Python 95 89

c is your function. Returns 0 if there is no inverse (i.e. when x is even).

p=lambda x,y,m:y and p(x,y/2,m)**2*x**(y&1)%m or 1
c=lambda x,n:[0,p(x,2**n-1,2**n)][x%2]

# Python, 29 bytes

lambda x,n:pow(x,2**n-1,2**n)

This returns 0 for even x. It uses Euler’s theorem, with the observation that 2^n − 1 is divisible by 2^(n − 1) − 1, via Python’s builtin fast modular exponentiation. This is plenty fast enough for n up to 7000 or so, where it starts taking more than about a second.

## Mathematica - 22

f=PowerMod[#,-1,2^#2]&

f[x,n] returns y with x*y=1 mod 2^n, otherwise x is not invertible modulo 2^n

## GolfScript (23 chars)

{:^((1${\.**2^?%}+*}:f; The sentinel result for a non-existent inverse is 0. This is a simple application of Euler's theorem. $x^{\varphi(2^n)} \equiv 1 \pmod {2^n}$, so $x^{-1} \equiv x^{2^{n-1}-1} \pmod {2^n}$ Unfortunately that's rather too big an exponential to compute directly, so we have to use a loop and do modular reduction inside the loop. The iterative step is $x^{2^k-1} = \left(x^{2^{k-1}-1}\right)^2 \times x$ and we have a choice of base case: either k=1 with {1\:^(@{\.**2^?%}+*}:f; or k=2 with {:^((1${\.**2^?%}+*}:f;

I'm working on another approach, but the sentinel is more difficult.

The key observation is that we can build the inverse up bit by bit: if $xy \equiv 1 \pmod{2^{k-1}}$ then $xy \in \{ 1, 1 + 2^{k-1} \} \pmod{2^k}$, and if $x$ is odd we have $x(y + xy-1) \equiv 1 \pmod{2^k}$. (If you're not convinced, check the two cases separately). So we can start at any suitable base case and apply the transformation $y' = (x+1)y - 1$ a suitable number of times.

Since $0x \equiv 1 \pmod {2^0}$ we get, by induction

$x\left(\frac{1 - (x+1)^n}{x}\right) \equiv 1 \pmod {2^n}$

where the inverse is the sum of a geometric sequence. I've shown the derivation to avoid the rabbit-out-of-a-hat effect: given this expression, it's easy to see that (given that the bracketed value is an integer, which follows from its derivation as a sum of an integer sequence) the product on the left must be in the right equivalence class if $x+1$ is even.

That gives the 19-char function

{1$)1$?@/~)2@?%}:f;

which gives correct answers for inputs which have an inverse. However, it's not so simple when $x$ is even. One potentially interesting option I've found is to add x&1 rather than 1.

{1$.1&+1$?@/~)2@?%}:f;

This seems to give sentinel values of either $0$ or $2^{n-1}$, but I haven't yet proved that.

Taking that one step further, we can ensure a sentinel of $0$ for even numbers by changing the expression $1 - (x+1)^n$ into $1 - 1^n$:

{1$.1&*)1$?@/~)2@?%}:f;

That ties with the direct application of Euler's theorem for code length, but is going to have worse performance for large $n$. If we take the arguments the other way round, as n x f, we can save one character and get to 22 chars:

{..1&*)2\$?\/~)2@?%}:f;

## Ruby - 88 characters

Use the function f.

def e a,b;a%b==0?[0,1]:(x,y=e(b,a%b);[y,x-(y*(a/b))])end
def f x,n;e(x,2**n)[0]*(x%2)end

Simply the recursive function from the linked wiki page, returns 0 on error.

• You can save some characters by inlining e: (e=->a,b{...})[x,2**n][0]. Can also save a character by testing a%b<1 instead of a%b==0. – histocrat Jun 12 '14 at 14:19

_!1=1
x!n|r<-x!div(n+1)2=(2-r*x)*rmod2^n

Using an algorithm based on Hensel’s lemma that doubles the number of digits in every iteration, this runs in under a second for n up to about 30 million!

# Pyth, 9 bytes

.^Et^2Q^2

Try it here!

Takes input in reverse order. Or, 9 bytes too: .^EtK^2QK.

# Explanation

.^Et^2Q^2  - Full program.

.^         - Pow function. The same in Python (pow).
E        - The second input.
^2Q    - And 2 ^ first input.
t       - Decremented.
^2  - And 2 ^ first input again.

# GAP, 39 bytes

f:=function(x,n)return 1/x mod 2^n;end;

f(x,n) returns the inverse of x modulo 2^n and gives an error message

Error, ModRat: for <r>/<s> mod <n>, <s>/gcd(<r>,<s>) and <n> must be coprime

if no inverse exists.