primes.el

;;; primes.el --- prime number support for Emacs Lisp library code

;; Copyright (C) 1999
;;        Free Software Foundation, Inc.

;; This file is part of GNU Emacs.

;; GNU Emacs is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2, or (at your option)
;; any later version.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
;; along with GNU Emacs; see the file COPYING.  If not, write to the
;; Free Software Foundation, Inc., 59 Temple Place - Suite 330,
;; Boston, MA 02111-1307, USA.

;; Author: Nelson H. F. Beebe <beebe@math.utah.edu>
;; Maintainer: Nelson H. F. Beebe <beebe@math.utah.edu>
;; Created: 27 February 1999
;; Version: 1.00
;; Keywords: prime numbers, primality testing

;; This file is part of GNU Emacs.

;;; Commentary:

;; A prime number is any integer greater than one which has no exact
;; integer divisors other than one and itself.
;;
;; Prime numbers have increasingly important practical applications
;; in cryptography, and are also useful in hashing, besides being of
;; fundamental importance in number theory.
;;
;; This is a small collection of functions for:
;;
;;  * testing integers for primality,
;;  * generating nearby primes,
;;  * finding the n-th prime,
;;  * generating lists of primes in a given range,
;;  * factoring a number into a product of primes,
;;  * finding the greatest common divisor of two numbers, and
;;  * finding the least common multiple of two numbers.
;;
;; The functions provided are:
;;
;;  (gcd n m)           [cost: O((12(ln 2)/pi^2)ln n)]
;;  (lcm n m)           [cost: O((12(ln 2)/pi^2)ln n)]
;;  (next-prime n)          [cost: O(sqrt(N))]
;;  (nth-prime n)           [cost: O(N*sqrt(N))]
;;  (prev-prime n)          [cost: O(sqrt(N))]
;;  (prime-factors n)       [cost: O(N)]
;;  (prime-p n)             [cost: O(sqrt(N))]
;;  (primes-between from to)    [cost: O((to - from + 1)*sqrt(N))]
;;  (this-or-next-prime n)      [cost: O(sqrt(N))]
;;  (this-or-prev-prime n)      [cost: O(sqrt(N))]
;;
;; The modest collection of functions here is likely to grow, and
;; perhaps may even be improved algorithmically.  The core of most of
;; these functions is the primality test, (prime-p N), whose running
;; time is O(sqrt(N)), which becomes excessive for large N.  Note that
;; sqrt(N) == 2^{(lg N)/2}, where lg N, the base-2 logarithm of N, is
;; the number of bits in N.  Thus O(sqrt(N)) means O(2^(bits in N)),
;; or O(10^(digits in N)).  That is, the running time increases
;; EXPONENTIALLY in the number of digits of N.
;;
;; Because knowledge of the cost of these functions may be critical to
;; the caller, each function's documentation string ends with a
;; bracketed cost estimate as a final paragraph.
;;
;; Faster algorithms capable of dealing with larger integers are
;; known.  For example, Maple V Release 5 (1997) implements a
;; probabilistic function, isprime(n), that is
;;
;;  ``very probably'' prime - see Knuth ``The Art of Computer
;;  Programming'', Vol 2, 2nd edition, Section 4.5.4, Algorithm P
;;  for a reference and H. Reisel, ``Prime numbers and computer
;;  methods for factorization''. No counter example is known and
;;  it has been conjectured that such a counter example must be
;;  hundreds of digits long.
;;
;; Robert Sedgewick also promises a fast prime test in Part 5 of his
;; book ``Algorithms in C'', not yet published at the time of writing
;; this in March 1999.
;;
;; Three algorithms for probabilistic primality tests for large
;; numbers are discussed in Bruce Schneier, ``Applied Cryptography'',
;; (Wiley, 1994, ISBN 0-471-59756-2), pp. 213--216.
;;
;; Other key references, described in more detail in the Emacs Lisp
;; Manual chapter for this library, include
;;
;;  Leonard M. Adleman, Algorithmic Number Theory --- The
;;  Complexity Contribution, Proc. 35th IEEE Symposium on the
;;  Foundations of Computer Science (FOCS'94), Shafi Goldwasser
;;  (Ed.), IEEE Computer Society Press (Silver Spring, MD),
;;  pp. 88--113, 1994, ISBN 0-8186-6582-3, ISSN 0272-5428.
;;
;;  Eric Bach and Jeffrey Shallit, Algorithmic Number Theory.
;;  Volume I: Efficient Algorithms, MIT Press (Cambridge, MA),
;;  1996, ISBN 0-262-02405-5.
;;
;;  Ronald L. Graham, Donald E. Knuth and Oren Patashnik, Concrete
;;  Mathematics, Addison-Wesley, Reading, MA, USA, 1989, ISBN
;;  0-201-14236-8.
;;
;;  Donald E. Knuth, Fundamental algorithms, The Art of Computer
;;  Programming, Volume 1, Third edition, Addison-Wesley (Reading,
;;  MA), 1997, ISBN 0-201-89683-4.
;;
;;  Donald E. Knuth, Seminumerical algorithms, The Art of Computer
;;  Programming, Volume 2, Third edition, Addison-Wesley (Reading,
;;  MA), 1997, ISBN 0-201-89684-2.
;;
;;  Steven S. Skiena, The Algorithm Design Manual, Springer-Verlag
;;  (New York, NY), 1998, ISBN 0-387-94860-0.
;;
;;; Code:

(provide 'primes)


(defconst primes-version "1.01"
  "Version number of primes library.")


(defconst primes-date "[11-Jan-2002]"
  "Revision date of primes library.")


(defun gcd (m n)
  "Return the Greatest Common Divisor of integers M and N, or nil if
they are invalid.

\[cost: O((12(ln 2)/pi^2)ln max(M,N)) == O(0.8427659... max(M,N))]"
  ;; For details of this 2300-year algorithm due to Euclid, see, e.g.
  ;; Ronald L. Graham, Donald E. Knuth and Oren Patashnik, ``Concrete
  ;; Mathematics'' (Addison-Wesley, Reading, MA, USA, 1989, ISBN
  ;; 0-201-14236-8), pp. 103--104.
  ;;
  ;; The complexity analysis is surprisingly difficult; see Donald
  ;; E. Knuth, ``Seminumerical algorithms, The Art of Computer
  ;; Programming, Volume 2'', third edition (Addison-Wesley, Reading,
  ;; MA, USA, 1997, ISBN 0-201-89684-2), pp. 356--373.

  (if (and (integerp m) (integerp n))   ; check for integer args
      (progn
        (setq m (abs m))        ; argument sign does not matter for gcd, so
        (setq n (abs n))        ; force positive for the algorithm below
        (cond
         ((and (= m 0) (= n 0)) 0)  ; gcd(0,0) ==> 0 by definition, for convenience
         ((and (= m 0) (> n 0)) n)  ; gcd(0,n) ==> n
         ((and (> m 0) (= n 0)) m)  ; gcd(m,0) ==> m
         ((and (> m n) (> n 0)) (gcd n m)) ; reinvoke with reversed arguments
         (t (gcd (% n m) m))))      ; else reduce recursively
    nil))               ; non-integer args: invalid


(defun lcm (m n)
  "Return the Least Common Multiple of integers M and N, or nil if
they are invalid, or the result is not representable (e.g., the
product M*N overflows).

\[cost: O((12(ln 2)/pi^2)ln max(M,N)) == O(0.8427659... max(M,N))]"
  (cond
   ((and (integerp m) (integerp n))     ; check for integer args
    (let ((mn) (the-gcd) (the-lcm))
      (if (or (= m 0) (= n 0))      ; fast special case: lcm(0,anything) == 0
          0
        ;; else compute lcm from (m * n) / gcd(m,n)
        ;;
        ;; Problem: GNU Emacs Lisp integer multiply does not detect or
        ;; trap overflow, which is a real possibility here, and it lacks
        ;; a double-length integer type to represent the product m * n.
        ;; Since the lcm may still be representable, we do the
        ;; intermediate computation in (double-precision)
        ;; floating-point, which is still not quite large enough to
        ;; represent all products of Emacs 28-bit integers stored in
        ;; 32-bit words, then convert back to integer results.  The
        ;; floor function will signal an error if the result is not
        ;; representable.  To try to avoid that, we first check that the
        ;; equality gcd * lcm = m * n is satisfied, and only if it is,
        ;; do we invoke floor.
        ;;
        ;; TO DO: find better algorithm without these undesirable
        ;; properties.
        (setq m (abs m))        ; argument sign does not matter for lcm, so
        (setq n (abs n))        ; force positive for the algorithm below
        (setq the-gcd (gcd m n))
        (setq mn (* (float m) (float n)))
        (setq the-lcm (/ mn the-gcd))
        (if (= (* the-gcd the-lcm) mn)  ; then got correct answer
            (floor the-lcm)
          nil))))           ; else out-of-range or invalid
   (t nil)))               ; non-integer args: invalid


(defun prime-factors (n)
  "Return a list of prime factors of N.

If N is prime, there are no factors, except the trivial one of N itself,
so the return value is the list (N).  Thus, if (length (prime-factors N))
is 1, N is prime.

Otherwise, if N is not an integer greater than 1, the return value is
nil, equivalent to an empty list.

\[cost: O(N)]"
  (interactive)
  (let ((result-list nil)
        (n-original n))
    (if (and (integerp n) (> n 1))
        (let ((limit (/ n 2))
              (divisor 2))
          (while (<= divisor limit)
            ;; To correctly handle factors of multiplicity > 1, we must
            ;; be careful to advance the divisor only when it is not a
            ;; factor!  When n is replaced by n/divisor, we can reset
            ;; limit, but only to n/2, not to n.  Consider
            ;; (prime-factors 15): the first factor found is 3, which
            ;; reduces n to 5, which will be the next prime factor
            ;; found, but would be lost if we reset limit to 5/2 == 2.
            ;;
            ;; If this divisor is rejected, as long as it is greater
            ;; than 2, and thus, odd, we can step it by 2, halving the
            ;; number of loop iterations.
            (if (= 0 (% n divisor))
                (setq n (/ n divisor)
                      limit n
                      result-list (append result-list (list divisor)))
              (if (= divisor 2)
                  (setq divisor 3)
                (setq divisor (+ divisor 2)))))
          ;; If we end the while loop with an empty result-list, then
          ;; the input N was prime, so set result-list to a one-element
          ;; list:
          (if (null result-list)
              (setq result-list (list n-original)))))
    result-list))


(defun prime-p (n)
  "Return N if it is a prime number, else nil.

Because Emacs integers are usually more limited in size than the host
word size would suggest, e.g.,

    \[-2^{27}, 2^{27} - 1] == [-134217728, 134217727]

on a 32-bit machine, avoid passing excessively large integers to this
function, otherwise you may experience this failure:

    \(next-prime 134217689)
    Arithmetic domain error: \"sqrt\", -134217728.0

While you may be able to use larger integers on some 64-bit machines,
the required run time for this function is then likely to be excessive.

\[cost: O(sqrt(N))]"
  (interactive)
  (if (integerp n)
      (cond ((< n 2) nil)       ;there are no primes below 2
            ((= n 2) 2)             ;special case for the only even prime
            ((= 0 (% n 2)) nil)         ;there are no other even primes
            (t (catch 'RESULT       ;else we have a positive odd candidate value
                 (let ((limit (floor (sqrt n)))
                       (divisor 2))
                   (while (<= divisor limit)
                     (if (= 0 (% n divisor))
                         (throw 'RESULT nil))
                     (setq divisor (1+ divisor)))
                   n))))
    nil))


(defun next-prime (n)
  "Return the next prime number after N, or else nil.

\[cost: O(sqrt(N))]"
  (if (integerp n)
      (cond
       ((> n 1)
        (let* ((k (if (= 0 (% n 2));start k at next odd number after n
                      (1+ n)
                    (+ n 2))))
          (while (not (prime-p k)) ;and loop upward over odd numbers
            (setq k (+ k 2)))
          k))
       (t 2))
    nil))


(defun nth-prime (n)
  "Return the N-th prime number, or else nil.

The first prime number is 2.

\[cost: O(N*sqrt(N))]"
  (if (integerp n)
      (cond
       ((<= n 0) nil)           ;non-positive args are invalid
       ((= n 1) 2)          ;special case: only even prime
       (t (let ((k 1))          ;n > 1, so test only odd values 3, 5, ...
            (while (> n 1)
              (setq k (+ k 2))
              (if (prime-p k)       ;count down each prime found
                  (setq n (1- n))))
            k)))            ;at loop exit, k is the desired nth-prime
    nil))


(defun prev-prime (n)
  "Return the prime number before (i.e., less than) N, or else nil.

\[cost: O(sqrt(N))]"
  (if (integerp n)
      (cond ((<= n 2) nil)      ;invalid: there are no primes before 2
            ((= n 3) 2)             ;special case: 2 is the only even prime
            (t              ;n > 3
             (let* ((k (if (= 0 (% n 2)) ;start k at prev odd number before n
                           (1- n)
                         (- n 2))))
               (while (not (prime-p k))     ;and loop downward over odd numbers
                 (setq k (- k 2)))
               k)))
    nil))


(defun primes-between (from to)
  "Return a list of prime numbers between FROM and TO, inclusive, else nil.

\[cost: O((to - from + 1)*sqrt(N)/2)]"
  (if (and (integerp from) (integerp to))
      (let ((k)
            (primes '()))
        (if (and (<= from 2) (<= 2 to))     ;handle special case of only even prime
            (setq primes '(2)))
        (setq k (if (= (% from 2) 0) (1+ from) from))
        ;; k is now odd, so we can loop over odd numbers only
        (while (<= k to)
          (if (prime-p k)
              (setq primes (append primes (list k))))
          (setq k (+ k 2)))
        primes)                 ;final successful list
    nil))               ;else invalid arguments


(defsubst this-or-next-prime (n)
  "Return N if it is prime, else return the next prime number after N,
else nil in N is invalid.

\[cost: O(sqrt(N))]"
  (or (prime-p n)
      (next-prime n)))


(defsubst this-or-prev-prime (n)
  "Return N if it is prime, else return the prime number before
(i.e., less than) N.

\[cost: O(sqrt(N))]"
  (or (prime-p n)
      (prev-prime n)))

(defun first-n-primes (n)
  "Return a list of the first N primes ordered least to greatest."
  (let ((curr 0) (res nil))
    (dotimes (_ n)
      (setq curr (next-prime curr))
      (setq res (cons curr res)))
    (nreverse res)))

;;; primes.el ends here