6

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

Все алгоритмы, что я нашёл в интернете, имеют переполнение, если числа близки к максимально поддерживаемым (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битные типы, поддерживаемые некоторыми компиляторами. Минус - не по стандарту.

15
  • Попробуйте использовать тип double или long double. Они вмещают гораздо большие значения, чем long long. Вот здесь можно посмотреть из диапазон
    – Стас
    29 июн 2020 в 8:09
  • @Стас тогда потеря точности будет, что тоже даст неверный результат. Да и по скорости надо потестить, думаю будет значительно медленнее
    – helloworld
    29 июн 2020 в 8:11
  • Главное - mod какое? такое, что mod*mod дает переполнение?...
    – Harry
    29 июн 2020 в 8:19
  • @Harry да, все 3 числа в квадрате могут давать переполнение. Однако, конкретно на моих данных, 2*a, 2*n, 2*mod не могут дать переполнение, т.е. 1 раз сложить по модулю само с собой можно любое число
    – helloworld
    29 июн 2020 в 8:21
  • Тогда надо искать 128-битную быструю реализацию вычисления остатка... Умножение пишется быстро, а вот деление... :(
    – Harry
    29 июн 2020 в 8:24

2 ответа 2

3

Всё очень просто ;) Быстрое возведение в степень. Умножение по модулю без переполнений. Сумма модульная без ошибок. И вуаля : быстро и правильно. Не используются операции умножения и деления. Сложность алгоритма 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)).

4
  • Разве это не вместо умножений аналогично бинарному возведению в степень использовать бинарное умножение, а значит тогда возведение в степень работает очень медленно?
    – вася
    2 июл 2020 в 16:23
  • Скорость очень хорошая. Посчитайте, проверьте. O(Log(степень)*Log(число)) @вася
    – AlexGlebe
    2 июл 2020 в 16:28
  • А у простого бинарного возведения в степень O(log(степень)).
    – вася
    2 июл 2020 в 16:31
  • Все проблемы начнутся при умножении больших чисел. Сложность будет расти порядка O((бит^2)/1000). Тогда время будет O(log(степень)*(бит^2)/1000). Практически разница с O(Log(степень)*бит) будет незаметна, но в перспективе и в лёгкости алгоритма (умножение займет всю страницу) логарифмическая сложность лучше. @вася
    – AlexGlebe
    3 июл 2020 в 4:47
0

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

#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);
  }
0

Ваш ответ

Нажимая на кнопку «Отправить ответ», вы соглашаетесь с нашими условиями использования и подтверждаете, что прочитали и поняли наши политику конфиденциальности и нормы поведения.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.