Модульное возведение в степень без ограничение диапазона


Типичный модульного возведения в степень может быть закодирована, используя следующий алгоритм.

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;
}

Переполнение может произойти с x * y или x * x. Это может только произойти, только если m*m превышает "максимальное значение + 1" типа. Чтобы справиться с этим, алгоритм дополнен тест, чтобы вызвать функцию, которая может обрабатывать большие значения m.

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

Дизайн

Ниже установить коды на C powmod(x,exp, m) функции для типов: unsigned, unsigned long, unsigned long long, uintmax_t.

Каждая функция вызывает "расширенная" функция, когда m большой. Должны широчайший математика оказаться недостаточно, более медленный бит-в-бит версия называется.


Цели обзора

  1. Функции интерфейса powmod.hархитектура соображения? (Как, проходя в mod_minus_1 вместо mod чтобы разрешить модулю спектр [1 ... _MAX+1] и [0 ... _MAX]) Каковы некоторые портативные улучшений?

  2. Функция внедрения в powmod.cкак толково и понятно? Кодирование условного *_THRESHOLD и *_NEXT макросы я нашел немного некрасиво и хотите улучшить. Был добавив u для простых констант полезное по Мишра? Стиль обзора ОК, но второстепенной задачей.


powmod.ч

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

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

Тестовый код

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


223
4
задан 10 февраля 2018 в 03:02 Источник Поделиться
Комментарии
1 ответ

Не нулевой модуль имеет смысла?

Рассмотрим

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

В двух случаях мы рассматриваем здесь. Я бы сказал, что mod == 0 это ошибка времени выполнения, и мы можем вернуть все, что нам нравится. Если mod == 1тогда остаток всегда будет ноль, и мы можем выбрать, что за недопустимый случае:

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

Это экономит нам выполнять арифметические операции в этих тривиальных случаях.

Тесты должны быть самоконтроль

Я считаю, что тесты, которые полагаются на проверки менее полезны, чем самопроверки тестов. Я рекомендую в том числе (точный) ожидаемые результаты в ходе тестов, мы можем потом молчать об успешных тестах (уменьшение выходной беспорядок), печатать неудачи стандартного потока ошибок, и, наконец, выход с ошибкой статус, если все тесты прошли успешно (return failure_count > 0).

Альтернативные запасной алгоритм

Мы можем перемножить два uintmax_t значения в пару uintmax_t результаты, маскируя каждого входного значения в верхней и нижней половине и выполнения четырех умножений отдельно. Тогда мы можем умножить на верхние результата (1+UINTMAX_MAX)%mod (или, скорее,(UINTMAX_MAX%mod+1)%mod) и добавьте ее в более низкий результат.

С осторожностью, мы могли бы даже быть в состоянии использовать uintmax_t[2] по всей выкладке, но я не считала, что достаточно подробно.

2
ответ дан 12 февраля 2018 в 11:02 Источник Поделиться