素因数分解パッケージ

習作。素数判定や、多倍長整数の素因数分解をする。

;; Package for math problems about prime numbers


(defpackage :math.prime
  (:use :common-lisp)
  (:export :primep
           :factor
           :afactor))


(in-package :math.prime)



(defparameter *primes*
  (make-array 25
              :fill-pointer t
              :adjustable t
              :initial-contents '( 2  3  5  7 11
                                  13 17 19 23 29
                                  31 37 41 43 47
                                  53 59 61 67 71
                                  73 79 83 89 97))
  "Array of prime numbers.")


(defparameter *squares* (make-hash-table)
  "Hash table of square numbers.")

(setf (gethash 2 *squares*) 4)



(defun square (n)
  "Return n * n."
  (let ((s (gethash n *squares*)))
    (if s s (square-extend n))))

(defun square-extend (n)
  "Extend *squares* by adding n, and return square n."
  (setf (gethash n *squares*) (* n n)))

(defun get-prime (i)
  "Return *primes*[i]."
  (if (< i (fill-pointer *primes*))
      (aref *primes* i)
      (prime-extend i)))

(defun prime-extend (i)
  "Extend *primes* to index i, and return *primes*[i]."
  (do ((j (fill-pointer *primes*) (1+ j)))
      ((> j i) (aref *primes* i))
    (vector-push-extend (find-prime-as-next j) *primes*)))

(defun find-prime-as-next (i)
  "Calculate and return i-th prime number by using *primes*[i-1]."
  (do ((n (+ 2 (aref *primes* (1- i)))
          (+ 2 n)))
      ((next-prime-p n) n)))

(defun next-prime-p (n)
  "Return t if n is a prime number, but n must be greater than 2, and before using this all prime numbers less than n exist in *primes*."
  (do* ((i 1 (1+ i))
        (p (aref *primes* i)
           (aref *primes* i)))
       ((< n (square p)) t)
    (if (zerop (mod n p))
        (return nil))))

(defun primep (n)
  "Return t if n is a prime number, otherwise return nil."
  (cond
    ((< n 2) nil)
    ((= n 2) t)
    ((zerop (/ n 2)) nil)
    ((position n *primes*) t)
    (t (do* ((i 1 (1+ i))
             (p (get-prime i)
                (get-prime i)))
            ((< n (square p)) t)
         (if (zerop (mod n p))
             (return nil))))))

(defun factor (n)
  "Decompose n into prime factors, make a list in ascending order, and return it.
Example: (factor 12) => (2 2 3)"
  (if (< n 2)
      nil
      (do* ((i 0 (1+ i))
            (p 2 (get-prime i))
            (stack nil))
           ((< n (square p))
            (nreverse (cons n stack)))
        (loop
           (if (zerop (mod n p))
               (progn (push p stack)
                      (setq n (/ n p)))
               (return)))
        (if (= n 1)
            (return (nreverse stack))))))

(defun afactor (n)
  "Decompose n into prime factors, make an association list ,each element is (factor . exponent), and return it."
  (if (< n 2)
      nil
      (do* ((i 0 (1+ i))
            (p 2 (get-prime i))
            (stack nil))
           ((< n (square p))
            (nreverse (cons (cons n 1) stack)))
        (when (zerop (mod n p))
          (push (cons p 1) stack)
          (loop
             (setq n (/ n p))
             (if (zerop (mod n p))
                 (incf (cdar stack))
                 (return))))
        (if (= n 1)
            (return (nreverse stack))))))

使用例：

CL-USER> (primep 1234567891)
T
CL-USER> (factor 12345678910111213)
(113 125693 869211457)
CL-USER> (afactor 987654321)
((3 . 2) (17 . 2) (379721 . 1))