Factoring Factorials

January 24, 2014

The trick is to build up the solution from the primes less than n instead of breaking down the solution starting from n!. Our solution is due to Will Ness. Consider 33!:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
2   2   2   2    2     2     2     2     2     2     2     2     2     2     2     2
    2       2          2           2           2           2           2           2
            2                      2                       2                       2
                                   2                                               2
                                                                                   2
  3     3     3        3        3        3        3        3        3        3        3
              3                          3                          3
                                                                    3
      5          5              5              5              5              5
                                                              5
          7                  7                    7                    7
                   11                               11                               11
                         13                                     13
                                     17    19          23                29    31

The pattern is obvious: For primes less than the square root of n, compute the number of times each prime and each prime power appears. For the remaining primes less than half of n, compute the number of times each prime appears (but not each prime power). Then each prime that still remains appears once. Like this:

33! = 2^( 33 div 2 + 33 div 4 + 33 div 8 + 33 div 16 + 33 div 32) *
      3^( 33 div 3 + 33 div 9 + 33 div 27) *
      5^( 33 div 5 + 33 div 25) *
      ----
      7^( 33 div 7) * 11^( 33 div 11) * 13^( 33 div 13) *
      ----
      17 * 19 * 23 * 29 * 31

It’s easy to reduce that calculation to Scheme:

(define (fact-fact n) ; prime factors of n factorial
  (let loop ((ps (primes n)) (fs (list)))
    (cond ((null? ps) (reverse fs))
          ((< (* (car ps) (car ps)) n)
            (let ((p (car ps)))
              (let ((k (let loop ((q p) (k 0))
                         (if (< n q) k
                           (loop (* q p) (+ k (quotient n q)))))))
                (loop (cdr ps) (cons (cons (car ps) k) fs)))))
          ((< (+ (car ps) (car ps)) n)
            (loop (cdr ps) (cons (cons (car ps) (quotient n (car ps))) fs)))
          (else (loop (cdr ps) (cons (cons (car ps) 1) fs))))))

For example:

> (fact-fact 33)
((2 . 31) (3 . 15) (5 . 7) (7 . 4) (11 . 3) (13 . 2) (17 . 1)
(19 . 1) (23 . 1) (29 . 1) (31 . 1))

You can run the program at http://programmingpraxis.codepad.org/nFdoreYl, where there is also a demonstration that the calculation of the factors of 33! is correct.

About these ads

Pages: 1 2

