--- 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")))))
--- 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))))))
--- 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))
--- /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)))