# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1741205374 18000
# Node ID dee136fa9a816c31ba2a46aacc49307c0b010e34
# Parent  d5710762069c01068d8194553bc455814ea67d79
Refactor, add pf, add units

diff -r d5710762069c -r dee136fa9a81 cacl.asd
--- a/cacl.asd	Wed Mar 27 10:22:11 2024 -0400
+++ b/cacl.asd	Wed Mar 05 15:09:34 2025 -0500
@@ -17,4 +17,5 @@
                (:module "src" :serial t
                 :components ((:file "package")
                              (:file "base")
-                             (:file "math")))))
+                             (:file "math")
+                             (:file "units")))))
diff -r d5710762069c -r dee136fa9a81 src/base.lisp
--- a/src/base.lisp	Wed Mar 27 10:22:11 2024 -0400
+++ b/src/base.lisp	Wed Mar 05 15:09:34 2025 -0500
@@ -24,7 +24,7 @@
   (pop *stack*))
 
 (defun pop-all! ()
-  (prog1 *stack* (setf *stack* nil)))
+  (prog1 (copy-list *stack*) (setf *stack* nil)))
 
 
 (defmacro with-args (symbols &body body)
@@ -290,9 +290,9 @@
            (setf *stack* ,old-stack))))))
 
 
-(defun read-input ()
+(defun read-input (stream)
   (let ((*read-default-float-format* 'double-float)
-        (line (read-line *standard-input* nil :eof nil)))
+        (line (read-line stream nil :eof nil)))
     (if (eq :eof line)
       (setf *running* nil)
       (read-all-from-string line))))
@@ -306,8 +306,8 @@
         (cons (apply 'special-form input)))
       (save-stack))))
 
-(defun handle-all-input ()
-  (mapc #'handle-input (read-input)))
+(defun handle-all-input (stream)
+  (mapc #'handle-input (read-input stream)))
 
 
 (defun print-stack (&key (stack *stack*) (decorated t))
@@ -325,18 +325,22 @@
   (force-output))
 
 
-(defun run (&key (batch nil))
+(defun run (&key input (batch nil) (slurp nil))
   (setf *running* t
         *stack* nil
         *previous* (list nil))
   (let ((*package* (find-package :cacl)))
+    (when slurp
+      (iterate (while *running*)
+               (handle-all-input slurp))
+      (setf *running* t))
     (iterate (while *running*)
              (progn
                (unless batch
                  (terpri)
                  (print-stack)
                  (print-prompt))
-               (handle-all-input)))
+               (handle-all-input input)))
     (when batch
       (print-stack :decorated nil)))
   (values))
@@ -376,7 +380,7 @@
     :result-key 'batch
     :help "run in interactive mode (the default)"
     :long "interactive"
-    :short #\i
+    :short #\B
     :initial-value nil
     :reduce (constantly nil)))
 
@@ -388,19 +392,28 @@
     :short #\b
     :reduce (constantly t)))
 
+(defparameter *o-no-inform*
+  (adopt:make-option 'no-inform
+    :result-key 'inform
+    :help "suppress printing of informational message at startup"
+    :long "no-inform"
+    :short #\i
+    :reduce (constantly nil)))
+
 (defparameter *o-inform*
   (adopt:make-option 'inform
     :help "print informational message at startup (the default)"
     :long "inform"
     :initial-value t
+    :short #\I
     :reduce (constantly t)))
 
