Главная

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

2017-04-28 в 21:03, admin, рубрики: primeswing, Алгоритмы, математика, простые числа, факториал, факторизация, метки: primeswing, факториал

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

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

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

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

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

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

$binom{n}{lfloor n/2rfloor}=frac{n!}{lfloor n/2rfloor!cdotlceil n/2rceil!}$

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

$n!=lfloor n/2rfloor!^2cdot nwr$

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

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

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

$l_p(nwr)=sum_{kgeqslant1}lfloorfrac{n}{p^k}rfloor:mod:2$

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

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

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

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

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

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

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

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

$nwr=prod_{pleqslant n}p^{l_p(nwr)}$

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

Можно показать, что вычисление $nwr$ имеет сложность $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;
}

Автор: ibessonov

Источник