28.04.2017 21:03
ibessonov 21 6900 Источник
Наткнувшись недавно на эту статью, я понял, что редко упоминаются способы вычисления факториала, отличные от банального перемножения последовательных чисел. Нужно эту ситуацию исправить.
Предлагаю рассмотреть «асимптотически наиболее быстрый» алгоритм вычисления факториала!
Для начала напомню, что факториал 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 на практике оказывается самым быстрым.

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;
}

Поделиться с друзьями
-->

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


  1. koldyr
    29.04.2017 07:32
    #10197154

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


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


    1. ibessonov
      29.04.2017 09:20
      #10197194

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


      1. koldyr
        29.04.2017 11:30
        #10197256

        Асимптотику распределения простых чисел я знаю.
        Я правильно понимаю, что придется хранить или генерировать таблицу простых чисел от 2 до n/2?
        Какова трудоемкость детерминированного алгоритма теста на простоту?


        1. ibessonov
          29.04.2017 11:39
          #10197266

          Не совсем, понадобится таблица простых чисел от 2 до n.
          Решето Эратосфена составляется за O(n * log log n), это меньше, чем время вычисления факториала. Отдельного способа проверки на простоту алгоритм не требует.


          1. koldyr
            29.04.2017 11:49
            #10197278

            Прекрасно. Как-то я упустил, что сложность решета ниже сложности вычисления факториала.


  1. choupa
    29.04.2017 11:31
    #10197258
    +1

    Мне кажется, самый практичный метод вычисления факториала

    array[0, 1, 2, 6, 24, 120, ...]

    Простите за серость, а зачем нужно вычислять большие факториалы? Криптография? Комбинаторика? Ряды Тейлора? Собеседование в Гугл?


    1. ibessonov
      29.04.2017 11:35
      #10197262

      Боюсь, что это вопрос не ко мне. Вероятнее всего различные алгоритмы вычисления факториала имеют соревновательный характер.
      Хранить массив может быть очень затратно по памяти, особенно если нужны значения порядка 1000000! (вдруг они кому-то нужны). На практике я такого не встречал и сам всегда хранил кэшированные значения в обычном массиве (во время решения задач по программированию).


    1. alexeykuzmin0
      03.05.2017 11:14
      #10201100

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


  1. mikhaelkh
    29.04.2017 11:46
    #10197274
    +1

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

    Предположим, что разделим числа на две половины, рекурсивно вычислим произведение и перемножим половины. Сложность у этого метода ~M(n*log(n))*log(n), что хуже ~M(n*log(n)) у PrimeSwing. Или как предлагается перемножать?


    1. ibessonov
      29.04.2017 12:12
      #10197296

      Извиняюсь, действительно у меня в этом моменте ошибка. Вычисление с помощью PrimeSwing имеет трудоёмкость как у простого перемножение двух чисел порядка n!, а не вычисления n! стандартным способом. Сейчас внесу обновление в статью.
      Спасибо!


  1. koldyr
    29.04.2017 11:53
    #10197288

    Стоит отметить, что статья имеет практическую направленность, но согласно вики:
    Метод Шёнхаге — Штрассена считался асимптотически быстрейшим методом с 1971 до 2007 годы, пока не был заявлен новый метод с лучшей оценкой сложности умножения.[3] На практике метод Шёнхаге — Штрассена начинает превосходить более ранние классические методы, такие как умножение Карацубы и алгоритм Тоома — Кука (обобщение метода Карацубы), начиная с целых чисел порядка 10^10000-10^40000.
    Таким образом становится актуальным вопрос о константе в оценке.


    1. ibessonov
      29.04.2017 12:23
      #10197300

      Интересно вышло, алгоритм PrimeSwing был опубликован в 2008-м, автор уже мог знать о более оптимальном способе умножения чисел.
      Если я верно понял, автор производил вычисления факториалов примерно до 1000000, то есть результаты были не настолько большими, чтобы использовать Алгоритма Фюрера (если верить википедии). А так да, теоретическую оценку наверняка можно улучшить.
      На практике я не заморачивался и использовал то умножение, какое было реализовано за меня в BigInteger.
      Спасибо за замечание.


      1. alexeykuzmin0
        03.05.2017 11:16
        #10201106

        В BigInteger (по крайней мере раньше было) реализовано умножение столбиком (за квадрат). И не забудьте еще расходы времени на перевод из десятичной в двоичную СС.


        1. vlanko
          04.05.2017 11:51
          #10202970

          В Java используются более быстрые алгоритмы. А перевод для вывода — если печатать в консоль — да, займет время :)


    1. mikhaelkh
      29.04.2017 12:41
      #10197320

      Хорошо про умножение с практической точки зрения написано здесь


  1. malishich
    29.04.2017 16:39
    #10197458

    Стать интересная, не знал что и здесь простые числа применяются. На практике мне не встречались задачи где нужно было вычислять просто факториалы чисел больше 20-30 (например в uint64_t влезает 20, а вот 21! уже нет, зато в double влезает, хоть и с погрешностью). Встречались отношения (дроби) больших факториалов (комбинаторные сочетания Cmn для близких чисел и размещения Amn для сильно разных чисел), порядка 300-т. Но в этих задачах я решал просто разложением на простые множители всех чисел произведения в числителе и знаменателе, подсчёту количества таких множителей в числителе и знаменателе, сокращению, и дальнейшему вычислению отношения (там обычно много чего сокращалось и части дроби «влазили» с малой погрешностью в тип double, но чаще вообще в uint64_t).


    1. MikailBag
      29.04.2017 22:08
      #10197754

      Зачем считать для цешек факториалы, если есть треугольник паскаля?
      Я могу оценить время на генерацию N строк как O(N^3).
      В Вашем случае (N = 300) — пара минут максимум.


      1. alexeykuzmin0
        03.05.2017 11:19
        #10201112

        Зачем использовать треугольник Паскаля, если есть замечательное соотношение C(n, k) = C(n — 1, k — 1) * (n / k)? Вычисляется, в отличие от треугольника Паскаля, за O(k) операций над числами, а из операций потребуется только умножение и деление на короткое число (причем деление всегда гарантированно нацело будет).


        1. MikailBag
          03.05.2017 16:28
          #10201745

          Ух ты ж.
          Спасибо, отличное соотношение!


  1. mikhaelkh
    29.04.2017 21:41
    #10197738

    В этом алгоритме мы по сути выделяем большой множитель, являющийся квадратом. То есть A[k] = A[k+1]^2 * B[k+1], A[0] = n!, для вычисления A[k] мы вычисляем A[k+1], B[k+1]. Но мы можем выделить ещё больший множитель, оставляя в B[k+1] только произведение простых, входящие в A[k] нечётных степенях. Кажется, это должно работать ещё быстрее, хотя разница должна быть невелика.


  1. nikitadanilov
    30.04.2017 02:26
    #10197888

    «в примарном разложении» — должно быть «в разложении на простые множители». «Примарное разложение» это немного другое (см. теорему Ласкера).