; factoring factorials

(define (primes n)
  (let ((sieve (make-vector (+ n 1) #t)))
    (let loop ((p 2) (ps (list)))
      (if (< n p) (reverse ps)
        (if (vector-ref sieve p)
            (do ((i (* p p) (+ i p)))
                ((< n i) (loop (+ p 1) (cons p ps)))
              (vector-set! sieve i #f))
            (loop (+ p 1) ps))))))

(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))))))

(display (fact-fact 33)) (newline)

(display
  (let loop ((xs (fact-fact 33)) (fact 1))
    (if (null? xs) fact
      (loop (cdr xs) (* fact (expt (caar xs) (cdar xs)))))))
(newline)

(define (fact n) (if (= n 1) 1 (* n (fact (- n 1)))))

(display (fact 33)) (newline)

((2 . 31) (3 . 15) (5 . 7) (7 . 4) (11 . 3) (13 . 2) (17 . 1) (19 . 1) (23 . 1) (29 . 1) (31 . 1))
8683317618811886495518194401280000000
8683317618811886495518194401280000000