8 Responses to “Factoring Factorials”

  1. #include
    #include

    // reducing numbers from biggest to 2
    // 16 -> 2*8 … 8 -> 2*4 … 4 -> 2*2 …

    void main(int argc, char **argv) {
    int n = atoi(argv[1]);
    int *p = (int *) malloc(sizeof(int) * (n + 1));
    int i, j, d;
    for(i = 0; i 1; i–)
    if(p[i]) {
    for(j = i + i, d = 2; j <= n; j += i, d++) {
    if(p[j]) {
    p[i] += p[j];
    p[d] += p[j];
    p[j] = 0;
    }
    }
    }
    printf("1");
    for(i = 2; i <= n; i++)
    if(p[i])
    printf(" * %i^%i", i, p[i]);
    printf("\n");
    }

  2. // (sorry...)
    #include <stdio.h>
    #include <stdlib.h>
    
    // reducing numbers from biggest to 2
    // 16 -> 2*8 ... 8 -> 2*4 ... 4 -> 2*2 ... 
    
    void main(int argc, char **argv) {
      int n = atoi(argv[1]);
      int *p = (int *) malloc(sizeof(int) * (n + 1));
      int i, j, d;
      for(i = 0; i <= n; i++)
        p[i] = 1;
      for(i = n; i > 1; i--)
        if(p[i]) {
          for(j = i + i, d = 2; j <= n; j += i, d++) {
            if(p[j]) {
              p[i] += p[j];
              p[d] += p[j];
              p[j] = 0;
            }
          }
        }
      printf("1");
      for(i = 2; i <= n; i++)
        if(p[i])
          printf(" * %i^%i", i, p[i]);
      printf("\n");
    }
    
  3. Recursivelly adding factorizations

    import Data.List
    import Data.Monoid
    import Data.Numbers.Primes

    facfac 1 = Factors [(1, 1)]
    facfac n = factorization n `mappend` facfac (n – 1)

    – Factorizations Monoid

    data Factorization a = Factors [(a, a)] deriving Show

    factorization = Factors . map (\ps -> (head ps, length ps)) . group . primeFactors

    instance Integral a => Monoid (Factorization a) where
    mempty = Factors []
    mappend (Factors []) fs = fs
    mappend fs (Factors []) = fs
    mappend xfs@(Factors (xf@(x, a):xs)) yfs@(Factors (yf@(y, b):ys)) =
    case x `compare` y of
    EQ -> Factors ((x, a + b):ms)
    LT -> Factors (xf:ns)
    GT -> Factors (yf:os)
    where Factors ms = mappend (Factors xs) (Factors ys)
    Factors ns = mappend (Factors xs) yfs
    Factors os = mappend xfs (Factors ys)

  4. (Excuse me, I can’t modify previous post)
    In C version, we can reduce main loop from “i = n” to “i = n >> 1″.
    That version has O(n log log n) and are not needed primes calculation nor mul nor div operations (space is O(n)).

  5. David said

    Did in Go using my Project Euler library. Since it’s designed to solve PE problems everything is done with type uint64.

    package main
    
    import (
        . "fmt"
        . "github.com/dgottner/euler"
        "os"
        "strconv"
    )
    
    func main() {
        var n, p uint64
    
        if len(os.Args) < 2 {
            Println("Usage: factfact n")
            return
        } else {
            val, err := strconv.ParseInt(os.Args[1], 10, 64)
            if err != nil {
                Println(err.Error())
                return
            }
            n = uint64(val)
        }
    
        Printf("%d! = 1", n)
        for p = 2; p <= n; p = NextPrime(p) {
            exp := n / p
            for fac := p * p; fac <= n; fac *= p {
                exp += n / fac
            }
    
            Printf(" * %d", p)
            if exp != 1 {
                Printf("^%d", exp)
            }
        }
    
        Println()
    }
    
  6. kernelbob said

    That’s a fun little problem. Here is my solution. It is quick and dirty, not efficient.

    from itertools import count, takewhile
    
    
    def primes():
        yield 2
        primes = [2]
        n = 3
        while True:
            if all(n % p for p in primes):
                yield n
                primes.append(n)
            n += 2
    
    
    def pff(n):
        """prime factorization of n!
           returns a list of pairs (prime, power).
        """
        
        up_to_n = lambda iter: takewhile(lambda p: p <= n, iter)
        return [(prime,
                 sum(n/w
                     for w in up_to_n(prime**power
                                      for power in count(1))))
                for prime in up_to_n(primes())]
                
    
    from math import factorial
    from operator import mul
    print factorial(100) == reduce(mul, (a**b for (a, b) in pff(100)))
    
  7. Will Ness said

    Assuming `primes` is already implemented, returning a stream of prime numbers in order, in Haskell,

    ff n = [(p, sum . takeWhile (> 0) . tail $ iterate (`div` p) n) | p <- takeWhile (<= n) primes]
    
    -- Prelude> ff 33
    -- [(2,31),(3,15),(5,7),(7,4),(11,3),(13,2),(17,1),(19,1),(23,1),(29,1),(31,1)]
    
  8. Paul said

    Another Python version. The optimization suggested by Will Ness did not work (in Python).

    def facfac(n):
        for p in primes(n):
            nfacs = 0
            f = n // p
            while f:
                nfacs += f
                f //= p
            yield p, nfacs
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 634 other followers

%d bloggers like this: