SavePearlHarbor

Быстрое вычисление факториала — PrimeSwing

—

от автора

admin

Наткнувшись недавно на эту статью, я понял, что редко упоминаются способы вычисления факториала, отличные от банального перемножения последовательных чисел. Нужно эту ситуацию исправить.
Предлагаю рассмотреть «асимптотически наиболее быстрый» алгоритм вычисления факториала!

Для начала напомню, что факториал n — это произведение всех натуральных чисел от 1 до n ($n!=1\cdot2\cdot\ldots\cdot n$), при этом $0!=1$;

1. Декомпозиция факториала

Введём функцию, именуемую swinging factorial, следующим образом:

$n\wr=\frac{n!}{\lfloor n/2\rfloor!^2}$

Данная дробь всегда будет целым числом по простой причине — она является делителем центрального биномиального коэффициента $\binom{n}{\lfloor n/2\rfloor}$, равного

$\binom{n}{\lfloor n/2\rfloor}=\frac{n!}{\lfloor n/2\rfloor!\cdot\lceil n/2\rceil!}$

Разворачивая определение swinging factorial, мы получим новую рекуррентную формулу факториала:

$n!=\lfloor n/2\rfloor!^2\cdot n\wr$

Она будет особенно хороша, если мы научимся эффективно вычислять значения $n\wr$.

2. Простые множители swinging factorial

Обозначим $l_p(n\wr)$ как степень простого числа $p$ в примарном разложении $n\wr$. Тогда будет справедлива следующая формула:

$l_p(n\wr)=\sum_{k\geqslant1}\lfloor\frac{n}{p^k}\rfloor\:mod\:2$

Доказательство

Воспользуемся теоремой Лежандра о простых множителях факториала:

$\begin{array}{ccl} l_p(n!/\lfloor n/2\rfloor!^2)&=&l_p(n!)-2l_p(\lfloor n/2\rfloor!)\\ &=&\sum_{k\geqslant1}\lfloor n/p^k\rfloor-2\sum_{k\geqslant1}\lfloor \lfloor n/2\rfloor/p^k\rfloor\\ &=&\sum_{k\geqslant1}(\lfloor n/p^k\rfloor-2\lfloor \lfloor n/p^k\rfloor/2\rfloor) \\\end{array}$

Для последнего выражения воспользуемся тем фактом, что $j-2\lfloor j/2\rfloor=j\:mod\:2$, и получим нужную нам формулу.

Как следствие, $l_p(n\wr)\leqslant log_p(n)$ и $p^{l_p(n\wr)}\leqslant n$. Если $p$ нечётно, то $l_p(p^a\wr)=a$. Другие частные случаи:

$$display$$\begin{array}{lrrl} (a)&\lfloor n/2\rfloor &< p \leqslant & n & \Rightarrow & l_p(n\wr)=1\\ (b)&\lfloor n/3\rfloor &< p \leqslant & \lfloor n/2\rfloor & \Rightarrow & l_p(n\wr)=0\\ (c)&\sqrt{n} &< p \leqslant & \lfloor n/3\rfloor & \Rightarrow & l_p(n\wr)=\lfloor n/p\rfloor\:mod\:2\\ (d)&2 &< p \leqslant & \sqrt{n} & \Rightarrow & l_p(n\wr) < log_2(n)\\ (e)& & p = & 2 & \Rightarrow & l_p(n\wr) =\sigma_2(\lfloor n/2\rfloor)\\ \end{array}$$display$$

$\sigma_2(n)$ здесь означает количество единиц в двоичном представлении числа $n$. Все эти факты могут быть использованы для дополнительной оптимизации в коде. Доказательства я приводить не буду, при желании вы легко сможете получить их самостоятельно.

Теперь, зная степени всех простых делителей $n\wr$, у нас есть способ вычисления swinging factorial:

$n\wr=\prod_{p\leqslant n}p^{l_p(n\wr)}$

3. Трудоёмкость алгоритма

Можно показать, что вычисление $n\wr$ имеет сложность $O(n(log\:n)^2log\:log\:n)$. Как ни странно, вычисление $n! $ имеет ту же сложность (в оценке используется алгоритм Шёнхаге-Штрассена, отсюда и такая интересная трудоёмкость; доказательства по ссылке в конце статьи).

Несмотря на то, что формально перемножение чисел от 1 до n имеет ту же трудоёмкость, алгоритм PrimeSwing на практике оказывается самым быстрым.

Ссылки и реализация

Реализация на Java

// main function public static BigInteger factorial(int n) {     return factorial(n, primes(n)); }  // recursive function with shared primes array private static BigInteger factorial(int n, int[] primes) {     if (n < 2) return BigInteger.ONE;     BigInteger f = factorial(n / 2, primes);     BigInteger ps = primeSwing(n, primes);     return f.multiply(f).multiply(ps); }  // swinging factorial function private static BigInteger primeSwing(int n, int[] primes) {     List<BigInteger> multipliers = new ArrayList<>();     for (int i = 0; i < primes.length && primes[i] <= n; i++) {         int prime = primes[i];         BigInteger bigPrime = BigInteger.valueOf(prime);         BigInteger p = BigInteger.ONE;         int q = n;         while (q != 0) {             q = q / prime;             if (q % 2 == 1) {                 p = p.multiply(bigPrime);             }         }         if (!p.equals(BigInteger.ONE)) {             multipliers.add(p);         }     }     return product(multipliers, 0, multipliers.size() - 1); }  // fast product for the list of numbers private static BigInteger product(List<BigInteger> multipliers, int i, int j) {     if (i > j) return BigInteger.ONE;     if (i == j) return multipliers.get(i);     int k = (i + j) >>> 1;     return product(multipliers, i, k).multiply(product(multipliers, k + 1, j)); }  // Eratosthenes sieve private static int[] primes(int upTo) {     upTo++;     if (upTo >= 0 && upTo < 3) {         return new int[]{};     }      int length = upTo >>> 1;     boolean sieve_bool[] = new boolean[length];     for (int i = 1, iterations = (int) Math.sqrt(length - 1); i < iterations; i++) {         if (!sieve_bool[i]) {             for (int step = 2 * i + 1, j = i * (step + 1); j < length; j += step) {                 sieve_bool[j] = true;             }         }     }      int not_primes = 0;     for (boolean not_prime : sieve_bool) {         if (not_prime) not_primes++;     }      int sieve_int[] = new int[length - not_primes];     sieve_int[0] = 2;     for (int i = 1, j = 1; i < length; i++) {         if (!sieve_bool[i]) {             sieve_int[j++] = 2 * i + 1;         }     }      return sieve_int; }

ссылка на оригинал статьи https://habrahabr.ru/post/323770/

Комментарии

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *