Оптимизация вычисления остатка по модулю

[Dmitriy40]

Камрады, всем привет, такой вопрос: а известна модификация алгоритма Баррета с линейным временем работы от длины делимого $n$ (при малом делителе $p$, помещающимся в регистр или встроенный тип данных)?

Обычно там же требуется умножение на длинный (длиной с делимое) обратный элемент к $p$, т.е. на $\lfloor 2^n/p \rfloor$, для получения частного, а это по идее квадратичная операция $(O(n^2)$, ну или как минимум $O(n \ln n)$ если применять сложные максимально быстрые методы умножения, да плюс потом и это частное ещё умножать на $p$ за $O(n)$ операций ($p$ мало и влезает в регистр), а у меня в процессе оптимизации вышло что квадрат или логарифм оттуда можно убрать и сократить второе умножение, да и вообще все вычисления упростились (достаточно $2n/64$ простых умножений формата например 64*64=128). Но разумеется обратный элемент должен быть известен априори/заранее.

Вряд ли это что-то прям новое, но упоминаний не видел вроде бы.
Или никому не нужно и все пользуются более общими методами деления для получения остатков?

[worm2]

Dmitriy40, вы имеете в виду разбиение $n$ на 64-битные части $n_k$ и предвычисление $p_k = 2^{64k} \pmod p$? Тогда остаток вычисляется как $\displaystyle{\sum\limits_k\left(p_k n_k \pmod p\right)}$ и вроде бы, действительно, достигается линейная сложность.

[mathpath]

в cado-nfs нашел файл mpz_poly.cpp, и там такой код:

/* r <- a mod m */
static void
mpz_mod_barrett (mpz_ptr r, mpz_srcptr a, mpz_srcptr m, mpz_srcptr invm)
{
  size_t n = mpz_sizeinbase (m, 2);
  mpz_srcptr r_or_a = a;

  while (mpz_sizeinbase (r_or_a, 2) > n + 1)
    {
      mpz_t a1;
      size_t sr = mpz_sizeinbase (r_or_a, 2);
      mpz_init (a1);
      /* if sr <= 2n we consider the sr-n most significant bits of r,
    otherwise we take the n most significant bits */
      mpz_tdiv_q_2exp (a1, r_or_a, (sr <= 2 * n) ? n : sr - n);
      mpz_mul (a1, a1, invm);
      mpz_tdiv_q_2exp (a1, a1, n);
      mpz_mul (a1, a1, m);
      /* if sr > 2*n we have to multiply by 2^(sr-2n) */
      if (sr >= 2 * n)
   mpz_mul_2exp (a1, a1, sr - 2 * n);
      mpz_sub (r, r_or_a, a1);
      r_or_a = r;
      mpz_clear (a1);
    }
  /* now r_or_a has at most n bits */
  while (mpz_cmpabs (r_or_a, m) >= 0)
    {
      if (mpz_cmp_ui (r_or_a, 0) > 0)
   mpz_sub (r, r_or_a, m);
      else
   mpz_add (r, r_or_a, m);
      r_or_a = r;
    }
  /* now |r_or_a| < m */
  if (mpz_cmp_ui (r_or_a, 0) < 0)
    {
      mpz_add (r, r_or_a, m);
      r_or_a = r;
    }
  /* now 0 <= r_or_a < m */
  if (r_or_a == a && r != a)
    mpz_set (r, r_or_a);
}

Это реализация алгоритма Барретта для длинных чисел (многоразрядных), он не требует однократного умножения всего a на invm, а делает это по частям, итеративно уменьшая длину остатка.

пока длина(r_or_a) > n+1:
    sr = длина(r_or_a) в битах
    if sr <= 2n:
        a1 = старшие n бит r_or_a
    else:
        a1 = старшие n бит r_or_a (сдвиг на sr-n)
    
    a1 = floor(a1 * invm / 2^n)   # оценка частного
    a1 = a1 * m                    # q_est * m
    
    if sr >= 2n:
        a1 = a1 * 2^{sr-2n}        # масштабируем обратно
    
    r = r_or_a - a1
    r_or_a = r

Здесь на каждом шаге мы уменьшаем длину остатка примерно на n бит, пока он не станет длины n+1 бит, после чего делаем финальную коррекцию.
Это для сильно больших чисел, не твой случай

-- 16.04.2026, 12:20 --

Для твоего случая дипсик наваял вот такое :

(Оффтоп)

#include <stdint.h>
#include <stddef.h>
#include <string.h>

// Предвычисленные данные для конкретного модуля p
typedef struct {
    uint64_t p;           // модуль (малое число, p < 2^63 для безопасности)
    uint64_t inv_p;       // обратное к p по модулю 2^64 (p * inv_p ≡ 1 mod 2^64)
    uint64_t p_correct;   // для коррекции при знаковых операциях
} barrett_small_t;

// Вычисление обратного по модулю 2^64 (метод Ньютона)
// Требует, чтобы p было нечётным (иначе обратного не существует)
static uint64_t mod_inverse_2exp64(uint64_t p) {
    // x = p^{-1} mod 2^64
    // Начальное приближение: x0 = p (неправильно, но метод Ньютона сойдётся)
    // Лучше: x0 = p (для нечётных p) или использовать x0 = 1
    
    // Метод Ньютона для x = 1/p mod 2^64:
    // x_{k+1} = x_k * (2 - p * x_k)
    uint64_t x = 1;  // начальное приближение (для нечётных p)
    
    for (int i = 0; i < 5; i++) {  // 5 итераций достаточно для 64 бит
        x = x * (2 - p * x);
    }
    return x;
}

// Инициализация для заданного модуля p (p должно быть нечётным)
void barrett_small_init(barrett_small_t *ctx, uint64_t p) {
    ctx->p = p;
    ctx->inv_p = mod_inverse_2exp64(p);
    // Вспомогательная константа для коррекции
    ctx->p_correct = p;  // для вычитаний
}

// Основная функция редукции: r = x mod p
// x представлен массивом 64-битных слов в little-endian порядке
// Результат возвращается в r (0 <= r < p)
uint64_t barrett_small_reduce(const barrett_small_t *ctx, 
                               const uint64_t *x_words, 
                               size_t n) {
    uint64_t r = 0;  // текущий остаток (64 бита)
    
    // Проходим по словам x от старшего к младшему
    for (size_t i = n; i-- > 0; ) {
        // Формируем 128-битное число: (r << 64) | x_words[i]
        // В C нет прямой 128-битной поддержки, используем __int128 (GCC/Clang)
        unsigned __int128 t = ((unsigned __int128)r << 64) | x_words[i];
        
        // Вычисляем очередную цифру частного (64 бита)
        // q_i = (t * inv_p) mod 2^64
        uint64_t q_i = (uint64_t)(t * ctx->inv_p);
        
        // Вычисляем новый остаток: r = t - q_i * p
        // Результат помещается в 128 бит, но мы берём младшие 64 бита
        // и корректируем, если нужно
        unsigned __int128 new_r = t - (unsigned __int128)q_i * ctx->p;
        
        // new_r всегда < p * 2^64, но может быть отрицательным?
        // В беззнаковой арифметике new_r корректно, но может быть >= 2^64
        // Безопасно: берём младшие 64 бита
        r = (uint64_t)new_r;
        
        // Коррекция: если r >= p, вычитаем p (максимум 1-2 раза)
        // На самом деле, в алгоритме деления с обратным по модулю 2^64,
        // r всегда < 2p, поэтому достаточно одного вычитания
        if (r >= ctx->p) {
            r -= ctx->p;
            // Возможно, потребуется второе вычитание (редко)
            if (r >= ctx->p) {
                r -= ctx->p;
            }
        }
    }
    
    return r;
}

// Альтернативная версия, работающая со словами от младшего к старшему
// (чаще используется при потоковой обработке)
uint64_t barrett_small_reduce_stream(const barrett_small_t *ctx,
                                      const uint64_t *x_words,
                                      size_t n) {
    uint64_t r = 0;
    
    for (size_t i = 0; i < n; i++) {
        unsigned __int128 t = ((unsigned __int128)r << 64) | x_words[i];
        uint64_t q_i = (uint64_t)(t * ctx->inv_p);
        unsigned __int128 new_r = t - (unsigned __int128)q_i * ctx->p;
        r = (uint64_t)new_r;
        
        // Коррекция
        if (r >= ctx->p) {
            r -= ctx->p;
            if (r >= ctx->p) {
                r -= ctx->p;
            }
        }
    }
    
    return r;
}

Дипсик говорит, что это не баррет, а кнут - это алгоритм деления Кнута (Knuth's Algorithm D) в базисе 2^64 с предвычисленным обратным элементом, хотя время также линейное

[Dmitriy40]

worm2 в сообщении #1722485 писал(а):

Dmitriy40, вы имеете в виду разбиение $n$ на 64-битные части $n_k$ и предвычисление $p_k = 2^{64k} \pmod p$? Тогда остаток вычисляется как $\displaystyle{\sum\limits_k\left(p_k n_k \pmod p\right)}$ и вроде бы, действительно, достигается линейная сложность.

Разбиение да, остальное нет: предвычисление другое, я же сразу написал, величины $d=\lfloor 2^{64k} / p \rfloor$ (фактически $1/p$ для примерно $p<2^{64}/k$) для Баррета, индексы в сумме у $n$ и $d$ должны идти навстречу друг другу. И главное модуль брать внутри суммы у меня не нужно (это же долго! и незачем, проще было само число разделить), только два умножения (и бесплатный сдвиг на 64 бита). Частное (для последующего умножения на делитель и вычитания из исходного числа, строго по Баррету, только укороченному ибо нефиг) при этом получается некорректным, но это нивелируется дополнительным шагом коррекции остатка - всего два коротких умножения с одним условным вычитанием (ровно одна итерация Баррета для 64 битных чисел) или несколькими (примерно $\log_2 k$) условными вычитаниями.

mathpath
Ужас. Даже смотреть не буду - у меня код раз в 10 короче. К тому же я не просил написать код, а лишь известна такая модификация алгоритма или это вдруг нечто новое и важное.

-- 16.04.2026, 14:44 --

mathpath
Деление многобитных чисел порциями по $2^{64}$ с оценкой старшей цифры по Кнуту - знаю, но это совсем другое. У меня частного не получается (можно, но чуть сложнее и заметно медленнее), только остаток, и я не претендую на делители больше примерно $2^{64}/k$ (тогда эффективность теряется), меня интересуют наоборот малые делители, сильно меньше размера слова.