Быстрый способ вычисления n! mod m, где m - простое число?

Ответ 1

Расширение моего комментария к ответу:

Да, есть более эффективные способы сделать это. Но они чрезвычайно беспорядочны.

Поэтому, если вам действительно не нужна эта дополнительная производительность, я не предлагаю попробовать их реализовать.

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

• Разделение на инвариантные целые числа с использованием умножения
• Сокращение Montgomery

Эти методы бывают быстрыми, потому что они существенно исключают модуль.

Только эти методы должны дать вам умеренное ускорение. Чтобы быть действительно эффективным, вам может понадобиться развернуть цикл, чтобы обеспечить лучший IPC:

Что-то вроде этого:

ans0 = 1
ans1 = 1
for i in range(1,(n+1) / 2):
    ans0 = ans0 * (2*i + 0) % modulus    
    ans1 = ans1 * (2*i + 1) % modulus    

return ans0 * ans1 % modulus

но принимая во внимание нечетное число итераций и объединяя его с одним из методов, с которыми я связан выше.

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

Обратите внимание, что хотя мой ответ является языковым агностиком, он предназначен в основном для C или С++.

Ответ 2

n может быть сколь угодно большим

Ну, n не может быть сколь угодно большим - если n >= m, то n! ≡ 0 (mod m) (поскольку m является одним из факторов, по определению факторного).

Предполагая n << m, и вам нужно точное значение, ваш алгоритм не может быть быстрее, насколько мне известно. Однако, если n > m/2, вы можете использовать следующее тождество (теорема Вильсона - Спасибо @Даниэль Фишер!)

чтобы ограничить количество умножений примерно на m-n

(m-1)! ≡ -1 (mod m)
1 * 2 * 3 * ... * (n-1) * n * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m)
n! * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m)
n! ≡ -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m)

Это дает нам простой способ вычисления n! (mod m) в m-n-1 умножениях, а также модульный обратный:

Мы можем перефразировать вышеупомянутое уравнение по-другому, что может или не может выполняться несколько быстрее. Используя следующее тождество:

мы можем перефразировать уравнение как

n! ≡ -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m)
n! ≡ -[(n+1-m) * ... * (m-2-m) * (m-1-m)]-1 (mod m)
       (reverse order of terms)
n! ≡ -[(-1) * (-2) * ... * -(m-n-2) * -(m-n-1)]-1 (mod m)
n! ≡ -[(1) * (2) * ... * (m-n-2) * (m-n-1) * (-1)(m-n-1)]-1 (mod m)
n! ≡ [(m-n-1)!]-1 * (-1)(m-n) (mod m)

Это может быть написано на Python следующим образом:

Если вам не нужно точное значение, жизнь становится немного проще - вы можете использовать приближение Стирлинга для вычисления приблизительного значения в O(log n) time (используя возведение в степень по квадрату).

Наконец, я должен упомянуть, что если это критически важно и вы используете Python, попробуйте перейти на С++. Из личного опыта вы должны ожидать примерно на порядок увеличения скорости и т.д., Просто потому, что это именно то, что связано с жестким циклом, связанным с процессором, который отличается от исходного кода (по какой-либо причине GMP кажется гораздо более тонко настроенный, чем Python Bignum).

Ответ 3

п! mod m можно вычислить в операциях O (n ^{1/2 + epsilon;}) вместо наивного O (n). Это требует использования многочленного умножения FFT и имеет смысл только для очень большого n, например. n > 10 ⁴.

Схема алгоритма и некоторые тайминги можно увидеть здесь: http://fredrikj.net/blog/2012/03/factorials-mod-n-and-wilsons-theorem/

Ответ 4

Если мы хотим вычислить M = a*(a+1) * ... * (b-1) * b (mod p), мы можем использовать следующий подход, если предположить, что мы можем добавить, вычесть и умножить быстро (mod p) и получить сложность времени выполнения O( sqrt(b-a) * polylog(b-a) ).

Для простоты предположим, что (b-a+1) = k^2, является квадратом. Теперь мы можем разделить наш продукт на k частей, т.е. M = [a*..*(a+k-1)] *...* [(b-k+1)*..*b]. Каждый из факторов этого продукта имеет вид p(x)=x*..*(x+k-1), для соответствующего x.

Используя алгоритм быстрого умножения многочленов, такой как алгоритм Schönhage-Strassen, можно разделить и покорить, можно найти коэффициенты многочлена p(x) in O( k * polylog(k) ). Теперь, по-видимому, существует алгоритм подстановки точек k в том же полиноме степени-k в O( k * polylog(k) ), что означает, что мы можем быстро вычислить p(a), p(a+k), ..., p(b-k+1).

Этот алгоритм подстановки многих точек в один многочлен описан в книге "Prime numbers" С. Померанса и Р. Крэндалла. В конце концов, когда у вас есть эти значения k, их можно умножить на O(k) и получить желаемое значение.

Обратите внимание, что все наши операции, в которых выполняется (mod p). Точное время работы O(sqrt(b-a) * log(b-a)^2 * log(log(b-a))).

Ответ 5

Расширяя мой комментарий, это занимает около 50% времени для всех n в [100, 100007], где m = (117 | 1117):

