We’ve previously explored the Extended Euclidean algorithm, and it’s easy to use a special case of it to implement the modular multiplicative inverse. We’ll start by reproducing the final function from an older post that derived a correct and efficient implementation for the Extended Euclidean algorithm, using inputs with an unsigned integral type:

```
template <class S, class U>
void extended_euclidean(U a, U b, U& gcd, S& x, S& y)
{
static_assert(std::numeric_limits<S>::is_integer);
static_assert(std::numeric_limits<S>::is_signed);
static_assert(std::numeric_limits<U>::is_integer);
static_assert(!(std::numeric_limits<U>::is_signed));
static_assert(std::is_same<typename std::make_signed<U>::type, S>::value);
S x1=1, y1=0, x0=0, y0=1;
U a1=a, a2=b, q=0;
while (a2 != 0) {
S x2 = x0 - static_cast<S>(q)*x1;
S y2 = y0 - static_cast<S>(q)*y1;
x0=x1; y0=y1;
U a0=a1;
x1=x2; y1=y2; a1=a2;
q = a0/a1;
a2 = a0 - q*a1;
}
x = x1;
y = y1;
gcd = a1;
}
```

*Figure 1: Extended Euclidean Algorithm, Unsigned Integer Inputs.*

There’s something particularly interesting about Figure 1, given that it computes the greatest common divisor **d** and the Bézout coefficients **x** and** y**. By starting from the definition of the Bézout coefficients and then taking the remainder of both sides modulo **a** (which requires **a** > 0), we can see that

a*x + b*y == gcd(a,b) == d

a*x + b*y == d

a*x + b*y ≡ d (mod a)

b*y ≡ d (mod a)

If **a** and **b** are coprime, then gcd(a,b) == d == 1, and continuing from above,

b*y ≡ 1 (mod a)

So under this condition, **y** is the modular multiplicative inverse (modulo **a**) of the input **b**.

Without knowing whether **a** and **b** are coprime, let’s assume some integer **z** exists such that

b*z ≡ 1 (mod a)

By the definition of congruence, this would mean that some integer **k** exists such that

b*z + a*k = 1

As above, let **d** be the greatest common divisor of **a** and **b**. By definition **d** divides both **a** and **b**, and so **d** divides b*z + a*k = 1, and thus **d** divides 1. Since **d** divides 1, **d **must be either -1 or 1. And since the greatest common divisor can not be negative, this means **d** is 1. Thus, if an integer **z** exists that is the modular multiplicative inverse (modulo **a**) of **b**, then **a** and **b** are coprime. Putting this together with the paragraph above, we know that the modular multiplicative inverse (modulo **a**) of **b** exists if and only if **a** and **b** are coprime.

So let us assume that **a** and **b** are coprime, and that therefore the inverse exists. We just showed for this case that **y** is the modular multiplicative inverse (modulo **a**) of the input **b**. And we can see that Figure 1 calculates and assigns the value of **y** upon completion. We can also see in Figure 1 that **x** and all the calculations of **x0**, **x1**, **x2** go completely unused when we are finding the modular multiplicative inverse. We noted that we must require **a** > 0 since we are using **a** for the modulo operation, but we will take the requirement further and assert that **a **> 1, since it makes our final work easier and a modulus of 1 would be pointless anyway (everything would be always congruent to 0). We can use these insights to change Figure 1 into a function to calculate the modular multiplicative inverse:

```
template <class S, class U>
void modular_multiplicative_inverse(U a, U b, U& gcd, S& y)
{
static_assert(std::numeric_limits<S>::is_integer);
static_assert(std::numeric_limits<S>::is_signed);
static_assert(std::numeric_limits<U>::is_integer);
static_assert(!(std::numeric_limits<U>::is_signed));
static_assert(std::is_same<typename std::make_signed<U>::type, S>::value);
assert(a > 1);
S y1=0, y0=1;
U a1=a, a2=b, q=0;
while (a2 != 0) {
S y2 = y0 - static_cast<S>(q)*y1;
y0=y1;
U a0=a1;
y1=y2;
a1=a2;
q = a0/a1;
a2 = a0 - q*a1;
}
y = y1;
gcd = a1;
}
```

*Figure 2: Modular Multiplicative Inverse (Can Be Improved).*

