I am trying to create a function prime-factors
that returns the prime factors of a number. To do so, I created is-prime
function, and prime-factors-helper
that will do a recursive check of the prime factors.
(defun is-prime (n &optional (d (- n 1)))
(if (/= n 1) (or (= d 1)
(and (/= (rem n d) 0)
(is-prime n (- d 1)))) ()))
(defun prime-factors-helper (x n)
(if (is-prime x)
(list x)
(if (is-prime n)
(if (AND (= (mod x n) 0) (<= n (/ x 2)))
(cons n (prime-factors-helper (/ x n) n))
(prime-factors-helper x (+ 1 n)))
(prime-factors-helper x (+ 1 n)))))
(defun prime-factors (x)
(prime-factors-helper x 2))
Question
I have a problem of optimisation. When I have a big number such as 123456789
, I get this error message Stack overflow (stack size 261120)
. I believe because since the correct answer is (3 3 3607 3803)
, my program once it constructs the list with the two first elements (3 3)
, it will take so long to find the next prime factor. How can I optimise my code?
CL-USER 53 > (prime-factors 512)
(2 2 2 2 2 2 2 2 2)
CL-USER 54 > (prime-factors 123456789)
Stack overflow (stack size 261120).
1 (abort) Return to level 0.
2 Return to top loop level 0.
Type :b for backtrace or :c <option number> to proceed.
Type :bug-form "<subject>" for a bug report template or :? for other options
copied from https://codereview.stackexchange.com/a/189932/20936:
There are several problems with your code.
Style
is-prime
is C/Java style. Lispers use primep
or prime-number-p
.
zerop
is clearer than (= 0 ...)
.
- Lispers use indentation, not paren counting, to read code. Your
code is thus virtually unreadable. Please use Emacs if you are
unsure how to format lisp properly.
Stack overflow
is-prime
is tail-recursive, so if you compile it, it should become a
simple loop and there should be no stack issues.
However, do not rush with it yet.
Algorithm
Number of iterations
Let us use trace
to see the
problems:
> (prime-factors 17)
1. Trace: (IS-PRIME '17)
2. Trace: (IS-PRIME '17 '15)
3. Trace: (IS-PRIME '17 '14)
4. Trace: (IS-PRIME '17 '13)
5. Trace: (IS-PRIME '17 '12)
6. Trace: (IS-PRIME '17 '11)
7. Trace: (IS-PRIME '17 '10)
8. Trace: (IS-PRIME '17 '9)
9. Trace: (IS-PRIME '17 '8)
10. Trace: (IS-PRIME '17 '7)
11. Trace: (IS-PRIME '17 '6)
12. Trace: (IS-PRIME '17 '5)
13. Trace: (IS-PRIME '17 '4)
14. Trace: (IS-PRIME '17 '3)
15. Trace: (IS-PRIME '17 '2)
16. Trace: (IS-PRIME '17 '1)
16. Trace: IS-PRIME ==> T
15. Trace: IS-PRIME ==> T
14. Trace: IS-PRIME ==> T
13. Trace: IS-PRIME ==> T
12. Trace: IS-PRIME ==> T
11. Trace: IS-PRIME ==> T
10. Trace: IS-PRIME ==> T
9. Trace: IS-PRIME ==> T
8. Trace: IS-PRIME ==> T
7. Trace: IS-PRIME ==> T
6. Trace: IS-PRIME ==> T
5. Trace: IS-PRIME ==> T
4. Trace: IS-PRIME ==> T
3. Trace: IS-PRIME ==> T
2. Trace: IS-PRIME ==> T
1. Trace: IS-PRIME ==> T
(17)
You do 17 iterations when only (isqrt 17) = 4
iterations are necessary.
Recalculations
Now compile is-prime
to turn recursion into a loop and see:
> (prime-factors 12345)
1. Trace: (IS-PRIME '12345)
1. Trace: IS-PRIME ==> NIL
1. Trace: (IS-PRIME '2)
1. Trace: IS-PRIME ==> T
1. Trace: (IS-PRIME '12345)
1. Trace: IS-PRIME ==> NIL
1. Trace: (IS-PRIME '3)
1. Trace: IS-PRIME ==> T
1. Trace: (IS-PRIME '4115)
1. Trace: IS-PRIME ==> NIL
1. Trace: (IS-PRIME '3)
1. Trace: IS-PRIME ==> T
1. Trace: (IS-PRIME '4115)
1. Trace: IS-PRIME ==> NIL
1. Trace: (IS-PRIME '4)
1. Trace: IS-PRIME ==> NIL
1. Trace: (IS-PRIME '4115)
1. Trace: IS-PRIME ==> NIL
1. Trace: (IS-PRIME '5)
1. Trace: IS-PRIME ==> T
1. Trace: (IS-PRIME '823)
1. Trace: IS-PRIME ==> T
(3 5 823)
You are checking the primality of the same numbers several times!
Extra optimization
primep
can find a divisor, not just check primality.
Optimized algorithm
(defun compositep (n &optional (d (isqrt n)))
"If n is composite, return a divisor.
Assumes n is not divisible by anything over d."
(and (> n 1)
(> d 1)
(if (zerop (rem n d))
d
(compositep n (- d 1)))))
(defun prime-decomposition (n)
"Return the prime decomposition of n."
(let ((f (compositep n)))
(if f
(nconc (prime-decomposition (/ n f))
(prime-decomposition f))
(list n))))
Note that one final optimization is possible -
memoization of
compositep
:
(let ((known-composites (make-hash-table)))
(defun compositep (n &optional (d (isqrt n)))
"If n is composite, return a divisor.
Assumes n is not divisible by anything over d."
(multiple-value-bind (value found-p) (gethash n known-composites)
(if found-p
value
(setf (gethash n known-composites)
(and (> n 1)
(> d 1)
(if (zerop (rem n d))
d
(compositep n (- d 1)))))))))
or, better yet, of prime-decomposition
:
(let ((known-decompositions (make-hash-table)))
(defun prime-decomposition (n)
"Return the prime decomposition of n."
(or (gethash n known-decompositions)
(setf (gethash n known-decompositions)
(let ((f (compositep n)))
(if f
(append (prime-decomposition (/ n f))
(prime-decomposition f))
(list n)))))))
note the use or append
instead of
nconc
.
Another interesting optimization is changing the iteration in
compositep
from descending to ascending.
This should speedup it up considerably as it would terminate early more
often:
(let ((known-composites (make-hash-table)))
(defun compositep (n)
"If n is composite, return a divisor.
Assumes n is not divisible by anything over d."
(multiple-value-bind (value found-p) (gethash n known-composites)
(if found-p
value
(setf (gethash n known-composites)
(loop for d from 2 to (isqrt n)
when (zerop (rem n d))
return d))))))