Быстрое вычисление факториала — PrimeSwing29.04.2017 10:48

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

Для начала напомню, что факториал 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 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 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;
}

Комментарии (2)

  • 29 апреля 2017 в 07:32

    0

    ↑

    ↓

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

    С какой сложностью нам достается это знание?

    • 29 апреля 2017 в 09:20

      0

      ↑

      ↓

      Простых чисел от 1 до n всего O (n / log n), вычисление lp происходит за O (log n).
      Сложностью же перемножения простых чисел в соответствующих степенях будет O (n (log n)2 log log n) — это уже не очевидная оценка, её доказательство есть по второй ссылке.

© Habrahabr.ru