; This file is ungoverned by the Cypherpunks' Anti-License.
; Do with it as you will.

; Lorraine Lee
; Warren, Michigan, A. S. U. (Alternative State of the Union)
; N8CHZ at YAHOO dot CA

; Kitchen tested on Gnu Common Lisp 2.2.2

; Miscellaneous number theory functions

(defun sopoi (k n) ; sum of k'th powers of integers 1 thru n
 (let ((s 0))
  (dotimes (i n s)
   (setf s (+ (expt (1+ i) k) s)))))

(defun ! (n)
 (let ((f 1))
  (dotimes (i n f)
   (setf f (* (1+ i) f)))))

(defun p (n r) ; permutations of n things taken r at a time
 (let ((p 1))
  (dotimes (i r p)
   (setf p (* (- n i) p)))))

(defun c (n r) ; combinations of n things taken r at a time
 (/ (p n r) (! r)))

; mod is in the most straightforward sense, i.e. using flooring rather than truncation

; truncate is in the forth sense, and returns 2 values (arrgh!)
; also, it uses truncation, so use with negative args is not mathematically sound.

; multiple-value-list is the function to list the values "returned" by truncate

; happily, gcd and lcm are built in

(defun factor (n &optional (test 3) (result nil))
 (let
  ((dm (multiple-value-list (truncate n test))))
  (cond
   ((evenp n) (factor (/ n 2) 3 (cons 2 result)))
   ((> (* test test) n) (cons n result))
   ((zerop (cadr dm)) (factor (car dm) test (cons test result)))
   (t (factor n (+ 2 test) result)))))

(defun phi (n) ; this is the 'totient' function
 (let ((result 0))
  (dotimes (i (- n 1) result)
   (if (= (gcd (1+ i) n) 1) (setf result (1+ result))))))

(defun divp (k n)
 (zerop (mod n k)))

(defun sigma (k n) ; sum of k'th powers of divisors of n
 (let ((s 0))
  (dotimes (i n s)
   (let ((j (1+ i)))
    (if (divp j n) (setf s (+ (expt j k) s)))))))

(defun power-residue (g b p)
 (mod (expt g b) p))

(defun partition (n &optional (k 1)) ; calculate number of ways to express n as sums of positive integers
 (cond
  ((> k n) 0)
  ((= k n) 1)
  (t (+ (partition (- n k) k) (partition n (1+ k))))))

(defun partitions (size) ; generate partition function table (much faster than recursive method above)
 (let* ((h (1+ (truncate size 3))) (p (make-array (list (1+ h) (1+ size)))))
  (do
   ((n 1 (1+ n)))
   ((> n size) p)
   (do
    ((k h (1- k)))
    ((< k 1))
    (setf (aref p k n)
     (cond
      ((> k n) 0)
      ((> (* 2 k) n) 1)
      ((> (* 3 k) n) (- (truncate n 2) k -2)) ; I have yet to prove this assertion.  Can you?
      (t (+ (aref p (1+ k) n) (aref p k (- n k))))))))))

(defun digit (place n &optional (base 10)) ; digit correspinding to base^place, i.e. 0,1,2...counting from right
 (do
  ((i place (1- i))
   (dm (multiple-value-list (truncate n base)) (multiple-value-list (truncate (car dm) base))))
  ((< i 1) (cadr dm))))

(defun squarep (n)
 (let ((s (truncate (sqrt n))))
  (if (= (* s s) n) s)))

(defun pq (n) ; Pell's coordinates
 (if
  (not (squarep n))
  (let ((p 2))
   (do
    ((q 1 (1+ q)))
    ((setf p (squarep p)) (list p (1- q)))
    (setf p (1+ (* n q q)))))))

(defun frac (x) (- x (floor x)))

(defun cf (x &optional (limit 40)) ; continued fraction
 (let ((r (list (floor x))))
  (if (/= (car r) x)
   (do
    ((y (/ (frac x)) (/ (frac y))))
    ((or
      (zerop (frac y))
      (and (cdr r) (= (* 2 (car (last r))) (car r))) ; comment in for pell coordinates
      (> (length r) limit))
     (reverse r))
    (setf r (cons (floor y) r))))))

(defun uncf (f) ; undo action of cf, above
 (let ((x (car (last f))))
  (do
   ((rest (cdr (reverse f)) (cdr rest)))
   ((null rest) x)
   (setf x (+ (car rest) (/ x))))))

(defun cftest (n) ; put cf through its paces
 (let ((r nil))
  (dotimes (i (1+ n) (reverse r))
   (setf r (cons (cf (sqrt i)) r)))))

(defun pell (n) ; return convergent of sqrt(n) corresponding to pell's coordinates of n
 (let ((r nil))
  (dotimes (i (1+ n) (reverse r))
   (setf r (cons (uncf (reverse (cdr (reverse (cf (sqrt i)))))) r)))))

(setf pie (cf 3141592653589793238462643383279419716939937510/1000000000000000000000000000000000000000000000))