# Modular exponentiation without range restriction

A typical modular exponentiation may be coded using the following algorithm.

powmod(x, expo, m) {
x = x mod m;
y = 1 mod m
while (expo > 0) {
if (is odd expo) {
y = (x * y) mod m;
}
expo /= 2;
x = (x * x) mod m;
}
return y;
}


Overflow may occur with x * y or x * x. This can only occur only if m*m is greater than the "maximum value + 1" of the type. To cope, the algorithm is amended with a test to call a function that can handle large values of m.

powmod(x, expo, m) {
if (m > square_root(max possible value + 1) {
return powmod_wider(x, expo, m);
}
x = x mod m;
....


# Design

The below set codes in C powmod(x,exp, m) functions for types: unsigned, unsigned long, unsigned long long, uintmax_t.

Each function calls a "wider" function when m is large. Should the widest math prove insufficient, a slower bit-by-bit version is called.

# Review goals

1. Function interface of powmod.h: Architecture considerations? (Like passing in mod_minus_1 instead of mod to allow a modulo range of [1 ... _MAX+1] vs. [0 ... _MAX]) What are some portable improvements?

2. Function implementation in powmod.c: How sensible and understandable? Coding the conditional *_THRESHOLD and *_NEXT macros I found a bit ugly and looking to improve. Was appending a u to simple constants useful concerning MISRA? Style review OK, but of secondary concern.

powmod.h

/*
* powmod.h
*  Created on: Jan 14, 2018
*      Author: chux
*/

#ifndef POWMOD_H_
#define POWMOD_H_

#include <stdint.h>

/*
* (x**expo) % mod
*
* The powmod functions compute x raised to the power expo, modded by mod.
*
* Valid for entire range of x, expo, mod, expect mod == 0.
*/
unsigned powmod(unsigned x, unsigned expo, unsigned mod);
unsigned long powmodl(unsigned long x, unsigned long expo, unsigned long mod);
unsigned long long powmodll(unsigned long long x, unsigned long long expo,
unsigned long long mod);
uintmax_t powmodmax(uintmax_t x, uintmax_t expo, uintmax_t mod);

#endif /* POWMOD_H_ */


powmod.c

/*
* powmod.c
*
*  Created on: Dec 1, 2018
*      Author: chux
*/

#include "powmod.h"
#include <limits.h>

/*
* Determine the function to call when wider math is warranted
*/
#if UINT_MAX <= ULONG_MAX/UINT_MAX
#define POWMOD_MOD_NEXT  powmodl
#elif UINT_MAX <= ULLONG_MAX/UINT_MAX
#define POWMOD_MOD_NEXT  powmodll
#elif UINT_MAX <= UINTMAX_MAX/UINT_MAX
#define POWMOD_MOD_NEXT  powmodmax
#else
#define POWMOD_MOD_NEXT  powmodmax_high
#endif

#if ULONG_MAX <= ULLONG_MAX/ULONG_MAX
#define POWMODL_MOD_NEXT  powmodll
#elif ULONG_MAX <= UINTMAX_MAX/ULONG_MAX
#define POWMODL_MOD_NEXT  powmodmax
#else
#define POWMODL_MOD_NEXT  powmodmax_high
#endif

#if ULLONG_MAX <= UINTMAX_MAX/ULLONG_MAX
#define POWMODLL_MOD_NEXT  powmodmax
#else
#define POWMODLL_MOD_NEXT  powmodmax_high
#endif

/*
* When mod > *_MOD_THRESHOLD, use a function that handles wider integer math.
* E. g. When UINT_MAX == 0xFFFFFFFF, POWMOD_MOD_THRESHOLD is 0x10000.
*/
#define POWMOD_MOD_THRESHOLD \
((UINT_MAX >> ((CHAR_BIT * sizeof (unsigned) + 1u)/2u)) + 1u)
#define POWMODL_MOD_THRESHOLD \
((ULONG_MAX >> ((CHAR_BIT * sizeof (unsigned long) + 1u)/2u)) + 1u)
#define POWMODLL_MOD_THRESHOLD \
((ULLONG_MAX >> ((CHAR_BIT * sizeof (unsigned long long) + 1u)/2u)) + 1u)
#define POWMODMAX_MOD_THRESHOLD \
((UINTMAX_MAX >> ((CHAR_BIT * sizeof (uintmax_t) + 1u)/2u)) + 1u)

/*
* (a+b)%mod
*/
static uintmax_t addmodmax(uintmax_t a, uintmax_t b, uintmax_t mod) {
uintmax_t sum = a + b;
if (sum < a) {
sum = (sum + 1u) % mod + UINTMAX_MAX % mod; // This addition does not overflow
}
return sum % mod;
}

/*
* (a*b)%mod
*/
static uintmax_t mulmodmax(uintmax_t a, uintmax_t b, uintmax_t mod) {
uintmax_t prod = 0;
while (b > 0) {
if (b % 2u) {
}
b /= 2u;
}
return prod;
}

/*
* power(a,b)%mod without resorting to wider math.
*/
static uintmax_t powmodmax_high(uintmax_t x, uintmax_t expo, uintmax_t mod) {
x %= mod;
uintmax_t y = mod > 1u; // 1u % mod;
while (expo > 0) {
if (expo % 2u) {
y = mulmodmax(x, y, mod);
}
expo /= 2u;
x = mulmodmax(x, x, mod);
}
return y;
}

/*
* See powmod.h
*/
unsigned powmod(unsigned x, unsigned expo, unsigned mod) {
if (mod > POWMOD_MOD_THRESHOLD) {
return (unsigned) POWMOD_MOD_NEXT(x, expo, mod);
}
x %= mod;
unsigned y = mod > 1u; // 1u % mod;
while (expo > 0) {
if (expo % 2u) {
y = (x * y) % mod;
}
expo /= 2u;
x = (x * x) % mod;
}
return y;
}

/*
* See powmod.h
*/
unsigned long powmodl(unsigned long x, unsigned long expo, unsigned long mod) {
if (mod > POWMODL_MOD_THRESHOLD) {
return (unsigned long) POWMODL_MOD_NEXT(x, expo, mod);
}
x %= mod;
unsigned long y = mod > 1u; // 1u % mod;
while (expo > 0) {
if (expo % 2u) {
y = (x * y) % mod;
}
expo /= 2u;
x = (x * x) % mod;
}
return y;
}

/*
* See powmod.h
*/
unsigned long long powmodll(unsigned long long x, unsigned long long expo,
unsigned long long mod) {
if (mod > POWMODLL_MOD_THRESHOLD) {
return (unsigned long long) POWMODLL_MOD_NEXT(x, expo, mod);
}
x %= mod;
unsigned long long y = mod > 1u; // 1u % mod;
while (expo > 0) {
if (expo % 2u) {
y = (x * y) % mod;
}
expo /= 2u;
x = (x * x) % mod;
}
return y;
}

/*
* See powmod.h
*/
uintmax_t powmodmax(uintmax_t x, uintmax_t expo, uintmax_t mod) {
if (mod > POWMODMAX_MOD_THRESHOLD) {
return powmodmax_high(x, expo, mod);
}
x %= mod;
uintmax_t y = mod > 1u; // 1u % mod;
while (expo > 0) {
if (expo % 2u) {
y = (x * y) % mod;
}
expo /= 2u;
x = (x * x) % mod;
}
return y;
}


Test code

#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include "powmod.h"

/*
* Test functions against basic math properties and against wider widths
* Test against a FP version.  Mis-matches here are not necessarily
* wrong for the integer function - just narrow the list of
* candidates to manually review.
*/
void powmod_test(unsigned x, unsigned expo, unsigned mod) {
unsigned y = powmod(x, expo, mod);
assert(y < mod);
assert(mod > 1 || y == 0);
unsigned long yl = powmodl(x, expo, mod);
assert(y == yl);
unsigned long long yll = powmodll(x, expo, mod);
assert(y == yll);
uintmax_t yj = powmodmax(x, expo, mod);
assert(y == yj);

long double yd = fmodl(powl(x, expo), mod);
if (isfinite(yd) && y != yd) {
printf("powmod(%x, %x, %x) --> %x %Le\n", x, expo, mod, y, yd);
fflush(stdout);
}
}

/*
* See powmod_test()
*/
void powmodmax_test(uintmax_t x, uintmax_t expo, uintmax_t mod) {
uintmax_t y = powmodmax(x, expo, mod);
fflush(stdout);
assert(y < mod);
assert(mod > 1u || y == 0);

long double yd = fmodl(powl(x, expo), mod);
if (isfinite(yd) && y != yd) {
printf("powmodmax(%jx, %jx, %jx) --> %jx %Le\n", x, expo, mod, y, yd);
fflush(stdout);
}
}

/*
* Exercise powmod() functions with various values.
*/
void powmod_tests(void) {
puts("Print out tests that failed checking against (double) math.");
puts("Inspect individually.");
unsigned u[] = {0, 1u, 2u, POWMOD_MOD_THRESHOLD - 1u, POWMOD_MOD_THRESHOLD,
POWMOD_MOD_THRESHOLD + 1u, UINT_MAX - 2u, UINT_MAX - 1u, UINT_MAX};
uintmax_t uj[] = {0, 1u, 2u, POWMODLL_MOD_THRESHOLD - 1u,
POWMODLL_MOD_THRESHOLD,
POWMODLL_MOD_THRESHOLD + 1u, UINTMAX_MAX - 2u, UINTMAX_MAX - 1u,
UINTMAX_MAX};
size_t n = sizeof u / sizeof u[0];
assert(n == sizeof uj / sizeof uj[0]);
for (size_t x = 0; x < n; x++) {
for (size_t e = 0; e < n; e++) {
for (size_t m = 0; m < n; m++) {
if (u[m]) {
powmod_test(u[x], u[e], u[m]);
}
if (uj[m]) {
powmodmax_test(uj[x], uj[e], uj[m]);
}
}
}
}
}

int main(void) {
powmod_tests();
puts("Done");
return 0;
}

/* Test results
Print out tests that failed checking against (double) math.
Inspect individually.
powmodmax(100000001, 2, 2) --> 1 0.000000e+00
powmodmax(100000001, 2, ffffffff) --> 4 3.000000e+00
powmodmax(100000001, 2, 100000000) --> 1 0.000000e+00
powmodmax(100000001, 2, 100000001) --> 0 4.294967e+09
powmodmax(100000001, 2, fffffffffffffffd) --> 200000004 8.589935e+09
powmodmax(100000001, 2, fffffffffffffffe) --> 200000003 8.589935e+09
powmodmax(100000001, 2, ffffffffffffffff) --> 200000002 8.589935e+09
powmodmax(fffffffffffffffd, 2, 2) --> 1 0.000000e+00
powmodmax(fffffffffffffffd, 2, ffffffff) --> 4 4.294967e+09
powmodmax(fffffffffffffffd, 2, 100000000) --> 9 0.000000e+00
powmodmax(fffffffffffffffd, 2, 100000001) --> 4 4.294967e+09
powmodmax(fffffffffffffffd, 2, fffffffffffffffd) --> 0 1.844674e+19
powmodmax(fffffffffffffffd, 2, fffffffffffffffe) --> 1 1.844674e+19
powmodmax(fffffffffffffffd, 2, ffffffffffffffff) --> 4 1.844674e+19
powmodmax(fffffffffffffffe, 2, ffffffff) --> 1 4.294967e+09
powmodmax(fffffffffffffffe, 2, 100000000) --> 4 0.000000e+00
powmodmax(fffffffffffffffe, 2, 100000001) --> 1 4.294967e+09
powmodmax(fffffffffffffffe, 2, fffffffffffffffd) --> 1 1.844674e+19
powmodmax(fffffffffffffffe, 2, fffffffffffffffe) --> 0 1.844674e+19
powmodmax(fffffffffffffffe, 2, ffffffffffffffff) --> 1 1.844674e+19
powmodmax(ffffffffffffffff, 2, 2) --> 1 0.000000e+00
powmodmax(ffffffffffffffff, 2, ffffffff) --> 0 4.294967e+09
powmodmax(ffffffffffffffff, 2, 100000000) --> 1 0.000000e+00
powmodmax(ffffffffffffffff, 2, 100000001) --> 0 4.294967e+09
powmodmax(ffffffffffffffff, 2, fffffffffffffffd) --> 4 3.000000e+00
powmodmax(ffffffffffffffff, 2, fffffffffffffffe) --> 1 0.000000e+00
powmodmax(ffffffffffffffff, 2, ffffffffffffffff) --> 0 1.844674e+19
Done
*/

• Nice question. This was the main motivation of a question I asked earlier (i.e. Accurate modular arithmetic with double precision. Commented Feb 10, 2018 at 15:16
• @JosephWood As C++, did not readily catch it, Will look into at as the C++ aspect is not the core part, but double math is. Commented Feb 10, 2018 at 15:25
• As I write more C++ than C these days, I read the whole code thinking how much duplication can disappear by using templates. My only criticism is that the tests aren't self-checking - I find that tests that rely on manual verification are less reliable. If I can find anything else to criticise, I'll write a proper review. Commented Feb 12, 2018 at 11:23
• @TobySpeight "tests aren't self-checking" --> As coded, most of the tests are self-checking - using simple code cross check and FP code. It is the remaining tests that are printed for manual review. I have found that writing test code can often exceed the code under test. Still yes, additional testing could be had. Commented Feb 12, 2018 at 15:39

# Does a zero modulus make sense?

Consider

  unsigned long y = mod > 1u; // 1u % mod;


There are two cases we're considering here. I'd say that mod == 0 is a run-time error, and we can return anything we like. If mod == 1, then the remainder will always be zero, and we can choose that for the invalid case too:

  if (mod <= 1u) {
return 0;
}


This saves us performing the arithmetic in these trivial cases.

# Tests should be self-checking

I find that tests that rely on manual verification less useful than self-checking tests. I recommend including the (exact) expected results in the test suite; we can then keep quiet about successful tests (reducing output clutter), print the failures to the standard error stream, and finally exit with a failure status unless all tests passed (return failure_count > 0).

# Alternative fall-back algorithm

We can multiply two uintmax_t values into a pair of uintmax_t results, by masking each input value into an upper and lower half and performing the four multiplications separately. Then we can multiply the upper result by (1+UINTMAX_MAX)%mod (or rather by(UINTMAX_MAX%mod+1)%mod) and add it to the lower result.

With care, we might even be able to use a uintmax_t[2] throughout the computation, but I haven't considered that in enough detail.

• Three multiplications, actually.
– vnp
Commented Feb 12, 2018 at 12:16
• @vnp: even without Karatsuba's algorithm, I should have recognised that we're computing a square, so we have (h+l)²==h² + 2hl + l² for only three multiplications. Commented Feb 12, 2018 at 12:55
• Concerning self-test: See comment. IAC, the test code is only an ancillary review part. Yet I do appreciate your thoughts about it. Commented Feb 12, 2018 at 15:43
• The Alternative fall-back algorithm has merit and is a good idea. The expansion of the idea may have missed some carry concerns, yet it is likely faster than the bit-wise loop of the original code. Commented Feb 12, 2018 at 15:46
• Perhaps I was too concise and cryptic there, @chux. Or perhaps I haven't thought it through completely. We do know that (1+UINTMAX_MAX)%mod must be less than mod`, but perhaps I was just pontificating on an idea without working it out properly. I might have to consider this a bit further. Commented Feb 19, 2018 at 9:34