5
\$\begingroup\$

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
*/
\$\endgroup\$
  • \$\begingroup\$ Nice question. This was the main motivation of a question I asked earlier (i.e. Accurate modular arithmetic with double precision. \$\endgroup\$ – Joseph Wood Feb 10 '18 at 15:16
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ – chux Feb 10 '18 at 15:25
  • \$\begingroup\$ 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. \$\endgroup\$ – Toby Speight Feb 12 '18 at 11:23
  • \$\begingroup\$ @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. \$\endgroup\$ – chux Feb 12 '18 at 15:39
3
\$\begingroup\$

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.

\$\endgroup\$
  • 2
    \$\begingroup\$ Three multiplications, actually. \$\endgroup\$ – vnp Feb 12 '18 at 12:16
  • \$\begingroup\$ @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. \$\endgroup\$ – Toby Speight Feb 12 '18 at 12:55
  • \$\begingroup\$ Concerning self-test: See comment. IAC, the test code is only an ancillary review part. Yet I do appreciate your thoughts about it. \$\endgroup\$ – chux Feb 12 '18 at 15:43
  • \$\begingroup\$ 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. \$\endgroup\$ – chux Feb 12 '18 at 15:46
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Toby Speight Feb 19 '18 at 9:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.