Возведение в степень по модулю без переполнений

Question

Как возвести число в степень по модулю без переполнений в умножениях со скоростью не хуже бинарного возведения в степень, не используя расширения компиляторов?

Все алгоритмы, что я нашёл в интернете, имеют переполнение, если числа близки к максимально поддерживаемым (64битный unsigned long long у меня). Например, отсюда:

int binpow (int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1)
            res *= a;
        a *= a;
        n >>= 1;
    }
    return res;
}

при максимальных типах выражения res = res * a % mod; и a = a * a % mod дают ULL "переполнение" и результат неверный.

Пример кода на этом сайте, выдающий неправильные результаты.

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

Ну или длинная арифметика, не тестил скорость, но думаю, что будет не очень быстро.

Быстрый + без переполнений вариант - использовать встроенные 128битные типы, поддерживаемые некоторыми компиляторами. Минус - не по стандарту.

Answer 1

Всё очень просто ;) Быстрое возведение в степень. Умножение по модулю без переполнений. Сумма модульная без ошибок. И вуаля : быстро и правильно. Не используются операции умножения и деления. Сложность алгоритма A ^ P mod M : O(Log(P)*Log(M)). Пример на Си:

// gcc-9 -Wall -Wextra -Wpedantic -Winline -std=c11 -Os pow_mod.c -o pow_mod
# include <stdio.h>
static  inline  unsigned  long  long  int sum_mod (
  unsigned  long  long  int x ,
  unsigned  long  long  int y ,
  unsigned  long  long  int m ) {
  if ( m - x > y )
    return  x + y ;
  return  y - ( m - x ) ; }

// ( x * y ) mod m , 0 < x < m , y < m
// ( 5 * y ) mod m = ( y mod m ) + ( ( ( y * 2 ) mod m ) * 2 ) mod m
unsigned  long  long  int mul_mod (
  unsigned  long  long  int x ,
  unsigned  long  long  int y ,
  unsigned  long  long  int m ) {
  if  ( x > y ) {
    unsigned  long  long  int tmp = x ;
    x = y ;
    y = tmp ; }
  // x <= y
  unsigned  long  long  int res = 0 ;
  unsigned  long  long  int iy  = y ;
  while ( x ) {
    if ( x & 1 )
      res = sum_mod ( res , iy , m ) ;
    iy  = sum_mod  ( iy , iy , m ) ;
    x >>= 1 ; }
  return  res ; }

// ( x ^ n ) mod m , x < m , n > 0 , m > 0
// ( x ^ 5 ) mod m == ( ( x mod m ) * ( ( ( ( ( x ^ 2 ) mod m ) ) ^ 2 ) mod m ) ) mod m
unsigned  long  long  int pow_mod (
  unsigned  long  long  int x ,
  unsigned  long  long  int n ,
  unsigned  long  long  int m ) {
  unsigned  long  long  int res = 1 ;
  while ( n ) {
    if ( n & 1 )
      res = mul_mod ( res , x , m ) ;
    x = mul_mod ( x , x , m ) ;
    n >>= 1 ; }
  return res  ; }

int main(){
  printf("sizeof(unsigned  long  long  int)*8=%zu\n",sizeof(unsigned  long  long  int)*8);
  unsigned  long  long  int res ;
  res = pow_mod(10,10,11);
  printf("10^10 mod 11 = %llu\n",res);
  res = pow_mod(100,100,111);
  printf("100^100 mod 111 = %llu\n",res);
  res = pow_mod(2147483647,2147483647,2147483648);
  printf("2147483647^2147483647 mod 2147483648 = %llu\n",res);
  res = pow_mod(12345678900987654321ULL,12345678900987654321ULL,12345678901234567890ULL);
  printf("12345678900987654321^12345678900987654321 mod 12345678901234567890 = %llu\n",res);
  }

Результат :

> ./pow_mod 
sizeof(unsigned  long  long  int)*8=64
10^10 mod 11 = 1
100^100 mod 111 = 100
2147483647^2147483647 mod 2147483648 = 2147483647
12345678900987654321^12345678900987654321 mod 12345678901234567890 = 566503238319246891

Проверяем :

Когда нужно будет использовать длинные числа (4096 битные например) то алгоритму нужны будут операции +,-,>(больше),=(set присвоение),&1(oddp нечётность), >>=1(деление пополам). Это простые функции линейной сложности. А если использовать другой приём: использовать умножение и потом деление, то это приведёт к громодкости кода (если оно будет ваше) и замедлению (в классических алгоритмах сложность порядка O(n^2)).

Answer 2

Можно проще и быстрее, учитывая что суммы не переполняются:

#include <cstdio>

unsigned  long  long   int mul_mod (
  unsigned  long  long  int x ,
  unsigned  long  long  int y ,
  unsigned   long  long  int m ) {
  if  ( x > y ) {
    unsigned  long  long  int tmp = x ;
    x = y ;
    y = tmp ; }
  unsigned  long  long  int res = 0 ;
  unsigned  long  long  int iy  = y ;
  while ( x ) {
    if ( x & 1 )
      res =  (res  + iy) % m ;
    iy  = (iy + iy) % m;
    x >>= 1 ; }
  return  res ; }

unsigned  long  long  int pow_mod (
  unsigned  long  long  int x ,
  unsigned   long  long  int n ,
  unsigned  long  long  int m ) {
  unsigned  long  long  int res = 1 ;
  while ( n ) {
    if ( n & 1 )
      res = mul_mod ( res ,   x , m ) ;
    x = mul_mod  ( x , x , m ) ;
    n >>= 1 ; }
  return res  ; }

int main(){
  unsigned  long  long  int res ;
  res = pow_mod(10,10,11);
  printf("10^10 mod 11 = %llu\n",res);
  res = pow_mod(100,100,111);
  printf("100^100 mod 111 = %llu\n",res);
  res = pow_mod(2147483647,2147483647,2147483648);
  printf("2147483647^2147483647 mod 2147483648 = %llu\n",res);
  }

БЛОГ НА HUSL

Возведение в степень по модулю без переполнений

Ответы (2 шт):