; uLisp Infinite Precision Package with Division - 11th March 2026
; see http://www.ulisp.com/show?1FHA
;

; Convert to a big number

(defun big (x) (if (stringp x) (bigstr x) (bign x)))

; Convert an integer n to a big number

(defun bign (n)
  (cond
   ((zerop n) nil)
   ((< n 100) (list n))
   (t (cons
       (mod n 100)
       (bign (truncate n 100))))))

; Convert a string s to a big number

(defun bigstr (s)
  (let ((ls (length s)) result)
    (dotimes (n (ceiling ls 2) (reverse result))
      (push 
       (read-from-string 
        (subseq s (max (- ls (* 2 n) 2) 0) (- ls (* 2 n))))
       result))))

; Print a big number

(defun pri (n &optional string)
  (let ((rn (reverse n)))
    (or 
     (format (not string) "~d~{~2,'0d~}~%" (car rn) (cdr rn)) 
     nothing)))

; Normalise - remove trailing zeros

(defun norm (a)
  (cond
   ((null a) nil)
   (t (let ((rest (norm (cdr a))))
        (cond
         ((and (null rest) (zerop (car a))) nil)
         (t (cons (car a) rest)))))))

; Add two big numbers a and b.

(defun add (a b)
  (cond
   ((null a) b)
   ((null b) a)
   (t (let ((n (+ (car a) (car b)))
            (c (add (cdr a) (cdr b))))
        (cond
         ((< n 100) (cons n c))
         (t (cons (mod n 100) (add (list 1) c))))))))

; Subtract big numbers a - b.  Assumes a >= b.

(defun sub (a b) (norm (sub* a b)))

(defun sub* (a b)
  (cond
   ((null b) a)
   (t (let ((n (- (car a) (car b)))
            (c (sub* (cdr a) (cdr b))))
        (cond
         ((>= n 0) (cons n c))
         (t (cons
             (+ n 100)
             (sub* c (list 1)))))))))

; Multiply two big numbers

(defun mul (a b)
  (cond
   ((or (null a) (null b)) nil)
   (t (add
       (bign (* (car a) (car b)))
       (cons 0
             (add
              (add
               (mul (list (car a)) (cdr b))
               (mul (cdr a) (list (car b))))
              (let ((c (mul (cdr a) (cdr b))))
                (when c (cons 0 c)))))))))

; Compare two big numbers and return :a>b, :a<b, or nil if they're equal.

(defun cmp (a b)
  (let ((la (length a)) (lb (length b)))
    (cond
     ((> la lb) :a>b)
     ((< la lb) :a<b)
     (t (cmp* a b)))))

(defun cmp* (a b)
  (cond
   ((null a) nil)
   ((cmp* (cdr a) (cdr b)))
   ((< (car a) (car b)) :a<b)
   ((> (car a) (car b)) :a>b)
   (t nil)))

; Divide bignum by 10

(defun div10 (a) (norm (div10* a)))

(defun div10* (a)
  (cond
   ((null a) nil)
   (t (let ((rest (div10* (cdr a))))
        (cons
         (+ (truncate (car a) 10) (* 10 (mod (or (cadr a) 0) 10)))
         rest)))))

; Divide a by b - returns a list of quotient and remainder

(defun divide (a b)
  (let ((cur (list 1))
        (den (copy-list b))
        (result nil)
        (remdr (copy-list a)))
    (loop
     (unless (eq (cmp a den) :a>b) (return))
     (push 0 den) (push 0 cur))
    (loop
     (when (null cur) (return))
     (dotimes (i 10)
       (when (eq (cmp remdr den) :a<b) 
         (setq result (add result (mul cur (bign i))))
         (return))
       (setq remdr (sub remdr den)))
     (setq den (div10 den)) (setq cur (div10 cur)))
    (list result remdr)))

(defun div (a b) (first (divide a b)))

(defun remdr (a b) (second (divide a b)))

; Applications

; Factorial n

(defun fac (n)
  (let ((f (big 1)))
    (dotimes (x n)
      (setq f (mul f (big (1+ x)))))
    f))

; Integer exponent - x^y

(defun iex (x y)
  (let ((e (big 1))
        (f (big x)))
    (loop
     (when (zerop y) (return))
     (when (oddp y) (setq e (mul e f)))
     (setq f (mul f f) y (ash y -1)))
    e))

; Factor - finds the smallest factor of a bignum n

(defun factor (n)
  (cond
   ((null (remdr n (big 2))) 2)
   (t 
    (let ((d (big 3)) (i (big 2)))
      (loop
       (when (eq (cmp (mul d d) n) :a>b) (return n))
       (when (null (remdr n d)) (return d))
       (setq d (add d i)))))))

; Factorize a bignum into a list of factors

(defun factorize (n)
  (let ((f (factor n)))
    (if (null (cmp n f))
        (list (pri f t)) 
      (cons (pri f t) (factorize (div n f))))))

; Calculate e to d digits

(defun euler (d)
  (let ((one (big 1))
        (term (iex 10 d))
        (denom nil)
        result)
    (loop
     (setq result (add result term))
     (setq denom (add denom one))
     (setq term (div term denom))
     (when (null term) (return result)))))

; Arctangent

(defun arctan (n x)
  (let* ((i 0) 
         (xn (div n x))
         (x2 (mul x x))
         (term xn)
         result)
    (loop
     (setq result 
           (if (evenp i) (add result term)
             (sub result term)))
     (incf i)
     (setq xn (div xn x2))
     (setq term (div xn (bign (1+ (* 2 i)))))
     (when (null term) (return result)))))

; Calculate pi to d digits using Machin's formula

(defun machin (d)
  (let ((n (iex 10 d)))
    (mul (big 4) (sub (mul (big 4) (arctan n (big 5))) (arctan n (big 239))))))