Можно ли использовать умножение Монтгомери для ускорения вычислений (большое число)! % (некоторые простые)

Question

Этот вопрос возникает из комментария, который я почти написал ниже этого вопроса, где Зак вычисляет факториал большого числа по модулю большого числа (которое мы будем считать основным для этого вопроса). Зак использует традиционное вычисление факториала, беря остаток при каждом умножении.

Я почти прокомментировал, что альтернативой для рассмотрения было умножение Монтгомери, но если больше думать об этом, я видел только эту технику, используемую для ускорения нескольких умножений на один и тот же множитель (в частности, для ускорения вычисленияⁿ mod p).

Мой вопрос: можно ли использовать умножение Монтгомери для ускорения вычисления n! мод р для больших п и р?

Answer 1

Наивно нет; вам нужно преобразовать каждое из n слагаемых произведения в "пространство Монтгомери", чтобы у вас было n полных сокращений mod m, аналогично "обычному" алгоритму.

Однако факториал - это не просто произвольное произведение n членов; это намного более структурировано. В частности, если у вас уже есть "Montgomerized" kr mod m, то вы можете использовать очень дешевое сокращение, чтобы получить (k+1)r mod m,

Так что это вполне осуществимо, хотя я не видел, чтобы это было сделано раньше. Я пошел дальше и написал быструю и грязную реализацию (очень непроверенную, я бы совсем не доверял ей):

// returns m^-1 mod 2**64 via clever 2-adic arithmetic (http://arxiv.org/pdf/1209.6626.pdf)
uint64_t inverse(uint64_t m) {
    assert(m % 2 == 1);
    uint64_t minv = 2 - m;
    uint64_t m_1 = m - 1;
    for (int i=1; i<6; i+=1) { m_1 *= m_1; minv *= (1 + m_1); }
    return minv;
}

uint64_t montgomery_reduce(__uint128_t x, uint64_t minv, uint64_t m) {
    return x + (__uint128_t)((uint64_t)x*-minv)*m >> 64;
}

uint64_t montgomery_multiply(uint64_t x, uint64_t y, uint64_t minv, uint64_t m) {
    return montgomery_reduce(full_product(x, y), minv, m);
}

uint64_t montgomery_factorial(uint64_t x, uint64_t m) {
    assert(x < m && m % 2 == 1);
    uint64_t minv = inverse(m); // m^-1 mod 2**64
    uint64_t r_mod_m = -m % m;  // 2**64 mod m
    uint64_t mont_term = r_mod_m;
    uint64_t mont_result = r_mod_m;
    for (uint64_t k=2; k<=x; k++) {
        // Compute the montgomerized product term: kr mod m = (k-1)r + r mod m.
        mont_term += r_mod_m;
        if (mont_term >= m) mont_term -= m;
        // Update the result by multiplying in the new term.
        mont_result = montgomery_multiply(mont_result, mont_term, minv, m);
    }
    // Final reduction
    return montgomery_reduce(mont_result, minv, m);
}

и сравнил его с обычной реализацией:

__uint128_t full_product(uint64_t x, uint64_t y) {
    return (__uint128_t)x*y;
}

uint64_t naive_factorial(uint64_t x, uint64_t m) {
    assert(x < m);
    uint64_t result = x ? x : 1;
    while (x --> 2) result = full_product(result,x) % m;
    return result;
}

и против обычной реализации с некоторым встроенным ассемблером, чтобы исправить незначительную неэффективность:

uint64_t x86_asm_factorial(uint64_t x, uint64_t m) {
    assert(x < m);
    uint64_t result = x ? x : 1;
    while (x --> 2) {
        __asm__("mov %[result], %%rax; mul %[x]; div %[m]"
                : [result] "+d" (result) : [x] "r" (x), [m] "r" (m) : "%rax", "flags");
    }
    return result;
}

На моем ноутбуке Haswell были получены следующие результаты:

implementation   speedup
---------------------------
naive            1.00x
x86_asm          1.76x
montgomery       5.68x

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

Это интересный подход для "скромных" х и м. Как только x станет большим, различные подходы, которые имеют сублинейную сложность в x, обязательно победят; Факториал имеет так много структуры, что этот метод не использует в своих интересах.