(defun is-factor (x y)
  (zerop (mod x y)))

(defun is-congruent (num a n)
  (= (mod num n) a))

(defun is-prime (n)
  (cond ((< n 4) (or (= n 2) (= n 3)))
        ((or (zerop (mod n 2)) (zerop (mod n 3))) nil)
        (t (loop for i from 5 to (floor (sqrt n)) by 6
                 never (or (is-factor n i)
                           (is-factor n (+ i 2)))))))

(defun first-factor (n)
  (loop for i from 2 to (floor (sqrt n))
        thereis (if (is-factor n i) i nil)))

(defun is-blum (n)
  (let ((factor (first-factor n)))
    (and factor
         (is-congruent factor 3 4)
         (/= factor (/ n factor))
         (is-congruent (/ n factor) 3 4)
         (is-prime factor)
         (is-prime (/ n factor)))))

(defun main ()
  (let ((n 0)
        (i 1)
        (blums '())
        (counts '()))
    (dotimes (i 10) (push (cons i 0) counts))
    (loop while (<= n 400000)
          do (setf i (+ i 2))
             (when (is-blum i)
               (incf n)
               (incf (cdr (assoc (mod i 10) counts)))
               (when (<= n 50)
                 (push i blums)
                 (when (= n 50)
                   (format t "First 50 Blum integers:~%")
                   (format t "~{~5d~5d~5d~5d~5d~5d~5d~5d~5d~5d~%~}~%" (reverse blums))))
               (if (find n '(26828 100000 200000 300000 400000))
                   (format t "The ~7:dth Blum integer is: ~9:d~%" n i))))
    (format t "~%% distribution of the first 400,000 Blum integers:~%")
    (dolist (x '(1 3 7 9))
      (format t "~3$% end in ~a~%" (/ (cdr (assoc x counts)) 4000) x))))