Simple infinite precision arithmetic
7th March 2026
This is a simple infinite precision arithmetic package I wrote as an example program for my Lisp interpreter for small microcontrollers, uLisp.
First a couple of disclaimers:
- This could probably be written more efficiently. However, one of the aims of the package was to keep it concise and understandable, to show how easy it is to write this type of application in Lisp.
- Secondly, bignums are a standard part of Common Lisp, so on that platform there's no point in writing a package like this. However, it's still interesting to see how it would be done if you needed to add it to a platform that doesn't have it.
Representing big numbers
Each big number is represented as a list of integers in the range 0 to 99, with the least-significant digits at the start of the list, so the number 1234567 would be represented as the list:
(67 45 23 1)
The function big converts a uLisp integer into big number notation:
(defun big (n)
(cond
((< n 100) (list n))
(t (cons
(mod n 100)
(big (truncate n 100))))))
Here's an example:
> (big 1234567) (67 45 23 1)
This size was chosen so that even on a 16-bit platform, two integer elements can be multiplied without overflow.
Big number from string
A restriction of big is that its argument can't be larger than the largest integer, which limits it to 2147483648 on 32-bit platforms. So I've also provided bigstr that allows you to create a big number from a string:
(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))))
For example:
> (bigstr "9876543210") (10 32 54 76 98)
Printing big numbers
The function pri uses format to print each pair of digits in the list representing a number in big notation:
(defun pri (n &optional string)
(let ((rn (reverse n)))
(format (not string) "~d~{~2,'0d~}" (car rn) (cdr rn))))
The first pair of digits is printed with any leading zero suppressed. For example:
> (pri '(67 45 23 1)) 1234567 NIL
Setting the second argument to t returns the number as a string:
> (pri '(67 45 23 1) t) "1234567"
Adding big numbers
The function add adds together two numbers in big notation:
(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))))))))
It recursively adds the digits a pair at a time, starting with the low-order digits which are at the start of the two lists, taking account of carries. Here's the calculation of 654321 + 987654 = 1641975:
> (add '(21 43 65) '(54 76 98)) (75 19 64 1)
Subtracting big numbers
The main function, sub*, works in a similar way to add:
(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)))))))))
The subtraction can result in a trailing zero on the end of the result. This doesn't affect its value, but can interfere with the operation of the other functions, so we remove trailing zeros with the norm (normalise) function:
(defun sub (a b) (norm (sub* a b)))
Here's the definition of norm:
(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)))))))
In this recursive function rest contains the normalised tail. If rest is nil and the current element is zero then this zero must be part of the trailing zeros, so we drop it; otherwise we keep it.
For example:
> (norm '(67 45 23 1 0 0)) (67 45 23 1)
Multiplying big numbers
To multiply two big numbers we use the fact that:
(a + 100b) * (c + 100d) = ac + 100 * (bc + ad) + 10000 * bd.
To multiply a big number b by 100 we simply need to do:
(cons 0 b)
Here's the full definition of mul:
(defun mul (a b)
(cond
((or (null a) (null b)) nil)
(t (add
(big (* (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)))))))))
Here's the calculation of 654321 * 987654 = 646242752934:
> (mul '(21 43 65) '(54 76 98)) (34 29 75 42 62 64)
Compare
For the division operation we need to be able to compare two big numbers. The compare function, cmp, returns :a>b or :a<b if one is larger than the other, 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)))))
If one of the lists is longer than the other we can immediately assume that it's larger. Otherwise we need to compare the list elements individually; this is handled by the recursive subfunction cmp*:
(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)))
In this case we need to compare elements starting from the ends of the lists, since they are the most significant. We achieve this by calling cmp* recursively on the cdr of each list. If that returns non-nil it's the result to be returned; otherwise the cdrs are equal and we compare the car of each list.
Divide
Rather than having separate functions for divide and remainder it's convenient to make the main divide function return a list containing both the quotient and remainder.
Divide works like the method of long division as taught in school. The procedure is as follows:
- Set cur (current) to 1, den (the denominator) to a copy of the second argument, result to zero, and remdr (the remainder) to a copy of the numerator.
- Normalise the operands by shifting den and cur up, until den is larger than or equal to the numerator.
Then repeat the following steps until cur is zero:
- While remdr is greater than or equal to den, repeatedly subtract den from remdr and add cur to result.
- Then shift cur and den down.
This leaves the decision of how large to make the shift after each subtraction of den from remdr. In school long division the shift is one decimal digit, so you potentially have to do up to nine subtractions of den from remdr. An alternative would be to make the shift one bignum list element, which would be easy to implement, but you would then have to do up to 99 subtractions. A third alternative would be to make the shift one bit, so you have to do at most one subtraction.
I wrote all three versions and compared their performance on an ESP32-S3 using the test:
(time (factor (big 1111111)))
based on the factor example below:
| Shifts | Time |
| 2 | 3.8 s |
| 10 | 3.5 s |
| 100 | 9.4 s |
Surprisingly the version using up to nine subtractions was the fastest. Here's that version:
(defun divide (a b)
(let ((cur (list 1))
(den (copy-list b))
(result (list 0))
(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 (big i))))
(return))
(setq remdr (sub remdr den)))
(setq den (div10 den)) (setq cur (div10 cur)))
(list result remdr)))
It needs the function div10, which divides a 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)))))
For convenience, here are functions that get the two results from divide separately: div that divides two bignums, and remdr that gives the remainder:
(defun div (a b) (first (divide a b))) (defun remdr (a b) (second (divide a b)))
Applications
Here are some examples using the infinite arithmetic package:
Calculating factorials
This routine calculates the factorial of a number n iteratively:
(defun fac (n)
(let ((f (big 1)))
(dotimes (x n)
(setq f (mul f (big (1+ x)))))
(pri f)
(print nothing)))
To calculate factorial 24:
> (fac 24) 620448401733239439360000
Integer exponents
Here's a second example, which will find the value of xy for integer x and y. It uses binary exponentiation, which drastically reduces the number of multiplications needed:
(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)))
(pri e)
(print nothing)))
For example to find 3100:
> (iex 3 100) 515377520732011331036461129765621272702107522001
Factorising
Here's a simple approach to finding the least prime factor of a number by testing candidate factors up to the square root of the number:
(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)))))))
For example:
> (pri (factor (bigstr "134913016999"))) 2999
If the number is prime, factor prints the number itself.
Finally here's a function factorize that uses the above factor function as the basis for a simple recursive routine to factorize a number into a list of its prime 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))))))
For example:
> (factorize (bigstr "84061014001"))
("3001" "4001" "7001")
Here's the whole infinite precision package: Infinite precision package.
blog comments powered by Disqus