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) 
 prod = addmodmax(prod, a, mod);
 
 b /= 2u;
 a = addmodmax(a, a, mod);
 
 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) 

/*
 * 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 

/*
 * 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
*/

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.