-(defparameter *o-no-inform*
-  (adopt:make-option 'no-inform
-    :result-key 'inform
-    :help "suppress printing of informational message at startup"
-    :long "no-inform"
-    :reduce (constantly nil)))
+(adopt:defparameters (*o-slurp* *o-no-slurp*)
+  (adopt:make-boolean-options 'slurp
+    :help "slurp input from stdin before continuing to a REPL"
+    :help-no "don't slurp input from stdin (the default)"
+    :long "slurp"
+    :short #\s))
 
 
 (defparameter *ui*
@@ -413,10 +426,12 @@
      (list *o-help*
            *o-rcfile*
            *o-no-rcfile*
+           *o-batch*
            *o-interactive*
-           *o-batch*
+           *o-no-inform*
            *o-inform*
-           *o-no-inform*)))
+           *o-slurp*
+           *o-no-slurp*)))
 
 
 (defun input-source (arguments)
@@ -435,7 +450,15 @@
         (print-version))
       (when-let ((rc (gethash 'rcfile options)))
         (load rc :if-does-not-exist nil))
-      (let ((*standard-input* (input-source arguments)))
-        (run :batch (gethash 'batch options))))))
+      (flet ((run% (slurp) ; todo this is gross 
+               (run :input (input-source arguments)
+                    :batch (gethash 'batch options)
+                    :slurp slurp)))
+        (if (gethash 'slurp options)
+          (with-open-file (human-input "/dev/tty" :direction :input)
+            (let ((slurp *standard-input*)
+                  (*standard-input* human-input))
+              (run% slurp)))
+          (run% nil))))))
 
 
diff -r d5710762069c -r dee136fa9a81 src/math.lisp
--- a/src/math.lisp	Wed Mar 27 10:22:11 2024 -0400
+++ b/src/math.lisp	Wed Mar 05 15:09:34 2025 -0500
@@ -1,5 +1,6 @@
 (in-package :cacl)
 
+;;;; Misc ---------------------------------------------------------------------
 (defun cube (number) (* number number number))
 
 (defun factorial (number)
@@ -13,6 +14,259 @@
            (multiplying (/ (- (1+ n) i) i))))
 
 
+;;;; Primes -------------------------------------------------------------------
+(declaim (ftype (function ((integer 0 #.array-dimension-limit))
+                          (simple-array bit (*)))
+                sieve%))
+
+(eval-dammit
+  (defun sieve% (limit)
+    "Return a bit vector of primality for all odd numbers below `limit`."
+    (iterate
+      (with length = (truncate limit 2))
+      (with numbers = (make-array length :initial-element 1 :element-type 'bit))
+      (for bit :in-vector numbers :with-index n :from 1)
+      (when (= 1 bit)
+        (iterate
+          (with step = (1+ (* 2 n)))
+          (for composite :from (+ n step) :by step :below length)
+          (setf (aref numbers composite) 0)))
+      (finally
+        (setf (aref numbers 0) 0)
+        (return numbers)))))
+
+(defun sieve (limit)
+  "Return a vector of all primes below `limit`."
+  (declare (optimize speed))
+  (check-type limit (integer 0 #.array-dimension-limit))
+  (if (< limit 2)
+    (vector)
+    (iterate
+      (declare (iterate:declare-variables))
+      (if-first-time
+        (collect 2 :into result :result-type vector))
+      (for bit :in-vector (sieve% limit) :with-index n :from 1)
+      (when (= 1 bit)
+        (collect (1+ (* 2 n)) :into result :result-type vector))
+      (finally (return result)))))
+
+
+(defun prime-factors (n)
+  "Return the prime factors of `n`.
+
+  The result will be a set of primes.  For example:
+
+    (prime-factors 60)
+    ; =>
+    (2 3 5)
+
+  "
+  (let ((result (list)))
+    (iterate (while (evenp n)) ; handle 2, the only even prime factor
+             (if-first-time
+               (push 2 result))
+             (setf n (/ n 2)))
+    (iterate
+      (for i :from 3 :to (sqrt n) :by 2) ; handle odd (prime) divisors
+      (iterate (while (dividesp n i))
+               (if-first-time
+                 (push i result))
+               (setf n (/ n i))))
+    (when (> n 2) ; final check in case we ended up with a prime
+      (push n result))
+    (nreverse result)))
+
+(defun prime-factorization (n)
+  "Return the prime factorization of `n` as a flat list.
+
+  The result will be a list of primes, with duplicates.  For example:
+
+    (prime-factorization 60)
+    ; =>
+    (2 2 3 5)
+
+  "
+  ;; from http://www.geeksforgeeks.org/print-all-prime-factors-of-a-given-number/
+  (let ((result (list)))
+    (iterate (while (evenp n)) ; handle 2, the only even prime factor
+             (push 2 result)
+             (setf n (/ n 2)))
+    (iterate
+      (for i :from 3 :to (sqrt n) :by 2) ; handle odd (prime) divisors
+      (iterate (while (dividesp n i))
+               (push i result)
+               (setf n (/ n i))))
+    (when (> n 2) ; final check in case we ended up with a prime
+      (push n result))
+    (nreverse result)))
+
+(defun prime-factorization-pairs (n &key include-zeros)
+  "Return the prime factorization of `n` as a list of prime/exponent conses.
+
+  The result will be a list of `(prime . exponent)` conses.  If `include-zeros`
+  is given, even primes whose exponent turns out to be zero will be included.
+
+  For example:
+
+    (prime-factorization-pairs 60)
+    ; =>
+    ((2 . 2) (3 . 1) (5 . 1))
+
+  "
+  (let ((result (list)))
+    (iterate (while (evenp n)) ; handle 2, the only even prime factor
+             (else (when include-zeros
+                     (push (cons 2 0) result)))
+             (if-first-time
+               (push (cons 2 1) result)
+               (incf (cdar result)))
+             (setf n (/ n 2)))
+    (iterate
+      (for i :from 3 :to (sqrt n) :by 2) ; handle odd (prime) divisors
+      (iterate (while (dividesp n i))
+               (else (when (and include-zeros (primep i))
+                       (push (cons i 0) result)))
+               (if-first-time
+                 (push (cons i 1) result)
+                 (incf (cdar result)))
+               (setf n (/ n i))))
+    (when (> n 2) ; final check in case we ended up with a prime
+      (push (cons n 1) result))
+    (nreverse result)))
+
+
+(defun expmod (base exp m)
+  "Return base^exp % m quickly."
+  ;; From SICP and
+  ;; https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Common_Lisp
+  ;;
+  ;; We want to avoid bignums as much as possible.  This computes (base^exp % m)
+  ;; without having to deal with huge numbers by taking advantage of the fact
+  ;; that:
+  ;;
+  ;;     (x * y) % m
+  ;;
+  ;; is equivalent to:
+  ;;
+  ;;     ((x % m) * (y % m)) % m
+  ;;
+  ;; So for the cases where `exp` is even, we can split base^exp into an x and
+  ;; y both equal to base^(exp/2) and use the above trick to handle them
+  ;; separately.  Even better, we can just compute it once and square it.
+  ;;
+  ;; We also make it tail recursive by keeping a running accumulator:
+  ;;
+  ;;    base^exp * acc
+  (labels
+      ((recur (base exp acc)
+         (cond
+           ((zerop exp) acc)
+           ((evenp exp)
+            (recur (rem (square base) m)
+                   (/ exp 2)
+                   acc))
+           (t
+            (recur base
+                   (1- exp)
+                   (rem (* base acc) m))))))
+    (recur base exp 1)))
+
+
+(defun factor-out (n factor)
+  "Factor the all the `factor`s out of `n`.
+
+  Turns `n` into:
+
+    factor^e * d
+
+  where `d` is no longer divisible by `n`, and returns `e` and `d`.
+
+  "
+  (iterate (for d :initially n :then (/ d factor))
+           (for e :initially 0 :then (1+ e))
+           (while (dividesp d factor))
+           (finally (return (values e d)))))
+
+(defun miller-rabin-prime-p (n &optional (k 11))
+  "Return whether `n` might be prime.
+
+  If `t` is returned, `n` is probably prime.
+  If `nil` is returned, `n` is definitely composite.
+
+  "
+  ;; https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Common_Lisp
+  (cond
+    ((< n 2) nil)
+    ((< n 4) t)
+    ((evenp n) nil)
+    (t (multiple-value-bind (r d)
+           (factor-out (1- n) 2)
+         (flet ((strong-liar-p (a)
+                  (let ((x (expmod a d n)))
+                    (or (= x 1)
+                        (iterate (repeat r)
+                                 (for y :initially x :then (expmod y 2 n))
+                                 (when (= y (1- n))
+                                   (return t)))))))
+           (iterate (repeat k)
+                    (for a = (random-range-exclusive 1 (1- n)))
+                    (always (strong-liar-p a))))))))
+
+(defun brute-force-prime-p (n)
+  "Return (slowly) whether `n` is prime."
+  (cond
+    ((or (= n 0) (= n 1)) nil)
+    ((= n 2) t)
+    ((evenp n) nil)
+    (t (iterate (for divisor :from 3 :to (sqrt n))
+                (when (dividesp n divisor)
+                  (return nil))
+                (finally (return t))))))
+
+
+;;; We precompute a bit vector of the primality of the first few odd prime
+;;; numbers to make checking primes faster.
+
+(defconstant +precomputed-primality-limit+ 10000000)
+
+(defparameter *precomputed-primality-bit-vector*
+  (sieve% +precomputed-primality-limit+))
+
+(deftype precomputed-primality-bit-vector ()
+  `(simple-array bit (,(truncate +precomputed-primality-limit+ 2))))
+
+(deftype integer-with-precomputed-primality ()
+  `(integer 0 (,+precomputed-primality-limit+)))
+
+(defun-inline precomputed-prime-p% (n)
+  (declare (optimize speed (debug 0) (safety 1))
+           (type integer-with-precomputed-primality n)
+           (type precomputed-primality-bit-vector *precomputed-primality-bit-vector*))
+  (not (zerop (aref *precomputed-primality-bit-vector* (truncate n 2)))))
+
+
+(defun primep (n)
+  "Return (less slowly) whether `n` is prime."
+  (cond
+    ;; short-circuit a few edge/common cases
+    ((< n 2) nil)
+    ((= n 2) t)
+    ((evenp n) nil)
+    ((< n +precomputed-primality-limit+) (precomputed-prime-p% n))
+    (t (miller-rabin-prime-p n))))
+
+(defun compositep (n)
+  "Return whether `n` is composite."
+  (and (not (= 1 n))
+       (not (primep n))))
+
+(defun coprimep (a b)
+  (= 1 (gcd a b)))
+
+
+
+
+;;;; Commands -----------------------------------------------------------------
 
 (define-constant-command e (exp 1.0d0))
 (define-constant-command pi pi)
@@ -122,6 +376,26 @@
   "Pop the entire stack, multiply everything together, and push the result."
   (push! (product (pop-all!))))
 
+(define-command (mean average) ()
+  "Pop the entire stack, compute the mean, and push the result."
+  (let* ((xs (pop-all!))
+         (n (length xs)))
+    (when (zerop n)
+      (error "Cannot compute mean of empty sequence."))
+    (push! (/ (summation xs) n))))
+
+(define-command median ()
+  "Pop the entire stack, compute the median, and push the result."
+  (let* ((xs (sort (pop-all!) #'<))
+         (n (length xs)))
+    (when (zerop n)
+      (error "Cannot compute median of empty sequence."))
+    (push! (if (oddp n)
+             (nth (floor n 2) xs)
+             (/ (+ (nth (1- (/ n 2)) xs)
+                   (nth (/ n 2) xs))
+                2)))))
+
 (define-command log (base x)
   (push! (log x base)))
 
@@ -133,50 +407,29 @@
   (push! (- x y))
   (push! (+ x y)))
 
-(define-command kb (bytes)
-  "Convert bytes to kilobytes."
-  (push! (coerce (/ bytes (expt 1024 1)) 'double-float)))
+(define-command primep (x)
+  (push! (primep x)))
 
-(define-command mb (bytes)
-  "Convert bytes to megabytes."
-  (push! (coerce (/ bytes (expt 1024 2)) 'double-float)))
-
-(define-command gb (bytes)
-  "Convert bytes to gigabytes."
-  (push! (coerce (/ bytes (expt 1024 3)) 'double-float)))
+(define-command coprimep (x y)
+  (push! (coprimep x y)))
 
-(define-command tb (bytes)
-  "Convert bytes to terabytes."
-  (push! (coerce (/ bytes (expt 1024 4)) 'double-float)))
+(define-command compositep (x)
+  (push! (compositep x)))
 
-(define-command pb (bytes)
-  "Convert bytes to petabytes."
-  (push! (coerce (/ bytes (expt 1024 5)) 'double-float)))
-
-(define-command eb (bytes)
-  "Convert bytes to exabytes."
-  (push! (coerce (/ bytes (expt 1024 6)) 'double-float)))
+(define-command (prime-factors factors) (x)
+  (dolist (f (prime-factors x))
+    (push! f)))
 
-(define-command kbb (bytes)
-  "Convert kilobytes to bytes."
-  (push! (* bytes (expt 1024 1))))
-
-(define-command mbb (bytes)
-  "Convert megabytes to bytes."
-  (push! (* bytes (expt 1024 2))))
-
-(define-command gbb (bytes)
-  "Convert gigabytes to bytes."
-  (push! (* bytes (expt 1024 3))))
+(define-command (prime-factorization factorization fz) (x)
+  (dolist (f (prime-factorization x))
+    (push! f)))
 
-(define-command tbb (bytes)
-  "Convert terabytes to bytes."
-  (push! (* bytes (expt 1024 4))))
-
-(define-command pbb (bytes)
-  "Convert petabytes to bytes."
-  (push! (* bytes (expt 1024 5))))
-
-(define-command ebb (bytes)
-  "Convert exabytes to bytes."
-  (push! (* bytes (expt 1024 6))))
+(define-command pf (x)
+  "Print x as a nicely-formatted number."
+  (multiple-value-bind (ipart fpart) (ftruncate x)
+    (let ((ipart (round ipart))
+          (fpart (if (zerop fpart)
+                   ""
+                   (subseq (format nil "~F" (abs fpart)) 1))))
+      (format t "~:D~A" ipart fpart)))
+  (values))
diff -r d5710762069c -r dee136fa9a81 src/units.lisp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/units.lisp	Wed Mar 05 15:09:34 2025 -0500
@@ -0,0 +1,57 @@
+(in-package :cacl)
+
+(define-command bkb (bytes)
+  "Convert bytes to kilobytes."
+  (push! (coerce (/ bytes (expt 1024 1)) 'double-float)))
+
+(define-command bmb (bytes)
+  "Convert bytes to megabytes."
+  (push! (coerce (/ bytes (expt 1024 2)) 'double-float)))
+
+(define-command bgb (bytes)
+  "Convert bytes to gigabytes."
+  (push! (coerce (/ bytes (expt 1024 3)) 'double-float)))
+
+(define-command btb (bytes)
+  "Convert bytes to terabytes."
+  (push! (coerce (/ bytes (expt 1024 4)) 'double-float)))
+
+(define-command bpb (bytes)
+  "Convert bytes to petabytes."
+  (push! (coerce (/ bytes (expt 1024 5)) 'double-float)))
+
+(define-command beb (bytes)
+  "Convert bytes to exabytes."
+  (push! (coerce (/ bytes (expt 1024 6)) 'double-float)))
+
+(define-command kbb (bytes)
+  "Convert kilobytes to bytes."
+  (push! (* bytes (expt 1024 1))))
+
+(define-command mbb (bytes)
+  "Convert megabytes to bytes."
+  (push! (* bytes (expt 1024 2))))
+
+(define-command gbb (bytes)
+  "Convert gigabytes to bytes."
+  (push! (* bytes (expt 1024 3))))
+
+(define-command tbb (bytes)
+  "Convert terabytes to bytes."
+  (push! (* bytes (expt 1024 4))))
+
+(define-command pbb (bytes)
+  "Convert petabytes to bytes."
+  (push! (* bytes (expt 1024 5))))
+
+(define-command ebb (bytes)
+  "Convert exabytes to bytes."
+  (push! (* bytes (expt 1024 6))))
+
+(define-command c2f (celsius)
+  "Convert Celsius to Fahrenheit."
+  (push! (+ (* 9/5 celsius) 32)))
+
+(define-command f2c (fahrenheit)
+  "Convert Fahrenheit to Celsius."
+  (push! (* (- fahrenheit 32) 5/9)))