;; -*- coding: utf-8 -*-
(eval-when (:compile-toplevel :load-toplevel :execute)
  (defparameter OPT
    #+swank '(optimize (speed 3) (safety 2))
    #-swank '(optimize (speed 3) (safety 0) (debug 0)))
  #+swank (progn (ql:quickload '(:cl-debug-print :fiveam))
                 (shadow :run)
                 (use-package :fiveam))
  #-swank (set-dispatch-macro-character #\# #\> (lambda (s c p) (declare (ignore c p)) (read s nil (values) t))))
#+swank (cl-syntax:use-syntax cl-debug-print:debug-print-syntax)

;; BEGIN_INSERTED_CONTENTS
;;;
;;; Real FFT
;;;
;;; Reference:
;;; http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn2_12.html#sec2_1_2
;;;

(deftype fft-float () 'single-float)

(declaim (inline power2-p))
(defun power2-p (x)
  "Checks if X is a power of 2."
  (zerop (logand x (- x 1))))

;; For FFT of fixed length, preparing the table of cos(i*theta) and sin
;; (i*theta) will be efficient.
(defun %make-trifunc-table (n)
  (declare (optimize (speed 3))
           ((integer 0 #.most-positive-fixnum) n))
  (assert (power2-p n))
  (let* ((cos-table (make-array (ash n -2) :element-type 'fft-float))
         (sin-table (make-array (ash n -2) :element-type 'fft-float))
         (theta (/ (coerce (* 2 pi) 'fft-float) n)))
    (dotimes (i (ash n -2))
      (setf (aref cos-table i) (cos (* i theta))
            (aref sin-table i) (sin (* i theta))))
    (values cos-table sin-table)))

(defparameter *cos-table* nil)
(defparameter *sin-table* nil)

(defmacro with-fixed-base (size &body body)
  "Makes FFT faster when the SIZE of target vectors is fixed in BODY. This macro
computes and holds the roots of unity for SIZE, which DFT! and INVERSE-DFT!
called in BODY automatically detects; they will signal an error when they
receive a vector of different size."
  (let ((s (gensym)))
    `(let ((,s ,size))
       (multiple-value-bind (*cos-table* *sin-table*) (%make-trifunc-table ,s)
         ,@body))))

(defun %dft-fixed-base! (f)
  (declare #.OPT
           ((simple-array fft-float (*)) f))
  (prog1 f
    (let* ((n (length f))
           (cos-table *cos-table*)
           (sin-table *sin-table*)
           (factor n))
      (declare ((integer 0 #.most-positive-fixnum) factor)
               ((simple-array fft-float (*)) cos-table sin-table))
      (assert (power2-p n))
      (assert (= (ash n -2) (length cos-table)))
      ;; bit-reverse ordering
      (let ((i 0))
        (declare ((integer 0 #.most-positive-fixnum) i))
        (loop for j from 1 below (- n 1)
              do (loop for k of-type (integer 0 #.most-positive-fixnum)
                          = (ash n -1) then (ash k -1)
                       while (> k (setq i (logxor i k))))
                 (when (< j i)
                   (rotatef (aref f i) (aref f j)))))
      (do* ((mh 1 m)
            (m (ash mh 1) (ash mh 1)))
           ((> m n))
        (declare ((integer 0 #.most-positive-fixnum) mh m))
        (let ((mq (ash mh -1)))
          (setq factor (ash factor -1))
          (do ((jr 0 (+ jr m)))
              ((>= jr n))
            (declare ((integer 0 #.most-positive-fixnum) jr))
            (let ((xreal (aref f (+ jr mh))))
              (setf (aref f (+ jr mh)) (- (aref f jr) xreal))
              (incf (aref f jr) xreal)))
          (do ((i 1 (+ i 1)))
              ((>= i mq))
            (declare ((integer 0 #.most-positive-fixnum) i))
            (let* ((index (the fixnum (* factor i)))
                   (wreal (aref cos-table index))
                   (wimag (- (aref sin-table index))))
              (do ((j 0 (+ j m)))
                  ((>= j n))
                (let* ((j+mh (+ j mh))
                       (j+m-i (- (+ j m) i))
                       (xreal (+ (* wreal (aref f (+ j+mh i)))
                                 (* wimag (aref f j+m-i))))
                       (ximag (- (* wreal (aref f j+m-i))
                                 (* wimag (aref f (+ j+mh i))))))
                  (declare ((integer 0 #.most-positive-fixnum) j+mh j+m-i))
                  (setf (aref f (+ j+mh i))
                        (+ (- (aref f (- j+mh i))) ximag))
                  (setf (aref f j+m-i)
                        (+ (aref f (- j+mh i)) ximag))
                  (setf (aref f (- j+mh i))
                        (+ (aref f (+ j i)) (- xreal)))
                  (incf (aref f (+ j i)) xreal))))))))))

(defun %inverse-dft-fixed-base! (f)
  (declare #.OPT
           ((simple-array fft-float (*)) f))
  (prog1 f
    (let* ((n (length f))
           (cos-table *cos-table*)
           (sin-table *sin-table*)
           (factor 1))
      (declare ((integer 0 #.most-positive-fixnum) factor)
               ((simple-array fft-float (*)) cos-table sin-table))
      (assert (power2-p n))
      (assert (= (ash n -2) (length cos-table)))
      (setf (aref f 0)
            (/ (aref f 0) 2))
      (setf (aref f (ash n -1))
            (/ (aref f (ash n -1)) 2))
      (do* ((m n mh)
            (mh (ash m -1) (ash m -1)))
           ((zerop mh))
        (declare ((integer 0 #.most-positive-fixnum) m mh))
        (let ((mq (ash mh -1)))
          (do ((jr 0 (+ jr m)))
              ((>= jr n))
            (declare ((integer 0 #.most-positive-fixnum) jr))
            (let ((xreal (- (aref f jr) (aref f (+ jr mh)))))
              (incf (aref f jr) (aref f (+ jr mh)))
              (setf (aref f (+ jr mh)) xreal)))
          (do ((i 1 (+ i 1)))
              ((>= i mq))
            (let* ((index (the fixnum (* factor i)))
                   (wreal (aref cos-table index))
                   (wimag (aref sin-table index)))
              (do ((j 0 (+ j m)))
                  ((>= j n))
                (let* ((j+mh (+ j mh))
                       (j+m-i (- (+ j m) i))
                       (xreal (- (aref f (+ j i)) (aref f (- j+mh i))))
                       (ximag (+ (aref f j+m-i) (aref f (+ j+mh i)))))
                  (declare ((integer 0 #.most-positive-fixnum) j+mh j+m-i))
                  (incf (aref f (+ j i)) (aref f (- j+mh i)))
                  (setf (aref f (- j+mh i))
                        (- (aref f j+m-i) (aref f (+ j+mh i))))
                  (setf (aref f (+ j+mh i))
                        (+ (* wreal xreal) (* wimag ximag)))
                  (setf (aref f j+m-i)
                        (- (* wreal ximag) (* wimag xreal))))))))
        (setq factor (ash factor 1)))
      ;; bit-reverse ordering
      (let ((i 0))
        (declare ((integer 0 #.most-positive-fixnum) i))
        (loop for j from 1 below (- n 1)
              do (loop for k of-type (integer 0 #.most-positive-fixnum)
                          = (ash n -1) then (ash k -1)
                       while (> k (setq i (logxor i k))))
                 (when (< j i)
                   (rotatef (aref f i) (aref f j))))))))

(declaim (inline dft!))
(defun dft! (f)
  (declare ((simple-array fft-float (*)) f))
  (if (zerop (length f))
      f
      (%dft-fixed-base! f)))

(declaim (inline inverse-dft!))
(defun inverse-dft! (f)
  (declare ((simple-array fft-float (*)) f))
  (prog1 f
    (let ((n (length f)))
      (unless (zerop n)
        (let ((factor (* 2 (/ (coerce n 'fft-float)))))
          (%inverse-dft-fixed-base! f)
          (dotimes (i n)
            (setf (aref f i) (* (aref f i) factor))))))))

(defmacro dbg (&rest forms)
  #+swank
  (if (= (length forms) 1)
      `(format *error-output* "~A => ~A~%" ',(car forms) ,(car forms))
      `(format *error-output* "~A => ~A~%" ',forms `(,,@forms)))
  #-swank (declare (ignore forms)))

(declaim (inline println))
(defun println (obj &optional (stream *standard-output*))
  (let ((*read-default-float-format* 'double-float))
    (prog1 (princ obj stream) (terpri stream))))

;; Body
(defconstant +dp-size+ (ash 1 18))

(defun main ()
  (declare #.OPT)
  (let* ((k (read))
         (dp (make-array +dp-size+ :element-type 'fft-float :initial-element 0f0))
         (multiplier (make-array +dp-size+ :element-type 'fft-float :initial-element 0f0)))
    (declare ((integer 0 100000) k)
             ((simple-array fft-float (#.+dp-size+)) dp multiplier))
    (loop for i = 1 then (rem (* i 10) k)
          while (= 0f0 (aref dp i))
          do (setf (aref dp i) 1f0
                   (aref multiplier i) 1f0))
    (with-fixed-base +dp-size+
      (dft! multiplier)
      (loop for i from 1 to 45
            do (when (= 1f0 (aref dp 0))
                 (println i)
                 (return-from main))
               (dft! dp)
               (setf (aref dp 0)
                     (* (aref dp 0) (aref multiplier 0)))
               (setf (aref dp (ash +dp-size+ -1))
                     (* (aref dp (ash +dp-size+ -1)) (aref multiplier (ash +dp-size+ -1))))
               (loop for i from 1 below (ash +dp-size+ -1)
                     for value1 = (- (* (aref dp i) (aref multiplier i))
                                     (* (aref dp (- +dp-size+ i)) (aref multiplier (- +dp-size+ i))))
                     for value2 = (+ (* (aref dp i) (aref multiplier (- +dp-size+ i)))
                                     (* (aref dp (- +dp-size+ i)) (aref multiplier i)))
                     do (setf (aref dp i) value1
                              (aref dp (- +dp-size+ i)) value2))
               (inverse-dft! dp)
               (dotimes (i k)
                 (if (or (> (aref dp i) 0.5f0)
                         (> (aref dp (+ i k)) 0.5f0))
                     (setf (aref dp i) 1f0)
                     (setf (aref dp i) 0f0))
                 (setf (aref dp (+ i k)) 0f0))
               (loop for i from (* 2 k) below +dp-size+
                     do (setf (aref dp i) 0f0))))))

#-swank(main)


Measured memory_bytes: 27988