オイラーの贈物とLisp(1.1.1)

オイラーの贈物 を読みながら、登場する数式などをCommon Lispで書いてみます。

1.1.1 自然数と素数(P6)より、有名なエラトステネスの篩と、 Wikipediaの 素数の項目 に載っていたウラムの螺旋を出力するコードです。 画像の作成には Vecto を利用しました。

(defun eratosthenes-sieve (n)
  "Return a sequence which indicates whether index is the prime number or not."
  ;; `n' is greater than or equal to 2
  (check-type n (integer 2 *))
  (let ((seq (make-array (1+ n) :initial-element t)))
    (setf (svref seq 0) nil) ; 0 is not a prime number.
    (setf (svref seq 1) nil) ; 1 is not a prime number.
    ;; The outer loop can stop at the square root of `n'.
    (loop :for i :from 2 :to (floor (sqrt n))
       :when (svref seq i)
       :do
       ;; The inner loop can start at the square of `i'.
       ;; (Multiples of `i' which are less than the square of `i' are already set to `nil'.)
       (loop :for j :from (* i 2) :to n :by i
          :do (setf (svref seq j) nil)))
    seq))


(ql:quickload "vecto")

(defun draw-uram-spiral (edge-length output-file-name &key (pixel 2))
  (assert (<= 1 edge-length))
  (let* ((limit (expt edge-length 2))
         (primes (eratosthenes-sieve (expt edge-length 2)))
         (picture-edge-length (* pixel edge-length)))
    (vecto:with-canvas
        (:width picture-edge-length :height picture-edge-length)
      (vecto:set-rgb-fill 1.0 1.0 1.0)
      (vecto:rectangle 0.0 0.0 picture-edge-length picture-edge-length)
      (vecto:fill-path)
      (vecto:set-rgb-fill 0.0 0.0 1.0)
      (let ((idx 2)
            (step 0)
            (x (floor edge-length 2))
            (y (floor edge-length 2)))
        (loop
           :while (<= idx limit)
           :do 
           (dotimes (_ (1+ (floor step 2)))
             (when (<= idx limit)
               (case (mod step 4)
                 (0 (incf x))
                 (1 (decf y))
                 (2 (decf x))
                 (3 (incf y)))
               (when (svref primes idx)
                 (vecto:rectangle (* x pixel) (* y pixel) pixel pixel)
                 (vecto:fill-path))
               (incf idx)))
           (incf step)))
      (vecto:save-png output-file-name))))

> (draw-uram-spiral 100 "uram.png")