6
\$\begingroup\$

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))

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 the correct answer of (prime-factors 123456789) is (3 3 3607 3803), my program will take so long to find the next prime factor once it constructs the list with the two first elements (3 3). 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
\$\endgroup\$
1
  • 1
    \$\begingroup\$ please fix indentation (use Emacs if unsure). your code is now unreadable. \$\endgroup\$
    – sds
    Commented Mar 19, 2018 at 3:00

2 Answers 2

Reset to default
6
\$\begingroup\$

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))))))

PS. This is identical to https://stackoverflow.com/a/49369108/850781

\$\endgroup\$
2
\$\begingroup\$

Your is-prime is very suboptimal in few aspects.

I. There is no sense to check for prime factors of a number ABOVE the integral part of that number's square root.

II. Aside from 2, all other prime factors of a number are odd (if any), so you only need to check at most sqrt(n)/2 numbers.

III. And your recursion does not look tail-optimizable, so you'll need to either reshape it so that you get a single tail-call (if your CL implementation is aware of tail optimizations), or transform it into an iterative construct (which, I believe, are highly available in CL, much more than in Scheme).

After all, you don't need to obtain a list of ALL prime factors of a number to say it's composite; and even less do you need to keep a stack frame per EACH number less that it.

prime-factors-helper is less sensitive (all the recursive calls are either in a tail context, or in a tail modulo cons-one), but still you don't need to check every number less than x (remember, only up to sqrt(x) and only half of them).

Additionally, concerning (= expr 0), doesn't CL standard library feature a primitive predicate to compare against zero, like zero? in Scheme?

\$\endgroup\$
1
  • \$\begingroup\$ Yes, zerop as sds suggests. \$\endgroup\$
    – bipll
    Commented Mar 20, 2018 at 11:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.