Modular exponentiation without range restriction

Question

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

Answer 1

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.