;;; primes.el
;;; Written by Russell Young (emacs@young-0.com) 2011
;;; 
;;; These are functions to handle prime numbers - finding and factoring. For efficiency,
;;; it assumes there will be multiple calls and so it caches a list of primes for
;;; future use.
;;;
;;; The functions can all be called programatically or interactively. If called 
;;; interactively they print out a message giving the result, otherwise they return it.
;;;
;;; Functions:
;;;  - nth-prime: finds the nth prime
;;;  - which-prime: returns the sequential number of the given number in the list of primes, or nil 
;;;    if it is not a prime
;;;  - factor: returns the prime factorization of the given number
;;;  - primes-less-than: gets the number of primes less than a given number
;;;  - phi-x: gets the ratio of the logs of a number to its greatest prime factor
;;;

(defvar *prime-cache-report-size* 100
  "Report on the prime cache size whenever it grows to a multiple of this.
If nil do not report")

(defvar *prime-cache* '(2.0 3.0)
  "Caches a list of prime numbers to save time
Needs to be salted with 3 to work")

(defvar *cache-increment-size* 50
  "When factoring, grow the prime cache by this size when it needs more primes")

(defun nth-prime (n)
  "Gets Nth prime number
Uses the prime caching system for efficiency"
  (interactive "nNth prime: ")
  (more-primes n t)
  (if (interactive-p) (message "Prime number %s is %s" n (nth n *prime-cache*)))
  (nth n *prime-cache*))

(defun which-prime (p)
  "Returns the sequential number of P in the series of primes, or nil if it is composite
Uses the prime caching system for efficiency."
  (interactive "nNumber: ")
  (setq p (* p 1.0))
  (more-primes p)
  (let* ((position (member p *prime-cache*)))
    (if position (setq position (- (length *prime-cache*) (length position))))
    (if (interactive-p) 
        (if position (message "%s is prime number %s" p position)
          (message "%s is not prime" p)))
    position))

(defun factor (n &optional start-pos)
  "Factors a number N. For large numbers be sure to enter as real rather than int.
Uses the prime caching system for efficiency. Does not grow the cache until
necessary, so factoring very large numbers does not spend time growing the cache
unless there is a very large prime. For numbers larger than emacs maximum int size
be sure to enter the number as a real."
  (interactive "nFactor: ")
  (let* ((factors ()))
    ;; This tests all the primes in the cache currently
    (some (lambda (x) (let ((count 0)) 
                        (while (= 0 (mod n x)) (setq count (1+ count) n (/ n x)))
                        (setq factors (case count
                                        (0 factors)
                                        (1 (cons (floor x) factors))
                                        (t (cons (cons (floor x) count) factors)))))
            (< n 2))
          (nthcdr (or start-pos 0) *prime-cache*))
    (if (> n 1)
        (if (< (sqrt n) (car (last *prime-cache*)))
            (setq factors (cons (floor n) factors))
          (let ((old-cache-len (length *prime-cache*)))
            (more-primes (+ *cache-increment-size* old-cache-len) t)
            (setq factors (append (factor n old-cache-len) factors)))))
    (if (interactive-p) (message "Factors are %s" factors))
    factors))

(defun primes-less-than (n)
  (interactive "nPrimes less than: ")
  (more-primes n)
  (let ((primes (remove nil (mapcar (lambda (x) (if (< x n) x)) *prime-cache*))))
    (if (interactive-p) (message "There are %s primes: %s" (length primes) primes))
    primes))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Neat function to find the "prime degree" of a number:
;; the ratio of the log of the number to the log of its
;; highest prime divisor (thus 0 < PD(x) <= 1, equality
;; for primes)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(defun phi-x (num) 
  (interactive "nNumber: ")
  (let ((result (/ (log10 (car (factor num))) (log10 num))))
    (if (interactive-p) (message (format "phi-x of %s is %s" num result)))
    result))

(defun phis (upto)
  (interactive "nGet phi for numbers up to: ")
  (let* ((a 2)
	 (numbers (mapcar (lambda (x) (setq a (1+ a))) (make-list upto t)))
	 (phis (mapcar (lambda (x) (cons x (phi-x x))) numbers)))
	(mapcar (lambda (x) (insert (format "%s\t%s\n" (car x) (substring "xxxxxxxxxxxxxxxxxxxP" 0 (floor (* 20 (cdr x))))))) phis)))


(defun more-primes (n &optional upto)
  "Internal function to increase the size of the prime cache list
If UPTO is true then the list will be expanded to have at least N primes.
If it is nil then the list will have primes at least up to N."
  (let* ((p (car (last *prime-cache*)))
         (last-cons (last *prime-cache*)))          ; for efficiency manipulate the list directly
    (while (<= (if upto (length *prime-cache*) p) n)
      (setq p (+ 2 p))
      (when (notany (lambda (x) (= (mod (/ p x) 1) 0)) *prime-cache*)
        (setf (cdr last-cons) (list p))
        (setq last-cons (cdr last-cons))
        (if (and *prime-cache-report-size* 
                 (= 0 (mod (length *prime-cache*) *prime-cache-report-size*)))
            (message "Prime cache size %s (up to %s)" (length *prime-cache*) (car last-cons)))
        ))))

;(setq *prime-cache* '(3.0))
;(setq a (factor 12345678901250.0))
;(apply '* a)
;(factor 11)
;(length *prime-cache*)

;;; diagnostic
(defmacro time (form)
  `(let* ((xstart (current-time))
          (result ,form))
     (message "Command took %s seconds" (date-subtract (current-time) xstart))
     (sleep-for 1)
     result))

(defun date-subtract (d1 d2 &optional ms)
  (let ((big (- (car d1) (car d2)))
         (middle (- (second d1) (second d2)))
         (little (- (third d1) (third d2))))
    (if (< little 0) 
        (setq little (+ little 65536)
              middle (1- middle)))
    (if (< middle 0) 
        (setq middle (+ middle 65536)
              big (1- big)))
    (+ (* 65536 big) middle (/ little 1000000.0))))

(provide 'primes)