Infinite precision arithmetic
11th March 2026
This is an updated version of the Infinite Precision Arithmetic package written in uLisp, adding support for division, and with several additional examples.
Introduction
One of the first example programs I wrote to test uLisp was the Infinite Precision Arithmetic package, which allowed you to calculate in uLisp with large numbers, beyond the precision of integers on the platform and limited only by the amount of workspace; hence "infinite". The original version dates back to 2016.
One of the aims of the package was to keep it simple and understandable, so users could extend it with additional functions. Also, it was designed to work on platforms with either 16-bit or 32-bit integers, allowing it to run on every platform, including AVR.
The original version only supported addition and multiplication. I've often thought of creating an update adding support for division, which would extend it to applications such as factorising large integers, or calculating the digits of π, and so here it is, together with some additional examples to demonstrate these extensions.
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 an integer or string into big number notation:
(defun big (x) (if (stringp x) (bigstr x) (bign x)))
If the argument is an integer it calls bign:
(defun bign (n)
(cond
((zerop n) nil)
((< n 100) (list n))
(t (cons
(mod n 100)
(bign (truncate n 100))))))
Here's an example:
> (big 1234567) (67 45 23 1)
A restriction of bign is that its argument can't be larger than the largest integer, which limits it to 2147483648 on 32-bit platforms. So bigstr allows you to create a big number from a string of any length:
(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:
> (big "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)))
(or
(format (not string) "~d~{~2,'0d~}~%" (car rn) (cdr rn))
nothing)))
The first pair of digits is printed with any leading zero suppressed. For example:
> (pri '(67 45 23 1)) 1234567
It returns nothing so the result isn't printed.
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
Here's the function sub to subtract big numbers:
(defun sub (a b) (norm (sub* a b)))
It uses a helper function, sub*, which 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 which doesn't affect its value, but can interfere with the operation of the other functions; for example:
> (sub* (big 123) (big 100)) (23 0)
The main sub function therefore removes trailing zeros with the norm (normalise) function:
(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 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 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)))
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)))))
f))
To calculate factorial 24:
> (pri (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)))
e))
For example to find 3100:
> (pri (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 returns the number itself.
Also, 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 as strings:
(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")
Calculating the digits of e
Euler's number e is transcendental, and can be calculated quite easily by the sum of an infinite series:
Here's the program:
(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)))))
This calculates e to d digits. For example:
> (pri (euler 30)) 2718281828459045235360287471339
Calculating π to d digits
Finally, here's a program that uses Machin's formula to calculate an arbitrary number of digits of π.
The formula is based on the following equation found by John Machin, an 18th century astronomer and mathematician:
If, like me, you're wondering where the values 5 and 239 come from, see the Wikipedia article [1].
The two arctangents can be calculated using the arctangent series:
Machin's formula converges extremely rapidly, making it a widely used method of calculating π to a large number of digits.
Here's the function arctan that iteratively calculates the arctangent of 1/x with a scale factor of n:
(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)))
(when (null term) (return result))
(incf i)
(setq xn (div xn x2))
(setq term (div xn (bign (1+ (* 2 i))))))))
So to calculate arctan 1/5 to 6 digits evaluate:
> (pri (arctan (big 1000000) (big 5))) 197397
Compare this with the calculation using floating-point arithmetic:
> (atan (/ 1 5)) 0.197396
Finally, the function machin calculates π 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))))))
Here's π to 30 digits:
> (pri (machin 30)) 3141592653589793238462643383280
Resources
Here's the updated version of the infinite precision package: Infinite precision package 2.
For the original article see: Infinite precision arithmetic.
- ^ Machin-like formula on Wikipedia.