fast Project-euler 7

;; Using the lambert w-function left-branch in order to approximate a pi-¹(x)
;; function which allows a sieve to be used to generate a target number of primes.
;; Since the pi-¹(x) function overestimates, it can be guaranteed to always generate enough
(defmacro while (test &body body)
  `(do ()
       ((not ,test))
     ,@body))

(defmacro with-gensyms (syms &body body)
  `(let ,(mapcar #'(lambda (x) `(,x (gensym))) syms)
     ,@body))


(defmacro collect-list (end-condition next-value var)
  "Generates a list based on the next-value provided."
  (with-gensyms (acc)
    `(let ((,acc nil))
       (while (not ,end-condition)
         (push ,var ,acc)
         (setq ,var ,next-value))
       (nreverse (cons ,var ,acc)))))

(defun divisible (n x)
  (zerop (rem n x)))

(defun take (n xs)
  (labels ((rec (n xs acc)
             (cond ((= n 0) (nreverse acc))
                   ((null xs) nil)
                   (t (rec (1- n) (cdr xs) (cons (car xs) acc))))))
    (rec n xs nil)))

(defun range (x y &optional (step 1))
  "Generates a range with an optional step parameter. It is non-inclusive."
    (if (> (abs (- y (+ step x))) (abs (- y x)))
          nil
              (nbutlast (collect-list (= x y) (+ x step) x))))
(defun square (x) (* x x))
(defun primes<n (n)
  (let ((upb (isqrt n)))
    (labels ((rec (p xs acc)
               (if (> p  upb)
                   (nconc (nreverse acc) xs)
                   (let ((nl (remove-if (lambda (x) (divisible x p)) xs)))
                     (rec (car nl) nl (cons p acc))))))
      (rec 2 (range 2 (1+ n)) nil))))
(defun W-1 (z)
  "SEE 'Numerical evaluation of the lambert w function and application
 to generation of generalized gaussian noise with exponent 1/2'"
  (cond
    ((and (< z -0.333) (>= z (/ -1 (exp 1))))
     (let ((p (- (sqrt (* 2 (1+ (* (exp 1) z)))))))
       (+
        -1
        p
        (- (* 1/3 (expt p 2)))
        (* 11/72 (expt p 3))
        (- (* 43/540 (expt p 4)))
        (* 769/17280 (expt p 5))
        (- (* 221/8505 (expt p 6))))))
    ((and (<= -0.333 z) (<= z -0.033))
     (/
      (+
       -8.0960
       (* 391.0025 z)
       (- (* 47.4252 (square z)))
       (- (* 4877.6330 (expt z 3)))
       (- (* 5532.7760 (expt z 4))))
      (+
       1
       (- (* 82.9423 z))
       (* 433.8688 (square z))
       (* 1515.3060 (expt z 3)))))
    ((and (< -0.033 z) (< z 0))
     (let ((l1 (log (- z))) (l2 (log (- (log (- z))))))
       (+
        l1
        (- l2)
        (/ l2 l1)
        (/ (* (+ -2 l2) l2)
           (* 2 (square l1)))
        (/ (* l2 (+ 6 (* -9 l2) (* 2 (square l2))))
           (* 6 (expt l1 3)))
        (/ (* l2 (+ -12 (* 36 l2) (* -22 (square l2)) (* 3 (expt l2 3))))
           (* 12 (expt l1 4)))
        (/ (* l2 (+ 60 (* -300 l2) (* 350 (square l2)) (* -125 (expt l2 3)) (* 12 (expt l2 4))))
           (* 60 (expt l1 5))))))))

(defun pi-1 (n)
  (round (exp (- (w-1 (/ -1 n))))))

(defun n-primes (n)
  (if (< n 15)
      (take n '(2 3 5 7 11 13 17 19 23 29 31 37 41 47))
      (let ((pl (primes<n (pi-1 n))))
        (take n pl))))

(princ (car (last (n-primes 10001))))
(princ "
")