Function facmod(n As Integer, m As Integer) As Integer
    Dim f As Integer = 1
    For i As Integer = 2 To n
        f = f * i
        If f > m Then
            f = f Mod m
        End If
    Next
    Return f
End Function

Ответ 6

Если n = (m - 1) для простого m, то http://en.wikipedia.org/wiki/Wilson 's_theorem n! mod m = (m - 1)

Также, как уже указывалось, n! mod m = 0, если n > m

Ответ 7

Я нашел эту следующую функцию на quora:
При f (n, m) = n! mod m;

function f(n,m:int64):int64;
         begin
              if n = 1 then f:= 1
              else f:= ((n mod m)*(f(n-1,m) mod m)) mod m;
         end;

Вероятно, бить с использованием цикла времени и умножать большое количество, хранящееся в строке. Кроме того, он применим к любому целому числу m.
Ссылка, в которой я нашел эту функцию: https://www.quora.com/How-do-you-calculate-n-mod-m-where-n-is-in-the-1000s-and-m-is-a-very-large-prime-number-eg-n-1000-m-10-9+7

Ответ 8

Этот алгоритм имеет O(n log log(n)) временную сложность предварительной обработки (из-за сита) и сложность пространства O (r (n)), где r(n) - prime счетная функция. После предварительной обработки запросов x! где x <= N - O(r(x) * log(x)).

К сожалению, это становится неосуществимым при 10 ¹⁰ так как r (1e9) - 455 052 511, для чего требуется почти 2 ГБ памяти для хранения простых целых чисел. Однако для 10 ⁹ нам понадобилось всего около 300 МБ памяти.

Предположим, что mod равен 1e9 + 7, чтобы избежать переполнения в С++, но это может быть что угодно.

Способ вычисления N! mod m fast

Чтобы вычислить N! быстро, мы можем разложить N на его простые множители. Можно заметить, что если р - простой множитель N!, то р должно быть &lt = N, так как N! является произведением целых чисел от 1 до N. Учитывая эти факты, мы можем вычислить N! используя только простые числа от 2 до N (включительно). Согласно Wikipedia, существует 50 847 534 грубо (~ 5 * 10 ^ 7) простых чисел, меньших или равных 10 ⁹. Это будет формула для использования:

v _{p (N!)} - это наибольшее число, такое, что p ^{v _p (N!)} делит N!. Формула Лежандра вычисляет v _p (N!) в O (log (N)) времени. Вот формула из Википедии.

int factorial(int n) {
    int ans = 1;
    for (int i = 0; i < primes.size() && primes[i] <= n; i++) {
        int e = 0;
        for (int j = n; j > 0; ) {
            j /= primes[i];
            e += j;
        }
        ans = ((u64) ans * fast_power(primes[i], e)) % mod;
    }
    return ans;
}

fast_power(x, y) возвращает x ^y% mod в O (log (y)) время, используя экспонирование путем возведения в квадрат.

Генерация простых чисел

Я генерирую простые числа, используя Сито Эратосфена. В моей реализации сита я использую битрейт для сохранения памяти. Моя реализация заняла 16-18 секунд, чтобы генерировать простые числа до 10 ^ 9, когда я скомпилирован с флагом -O3 с процессором второго поколения i3. Я уверен, что существуют лучшие алгоритмы/реализации для генерации простых чисел.

const int LIMIT = 1e9;
bitset<LIMIT+1> sieve;
vector<int> primes;

void go_sieve() {
    primes.push_back(2);
    int i = 3;
    for (int i = 3; i * i <= LIMIT; i+=2) {
        if (!sieve.test(i)) {
            primes.push_back(i);
            for (int j = i * i; j <= LIMIT; j += i) {
                sieve.set(j);
            }
        }
    }
    while (i <= LIMIT) {
        if (!sieve.test(i)) primes.push_back(i);
        i+=2;
    }
}

Потребление памяти

Согласно Wikipedia, существует 50 847 534 простых числа, которые меньше или равны 10 ^ 9. Так как 32-битовое целое число составляет 4 байта, основной вектор требует 203,39 МБ (50,847,534 * 4 байта). Битовый набор сит требует (125 МБ) 10 ^ 9 бит.

Производительность

Я смог вычислить 1,000,000,000! через 1,17 секунды после обработки.

Примечание. Я использовал алгоритм с флагом -O3.

Ответ 9

Предполагая, что оператор "mod" на выбранной вами платформе достаточно быстр, вы ограничены в первую очередь скоростью, с которой вы можете вычислить n! и пространство, в котором вы можете его вычислить.

Тогда это по существу двухэтапная операция:

• Вычислить n! (есть много быстрых алгоритмов, поэтому я не буду повторять здесь)
• Возьмем mod результата

Нет необходимости усложнять вещи, особенно если скорость является критическим компонентом. В общем, выполняйте как можно меньше операций внутри цикла.

Если вам нужно многократно рассчитать n! mod m, вы можете захотеть сохранить значения, исходящие из функции, выполняющей вычисления. Как всегда, это классическая комбинация пространства/времени, но таблицы поиска очень быстрые.

Наконец, вы можете комбинировать memoization с рекурсией (и, если необходимо, с батутами), чтобы быстро получить вещи.