How can I optimise my recursive Lisp function

Question

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

Answer 1

copied from https://codereview.stackexchange.com/a/189932/20936:

There are several problems with your code.

Style

1. is-prime is C/Java style. Lispers use primep or prime-number-p.
2. zerop is clearer than (= 0 ...).
3. 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))))))