We can improve the function above by returning the modular multiplicative inverse rather than setting **y** to its value. We can mostly remove the gcd as well; its main value to a caller is that it lets a caller know if the input **a** and the **modulus** are coprime (if not, the modular multiplicative inverse of **a** would not exist). Instead, as a convenience to a client of the function we can return zero if the inverse does not exist, since a valid modular multiplicative inverse can never be zero. After making these changes, we have the following:

```
template <class S, class U>
S modular_multiplicative_inverse(U a, U b)
{
static_assert(std::numeric_limits<S>::is_integer);
static_assert(std::numeric_limits<S>::is_signed);
static_assert(std::numeric_limits<U>::is_integer);
static_assert(!(std::numeric_limits<U>::is_signed));
static_assert(std::is_same<typename std::make_signed<U>::type, S>::value);
assert(a > 1);
S y1=0, y0=1;
U a1=a, a2=b, q=0;
while (a2 != 0) {
S y2 = y0 - static_cast<S>(q)*y1;
y0=y1;
U a0=a1;
y1=y2;
a1=a2;
q = a0/a1;
a2 = a0 - q*a1;
}
S y = y1;
U gcd = a1;
if (gcd == 1)
return y;
else
return 0;
}
```

*Figure 3: Modular Multiplicative Inverse (Almost Done).*

This still can be improved. We might as well rename **a** to **modulus**, since that’s what it is. We will also rename **b** to **val** in the function interface. And returning a signed value is a bit strange, since it means we might be returning a negative value. We can note in passing that a negative return value can be correct despite its strangeness – so long as the returned value belongs to the congruence class of the answer, modulo the modulus, it is technically correct. Nevertheless, it would be more convenient for a client of the function if we return the smallest non-negative integer belonging to that congruence class and return it as type **U**. And since the function interface won’t use template type **S** anymore, we can drop it from the template parameters. These changes give us the final implementation:

```
template <class U>
U modular_multiplicative_inverse(U modulus, U val)
{
static_assert(std::numeric_limits<U>::is_integer);
static_assert(!(std::numeric_limits<U>::is_signed));
using S = typename std::make_signed<U>::type;
assert(modulus > 1);
S y1=0, y0=1;
U a1=modulus, a2=val, q=0;
while (a2 != 0) {
S y2 = y0 - static_cast<S>(q)*y1;
y0=y1;
U a0=a1;
y1=y2;
a1=a2;
q = a0/a1;
a2 = a0 - q*a1;
}
S y = y1;
U gcd = a1;
if (gcd == 1)
return (y >= 0) ? static_cast<U>(y) : static_cast<U>(y + modulus);
else
return 0;
}
```

*Figure 4: Final Modular Multiplicative Inverse.*

We do need to explain the conditional at the end, which returns either **y** or **y**+**modulus** depending on whether **y** is negative or not. How do we know that **y**+**modulus** will be non-negative? And how do we know that the returned value will be the smallest non-negative integer that belongs to the answer’s congruence class? An earlier blog post on bounds for the Extended Euclidean Algorithm had a section called Final Bounds, which referred to combined essential assertion results. One of those assertion results wasn’t mentioned in the blog entry since it was only a step toward the final bounds discussed by the blog entry, but for our purposes now that assertion result gives us the information we need. It showed that at the end of the Extended Euclidean function,**abs(y0) <= max(1,a/2)**

Figure 4 is derived from the unsigned version of the Extended Euclidean algorithm, which in turn is derived from the normal Extended Euclidean implementation. So this bound applies to it, but we need to account for changes in variable names. **y0** from the essential asserts file corresponds to **y** in Figure 4, and **a** corresponds to **modulus**, and so we have in Figure 4 the final bounds**abs(y) <= max(1,modulus/2)**

Figure 4 asserts at its start that **modulus** > 1, which means **modulus/2** will be greater than or equal to 1, and thus we can omit the max() portion, giving us our final bounds as**abs(y) <= modulus/2**

Another way to state this is

–**modulus/2** <= **y <= modulus/2**

It is easy to see that if **y** is negative, adding **modulus** to it will provide a non-negative result that is less than **modulus**. And if **y** is non-negative, then **y** is obviously less than **modulus**. Thus our returned result will be the smallest non-negative integer that belongs to the answer’s congruence class.

We’ve now got a nice implementation of the modular multiplicative inverse, derived from the unsigned version of the Extended Euclidean algorithm.

Pingback: A Faster Multiplicative Inverse (Mod A Power Of 2) | Jeffrey Hurchalla