April 29, 2017 at 00:03

Fast Factorial Computation - PrimeSwing

    Having recently stumbled upon this article , I realized that ways of calculating factorial are rarely mentioned, other than the banal multiplication of consecutive numbers. It is necessary to correct this situation.
    I propose to consider the "asymptotically fastest" factorial algorithm!
    To begin with, I recall that the factorial n is the product of all natural numbers from 1 to n ($ n! = 1 \ cdot2 \ cdot \ ldots \ cdot n $), wherein $ 0! = 1 $;

    1. Factorial decomposition


    We introduce a function called swinging factorial as follows:

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


    This fraction will always be an integer for a simple reason - it is a multiple of the central binomial coefficient $ \ binom {n} {\ lfloor n / 2 \ rfloor} $which is equal

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


    Expanding the definition of swinging factorial , we get a new recurrent factorial formula:

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


    It will be especially good if we learn how to efficiently calculate values $ n \ wr $.

    2. Simple factors swinging factorial


    We denote $ l_p (n \ wr) $ as a prime power $ p $ in primary decomposition $ n \ wr $. Then the following formula will be valid:

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


    Evidence
    We use the Legendre theorem on prime factors of the factorial :

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


    For the last expression, we use the fact that $ j-2 \ lfloor j / 2 \ rfloor = j \: mod \: 2 $, and we get the formula we need.

    Consequently, $ l_p (n \ wr) \ leqslant log_p (n) $ and $ p ^ {l_p (n \ wr)} \ leqslant n $. If a$ p $ odd then $ l_p (p ^ a \ wr) = a $. Other special cases:

    $$ 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) $ here means the number of units in binary representation of a number $ n $. All these facts can be used for additional optimization in the code. I will not give evidence, if you wish, you can easily get it yourself.

    Now, knowing the degrees of all prime divisors$ n \ wr $, we have a way to calculate swinging factorial :

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


    3. The complexity of the algorithm


    It can be shown that the calculation $ n \ wr $ has difficulty $ O (n (log \: n) ^ 2log \: log \: n) $. Oddly enough, computing$ n! $it has the same complexity ( the Schoenhage-Strassen algorithm is used in the assessment , hence such an interesting laboriousness; evidence by reference at the end of the article).

    Despite the fact that formally multiplying numbers from 1 to n has the same complexity, the PrimeSwing algorithm in practice is the fastest.

    UPDATE: as was noted in this comment , here I was mistaken, multiplying numbers from 1 to n has a lot of complexity.

    References and implementation



    Java implementation
    // 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;
}

    Tags:

    Also popular now: