From 236458737ec0ff74c74950a59569d60d070fcaeb Mon Sep 17 00:00:00 2001
From: Zachary Huang <zack466@gmail.com>
Date: Fri, 26 Aug 2022 22:30:16 -0400
Subject: [PATCH] added some extra integer functions

---
 README.md                      |   8 ++
 packages.lisp => builtins.lisp |  20 +++-
 main.lisp                      |   2 +-
 primes.lisp                    | 189 +++++++++++++++++++++++++++++++++
 test.lisp                      |  19 +++-
 util.lisp                      |  52 ++++++++-
 zpcalc.asd                     |   3 +-
 7 files changed, 282 insertions(+), 11 deletions(-)
 rename packages.lisp => builtins.lisp (93%)
 create mode 100644 primes.lisp

diff --git a/README.md b/README.md
index fc52135..87c66a8 100644
--- a/README.md
+++ b/README.md
@@ -271,6 +271,13 @@ It will run code within a package and then return to the previous package.
  - `random` - returns a random integer between 0 (inclusive) and the top stack element (exclusive)
  - `rand` - returns a random float between 0 and 1
  - `isqrt` - returns the integer square root (rounded down) of the top stack element (which must be an integer)
+ - `fib` - returns the nth fibonacci number, where n is the top stack element
+ - `fact` - returns the factorial of the top stack element
+ - `prime` - returns the nth prime number, where n is the top stack element
+ - `primep` - returns 1 if the top stack element is prime, otherwise 0
+ - `totient` - returns phi of the top stack element (also known as the totient function)
+ - `choose` - returns n choose k, where n and k are on top of the stack (also known as the binomial coefficient)
+ - `permute` - returns n permute k, where n and k are on top of the stack
 
 ### Irrational Operations
  - `pow`  - returns a^b, where a and b are the top two stack elements. Tries to preserve exactness
@@ -354,6 +361,7 @@ It will run code within a package and then return to the previous package.
  - `(if)` - a conditional construct that allows for branched execution (see [Conditionals](#conditionals))
  - `(while)` - a construct that allows for looping (see [Looping](#looping))
  - `(in-package)` - enters a package (see [Packages](#packages))
+ - `(with-package)` - executes code within a package (see [Packages](#packages))
 
 ### Top-Level Actions (cannot be evaluated)
  - `quit` - quits the calculator
diff --git a/packages.lisp b/builtins.lisp
similarity index 93%
rename from packages.lisp
rename to builtins.lisp
index 39382b9..7fec98e 100644
--- a/packages.lisp
+++ b/builtins.lisp
@@ -1,6 +1,6 @@
 (in-package :cl-user)
 
-(defpackage zpcalc/packages
+(defpackage zpcalc/builtins
   (:use :cl)
   (:import-from
     #:zpcalc/util
@@ -16,7 +16,14 @@
     #:sqrt-exact
     #:->double
     #:approx-equal
-    #:bool->int)
+    #:bool->int
+    #:factorial
+    #:choose
+    #:permute
+    #:get-nth-prime
+    #:primep
+    #:phi
+    #:fib)
   (:import-from
     #:zpcalc/actions
     #:apply-unary!
@@ -32,7 +39,7 @@
   (:export
     #:*all-packages*
     #:*builtins*))
-(in-package :zpcalc/packages)
+(in-package :zpcalc/builtins)
 
 ;; environment :: symbol -> action
 (defvar *builtins* (make-hash-table))
@@ -113,6 +120,13 @@
 (setf (gethash :SQUARE *builtins*) (apply-unary! #'square))
 (setf (gethash :CUBE *builtins*) (apply-unary! (lambda (x) (* x x x))))
 (setf (gethash :ISQRT *builtins*) (apply-unary! #'isqrt))
+(setf (gethash :FACT *builtins*) (apply-unary! #'factorial))
+(setf (gethash :CHOOSE *builtins*) (apply-binary! #'choose))
+(setf (gethash :PERMUTE *builtins*) (apply-binary! #'permute))
+(setf (gethash :PRIME *builtins*) (apply-unary! #'get-nth-prime))
+(setf (gethash :PRIMEP *builtins*) (apply-unary! (compose #'bool->int #'primep)))
+(setf (gethash :TOTIENT *builtins*) (apply-unary! #'phi))
+(setf (gethash :FIB *builtins*) (apply-unary! #'fib))
 
 ;; irrational operations that try to preserve exactness
 (setf (gethash :POW *builtins*) (apply-binary! #'expt-exact))
diff --git a/main.lisp b/main.lisp
index 4df7432..4455aa2 100644
--- a/main.lisp
+++ b/main.lisp
@@ -7,7 +7,7 @@
     #:zpcalc/util
     #:zpcalc/env
     #:zpcalc/history
-    #:zpcalc/packages)
+    #:zpcalc/builtins)
   (:export
     #:Calc
     #:calc-stack
diff --git a/primes.lisp b/primes.lisp
new file mode 100644
index 0000000..36d6231
--- /dev/null
+++ b/primes.lisp
@@ -0,0 +1,189 @@
+;; ALL OF THE CREDIT FOR ALL THIS CODE GOES TO Yang Xiaofeng
+;; Taken from https://github.com/nakrakiiya/cl-prime-maker
+;; I just wanted to avoid adding quicklisp as a dependency
+
+(in-package :zpcalc/util)
+
+;; Many applications require the use of large prime numbers. The this program can be used to generate large primes and for primality testing.
+;; only big enough numbers are supported
+
+(declaim (optimize (speed 3)))
+
+;; for small prime numbers
+(eval-when (:compile-toplevel :load-toplevel :execute)
+  (defun make-prime-list-for-range (maximum)
+    (declare (type fixnum maximum))
+    (let ((result-array (make-array (list (1+ maximum)) :initial-element t)))
+      ;; init for 0, 1
+      (setf (aref result-array 0) nil)
+      (setf (aref result-array 1) nil)
+      ;; process the rest
+      (loop
+         for base-num from 2 below (1+ (floor maximum 2))
+         do
+           (let ((n (* base-num 2)))
+             (loop
+                (if (<= n maximum)
+                    (progn
+                      (setf (aref result-array n) nil)
+                      (incf n base-num))
+                    (return)))))
+      result-array)))
+
+(defparameter +primes-below-65535+ #.(make-prime-list-for-range 65535))
+
+(defun pow (a b m)  
+  "Computes V = (A^B) mod M. It's much faster than (mod (expt a b) m)."
+  (declare (type integer a b m))
+  (cond
+    ((= b 1)
+     (mod a m))
+    ((= b 2)
+     (mod (* a a) m))
+    (t
+     (let* ((b1 (truncate (/ b 2)))
+            (b2 (- b b1))
+            ;; B2 = B1 or B1+1
+            (p (pow a b1 m)))
+       (if (= b2 b1)
+           (mod (* p p) m)
+           (mod (* p p a) m))))))
+
+;; random:uniform
+(declaim (inline random-uniform))
+(defun random-uniform (n)
+  (declare (type integer n))
+  (1+ (random n)))
+
+;; new_seed
+(declaim (inline new-seed))
+(defun new-seed ()
+  (setq *random-state* (make-random-state t)))
+
+;(declaim (inline make/2))
+(defun make/2 (n d)
+  (declare (type integer n d))
+  (if (= n 0)
+      d
+      (make/2 (1- n) (+ (* 10 d) (1- (random-uniform 10))))))
+
+(defun make (n)
+  "make(n) -> I: Generates a random integer I with N decimal digits. "
+  (new-seed)
+  (make/2 n 0))
+
+;; Fermat's little theorem states that if N is prime then A^N mod N = A. So
+;; to test if N is prime we choose some random A which is less than N and
+;; compute A^N mod N. If this is not equal to A then N is definitely not a
+;; prime. If the test succeeds then A might be a prime (certain composite
+;; numbers pass the Fermat test, these are called pseudo-primes), if we
+;; perform the test over and over again then the probability of mis-classifying
+;; the number reduces by roughly one half each time we perform the test. After
+;; (say) one hundred iterations the probability of mis-classifying a number
+;; is approximately 2^-100. So we can be fairly sure that the classification
+;; is correct.
+
+;(declaim (inline primep/2 primep/3))
+(defun primep/3 (ntest n len)
+  (declare (type integer ntest n len))
+  (if (= ntest 0)
+      t
+      (let* ((k (random-uniform len))
+             ;; A is a random number less than N
+             (a (make k)))
+        (if (< a n)
+            (when (= a
+                     (pow a n n))
+              (primep/3 (1- ntest) n len))
+            (primep/3 ntest n len)))))
+
+(defun primep/2 (d ntests)
+  (declare (type integer d ntests))
+  (let ((n (1- (length (write-to-string d)))))
+    (primep/3 ntests d n)))
+
+(defun primep (n)
+  "Tests if N is a prime number. Returns T if N is a prime number. Returns NIL otherwise.
+NOTES:
+* If n <= 65535, the detection of whether a number is prime can always get the correct answer.
+* If n > 65535, the detection of whether a number is prime is based on the Fermat's little theorem.
+"
+  (declare (type integer n))
+  (if (<= n 1)
+      nil
+      (if (<= n 65535)
+          (aref +primes-below-65535+ n)
+          (progn (new-seed)
+                 (primep/2 n 100)))))
+
+;(declaim (inline make-prime/2))
+(defun make-prime/2 (k p)
+  (if (= k 0)
+      (error "impossible")
+      (if (primep p)
+          p
+          (make-prime/2 (1- k) (1+ p)))))
+
+(defun make-prime (k)
+  "Generates a random prime P with at least K decimal digits. Returns nil when k <= 0. Returns NIL otherwise. K should be an INTEGER. "
+  (declare (type integer k))
+  (when (> k 0)
+    (new-seed)
+    (let ((n (make k)))
+      (if (> n 3)
+          (let* ((max-tries (- n 3))
+                 (p1 (make-prime/2 max-tries (1+ n))))
+            p1)
+          (make-prime k)))))
+
+
+;; Generating the Nth prime. (S.M.Ruiz 2000)
+;; According to http://zh.wikipedia.org/wiki/%E7%B4%A0%E6%95%B0%E5%85%AC%E5%BC%8F  .
+
+(defvar *ruiz-pis* (make-hash-table))   ; cache pi(k)
+(defvar *ruiz-pis-part1* (make-hash-table)) ; cache a sub-part of pi(k)
+(defvar *ruiz-results* (make-hash-table))   ; cache the results of p(n)
+
+(defun compute-ruiz-pis-part1 (j)
+  (let ((result-from-hash (gethash j *ruiz-pis-part1*)))
+    (if (null result-from-hash)
+        (let* ((s-max (floor (sqrt j)))
+               (sigma1 (loop
+                          for s from 1 to s-max
+                          summing
+                            (- (floor (/ (1- j) s))
+                               (floor (/ j s)))))
+               (result (floor (* (/ 2 j)
+                                 (1+ sigma1)))))
+          (setf (gethash j *ruiz-pis-part1*) result)
+          result)              
+        result-from-hash)))
+
+(defun compute-ruiz-pi (k)
+  (let ((result-from-hash (gethash k *ruiz-pis*)))
+    (if (null result-from-hash)
+        (let ((result (cond
+                        ((= k 1) 0)
+                        ((= k 2)
+                         (1+ (compute-ruiz-pis-part1 2)))
+                        (t (+ 1
+                              (compute-ruiz-pis-part1 k)
+                              (compute-ruiz-pi (1- k)))))))
+          (setf (gethash k *ruiz-pis*) result)
+          result)
+        result-from-hash)))
+
+(defun get-nth-prime (n)
+  (declare (type integer n))
+  "Generate the Nth prime number when N >= 1. Otherwise this function always returns 2."
+  (if (>= n 1)
+      (let ((result-from-hash (gethash n *ruiz-results*)))
+        (if (null result-from-hash)
+            (let ((result (1+ (loop
+                                 for k from 1 to (* 2 (1+ (floor (* n (log n)))))
+                                 summing (- 1 (floor (/ (compute-ruiz-pi k)
+                                                        n)))))))
+              (setf (gethash n *ruiz-results*) result)
+              result)
+            result-from-hash))
+      2))
diff --git a/test.lisp b/test.lisp
index 1b90fa5..a909368 100644
--- a/test.lisp
+++ b/test.lisp
@@ -109,6 +109,12 @@
 (test-stack square '(-2 square 4 square) '(4 16))
 (test-stack cube '(-2 cube 4 cube) '(-8 64))
 (test-stack isqrt '(4 isqrt 22 isqrt) '(2 4))
+(test-stack fib '(1 fib 10 fib) '(1 55))
+(test-stack fact '(10 fact) '(3628800))
+(test-stack prime '(10 prime) '(29))
+(test-stack totient '(3198 totient) '(960))
+(test-stack choose '(10 2 choose) '(45))
+(test-stack permute '(10 2 permute) '(90))
 
 ;; irrational, tries to preserve exactness
 (test-stack pow '(2 3 pow 2.5 2.5 pow) `(8 ,(expt 2.5d0 2.5d0)))
@@ -136,9 +142,6 @@
 (test-stack clear '(1 2 3 clear) '())
 (test-stack id '(1 2 3 id) '(1 2 3))
 (test-stack sto/rcl '(1 2 sto 3 rcl) '(1 2 3 2))
-(test-stack package '(package) '(:user))
-(test-stack package-enter '(':foo package-enter package) '(:foo))
-(test-stack package-exists '(':bar package-exists) '(0))
 
 ;; top level actions and special constructs
 (deftest undo
@@ -191,8 +194,14 @@
 
 (deftest in-package
   (let ((state (make-instance 'Calc)))
-    (calc-interact state '((in-package baz) package))
-    (assert (equal '(:baz) (reverse (calc-stack state))))))
+    (calc-interact state '((in-package baz)))
+    (assert (equal :baz (calc-package state)))))
+
+(deftest with-package
+  (let ((state (make-instance 'Calc)))
+    (calc-interact state '((with-package foo (def bar 11))))
+    (calc-interact state '(foo.bar))
+    (assert (equal '(11) (reverse (calc-stack state))))))
 
 ;; TODO: test example code
 
diff --git a/util.lisp b/util.lisp
index 87207bb..097147d 100644
--- a/util.lisp
+++ b/util.lisp
@@ -24,7 +24,13 @@
     #:acond
     #:package-designator
     #:approx-equal
-    ))
+    #:factorial
+    #:choose
+    #:permute
+    #:prime
+    #:primep
+    #:phi
+    #:fib))
 (in-package :zpcalc/util)
 
 (defun symbol= (a b)
@@ -94,6 +100,50 @@ before comparison for floats with different exponents."
 (defun truthy (x)
   (not (falsy x)))
 
+(defun factorial (x)
+  (declare (type integer x))
+  (assert (>= x 0))
+  (labels ((rec (x acc)
+                (if (<= x 1)
+                  acc
+                  (rec (1- x) (* x acc)))))
+    (rec x 1)))
+
+;; binomial coefficient
+;; aka number of ways to choose k items out of n
+(defun choose (n k)
+  (assert (>= n k))
+  (/ (factorial n) (* (factorial k) (factorial (- n k)))))
+
+;; number of ways to permute k items out of n
+(defun permute (n k)
+  (assert (>= n k))
+  (/ (factorial n) (*  (factorial (- n k)))))
+
+;; nth fibonacci number using binet's formula
+(defun fib (n)
+  (assert (>= n 0))
+  (round (/ (- (expt (/ (+ 1 (sqrt 5.0d0)) 2.0d0) n)
+               (expt (/ (- 1 (sqrt 5.0d0)) 2.0d0) n)) (sqrt 5.0d0))))
+
+;; totient function
+;; copied algorithm from https://cp-algorithms.com/algebra/phi-function.html
+(defun phi (n)
+  (assert (>= n 1))
+  (let ((result n)
+        (i 2))
+    (loop while (<= (* i i) n)
+          do (progn
+               (when (= 0 (mod n i))
+                 (loop while (= 0 (mod n i))
+                       do (progn
+                            (setf n (truncate (/ n i)))
+                            (decf result (truncate (/ result i))))))
+               (incf i)))
+    (when (> n 1)
+      (decf result (truncate (/ result n))))
+    result))
+
 ;; taken from StackOverflow
 (defun make-keyword (name) (values (intern (string-upcase name) "KEYWORD")))
 
diff --git a/zpcalc.asd b/zpcalc.asd
index 3495970..1ce9c28 100644
--- a/zpcalc.asd
+++ b/zpcalc.asd
@@ -5,12 +5,13 @@
   :depends-on ()
   :serial t
   :components ((:file "util")
+               (:file "primes")
                (:file "env")
                (:file "state")
                (:file "conditions")
                (:file "history")
                (:file "actions")
-               (:file "packages")
+               (:file "builtins")
                (:file "main"))
   :description ""
   :build-operation "program-op"