ibessonov Apr 29 2017 at 00:03

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

3 min

16K

Наткнувшись недавно на эту статью, я понял, что редко упоминаются способы вычисления факториала, отличные от банального перемножения последовательных чисел. Нужно эту ситуацию исправить.
Предлагаю рассмотреть «асимптотически наиболее быстрый» алгоритм вычисления факториала!
Для начала напомню, что факториал n — это произведение всех натуральных чисел от 1 до n ( $n!=1\cdot2\cdot\ldots\cdot n$ ), при этом $inline$ ;

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)$ как степень простого числа $inline$ в примарном разложении $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$ . Если $inline$ нечётно, то $l_p(p^a\wr)=a$ . Другие частные случаи:

$\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}$

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

Теперь, зная степени всех простых делителей $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)$ . Как ни странно, вычисление $inline$ имеет ту же сложность (в оценке используется алгоритм Шёнхаге-Штрассена, отсюда и такая интересная трудоёмкость; доказательства по ссылке в конце статьи).

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

UPDATE: как было замечено в этом комментарии, тут я ошибся, перемножение чисел от 1 до n имеет большую трудоёмкость.

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

Реализация на 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;
}

Tags:

Hubs: