Java inverse modulo 2**64

Question

Given an odd long x, I'm looking for long y such that their product modulo 2**64 (i.e., using the normal overflowing arithmetic) equals to 1. To make clear what I mean: This could be computed in a few thousand year this way:

for (long y=1; ; y+=2) {
    if (x*y == 1) return y;
}

I know that this can be solved quickly using the extended Euclidean algorithm, but it requires the ability to represent all the involved numbers (ranging up to 2**64, so even unsigned arithmetic wouldn't help). Using BigInteger would surely help, but I wonder if there's a simpler way, possibly using the extended Euclidean algorithm implemented for positive longs.

Hacker's Delight [suggests](http://www.hackersdelight.org/HDcode/mulinv.c.txt) an algorithm for mod 2^32. I'd try some variant of that, though I'd test heavily, of course. (Perhaps that might be worth including in Guava...) — Louis Wasserman, Jul 28 '12 at 16:22
@Louis Wasserman: Nice link... in the meantime I think I've got something 3x faster -- I'll post my results later. Btw., I needed `pow(long, long)` (missing from `LongMath`) for one approach. — maaartinus, Jul 29 '12 at 00:11
`pow(long, long)` is deliberately missing because taking anything to a power that doesn't fit into an `int` is essentially guaranteed to overflow. (Though I guess that's what you want in your case.) — Louis Wasserman, Jul 29 '12 at 08:37
Yes. For things like `checkedPow` or `saturatedPow`, I agree that long exponent makes no sense, for `pow` I do not. I published [benchmarks and tests](https://dl.dropbox.com/u/4971686/published/maaartin/math/index.html) for 6 different solution possibilities. Concerning Guava, I'd suggest including things from linked Math64. — maaartinus, Jul 29 '12 at 13:44
The rule we ended up going with is that we assume you don't want to deliberately cause overflow. The distinction between `pow` and `checkedPow` is not whether or not it's expected to overflow, but whether or not you want to pay the overhead of checking. — Louis Wasserman, Jul 29 '12 at 14:07

score 2 · Answer 1 · answered Jul 28 '12 at 17:57

In the meantime I've recalled/reinvented a very simple solution:

public static int inverseOf(int x) {
    Preconditions.checkArgument((x&1)!=0, "Only odd numbers have an inverse, got " + x);
    int y = 1;
    for (int mask=2; mask!=0; mask<<=1) {
        final int product = x * y;
        final int delta = product & mask;
        y |= delta;
    }
    return y;
}

It works because of two things:

in each iteration if the corresponding bit of product is 1, then it's wrong, and the only way to fix is is by changing the corresponding bit of y
no bit of y influences any less significant bit of product, so no previous work gets undone

I started with int since for long it must work too, and for int I could run an exhaustive test.

Another idea: there must a a number n>0 such that x**n == 1, and therefore y == x**(n-1). This should probably be faster, I just can't recall enough math to compute n.

+1 interesting answer. I think the value of `n` you mention is `2**62`. — Luke Woodward, Jul 28 '12 at 21:19
@Luke Woodward: Using [Euler's totient function](http://en.wikipedia.org/wiki/Euler%27s_totient_function) I get `2**63`, but it seems to work with `2**62` as well. And it seems to lead to the fastest method. — maaartinus, Jul 28 '12 at 23:46

score 1 · Accepted Answer · answered Jul 28 '12 at 17:28

Here's one way of doing it. This uses the extended Euclidean algorithm to find an inverse of abs(x) modulo 2⁶², and at the end it 'extends' the answer up to an inverse modulo 2⁶⁴ and applies a sign change if necessary:

public static long longInverse(long x) {

    if (x % 2 == 0) { throw new RuntimeException("must be odd"); }

    long power = 1L << 62;

    long a = Math.abs(x);
    long b = power;
    long sign = (x < 0) ? -1 : 1;

    long c1 = 1;
    long d1 = 0;
    long c2 = 0;
    long d2 = 1;

    // Loop invariants:
    // c1 * abs(x) + d1 * 2^62 = a
    // c2 * abs(x) + d2 * 2^62 = b 

    while (b > 0) {
        long q = a / b;
        long r = a % b;
        // r = a - qb.

        long c3 = c1 - q*c2;
        long d3 = d1 - q*d2;

        // Now c3 * abs(x) + d3 * 2^62 = r, with 0 <= r < b.

        c1 = c2;
        d1 = d2;
        c2 = c3;
        d2 = d3;
        a = b;
        b = r;
    }

    if (a != 1) { throw new RuntimeException("gcd not 1 !"); }

    // Extend from modulo 2^62 to modulo 2^64, and incorporate sign change
    // if necessary.
    for (int i = 0; i < 4; ++i) {
        long possinv = sign * (c1 + (i * power));
        if (possinv * x == 1L) { return possinv; }
    }

    throw new RuntimeException("failed");
}

I found it easier to work with 2⁶² than 2⁶³, mainly because it avoids problems with negative numbers: 2⁶³ as a Java long is negative.

I've accepted your solution since your comment to mine lead to the fastest one, see [here](https://dl.dropbox.com/u/4971686/published/maaartin/math/index.html). — maaartinus, Jul 29 '12 at 13:46

Java inverse modulo 2**64

2 Answers2

Linked