prime.lisp

;;; prime factorization
;;; Copyright (C) 2015 SUZUKI Shingo ([email protected])
;;; Lisence: MIT

(defun prime-factorization (num)
  (cond ((= num 0) (list 0))
        ((= num 1) (list 1))
        ((= num -1) (list -1))
        ((< num 0) (cons -1 (prime-factorization (- num))))
        (t (reform-factor-list
             (prime-factorization-prime num 2 nil -1 nil))) ))

(defun product-factor-list (lst)
  (loop with result = 1
        for elem in lst
        when (numberp elem) do (setf result (* result elem))
        when (expt-form-p elem)
        do (setf result (* result (apply #'cl:expt (cdr elem))))
        finally (return result) ))

(defun prime-factorization-prime (org-num prime table index acc)
  (loop for prev-num = org-num then num
        and num = (/ org-num prime) then (/ num prime)
        while (integerp num)
        collect prime into prime-list
        finally (return (prime-factorization-in
                          prev-num table (1+ index)
                          (if prime-list (nconc prime-list acc) acc) ))))

(defun prime-factorization-in (num table index acc)
  (cond ((= num 1) acc)
        ((null table)
         (prime-factorization-in
           num
           (make-prime-table (max-prime-factor-candidate num))
           index acc) )
        ((>= index (length table))
         (cons num acc) )
        (t (prime-factorization-prime
             num (svref table index) table index acc ))))

(defun max-prime-factor-candidate (num)
  (floor (expt num 1/2)) )

(defun make-prime-table (max-num)
  (let* ((table (make-array (truncate (1- max-num) 2) :element-type 'integer))
         (sieve-max-idx
           (truncate (- (max-prime-factor-candidate max-num) 3) 2)) )
    (loop for idx from 0 below (length table)
          do (setf (svref table idx) (+ (* idx 2) 3)) )
    (loop for idx from 0 to sieve-max-idx
          for num across table
          when (integerp num)
            do (loop for idx2 from (1+ idx) below (length table)
                     for num2 = (svref table idx2)
                     when (and num2 (integerp (/ num2 num)))
                       do (setf (svref table idx2) nil) ))
    (minimize-table table) ))

(defun minimize-table (sparse-table)
  (coerce (loop for elem across sparse-table
                when elem collect elem)
          'simple-vector ))

(defun reform-factor-list (rev-lst)
  (let* ((lst (reverse rev-lst))
         (negp (and lst (= (first lst) -1)))
         (result (reform-factor-list-in (if negp (cdr lst) lst) nil)) )
    (if negp (cons -1 result) result) ))

(defun reform-factor-list-in (lst acc)
  (cond ((null lst) (nreverse acc))
        ((null acc) (reform-factor-list-in (cdr lst) (list (car lst))))
        (t (let* ((curr (car lst))
                  (tgt (car acc))
                  (exp-form (build-expt-form curr tgt)) )
             (reform-factor-list-in
               (cdr lst)
               (if exp-form
                 (cons exp-form (cdr acc))
                 (cons curr acc) ))))))

(defun expt-form-p (form)
  (and (listp form)
       (eq (first form) 'cl:expt)
       (= (length form) 3)
       (numberp (second form))
       (numberp (third form)) ))

(defun build-expt-form (num target-form)
  (if (integerp target-form)
    (when
      (= num target-form) `(cl:expt ,num 2))
    (when (and (expt-form-p target-form)
               (= (second target-form) num) )
      `(cl:expt ,num ,(1+ (third target-form))) )))