;; -*- 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
;;;
;;; In-place FFT
;;;
;;; Reference:
;;; http://www.prefield.com/algorithm/math/fft.html
;;; http://techtipshoge.blogspot.com/2011/08/fft4.html
;;;

(defconstant +dp-size+ (ash 1 18))
(deftype fft-float () 'single-float)

;; For FFT of fixed length, it is faster preparing the table of exp(i*theta).
(defun %make-cis-table (sign)
  (declare #.OPT
           ((integer -1 1) sign))
  (let* ((table (make-array #.(ash +dp-size+ -1) :element-type '(complex fft-float)))
         (theta (* sign #.(/ (coerce (* 2 pi) 'fft-float) +dp-size+))))
    (dotimes (i #.(ash +dp-size+ -1))
      (setf (aref table i) (cis (* i theta))))
    table))

(defconstant +cis+table+ (%make-cis-table 1))
(defconstant +cis-table+ (%make-cis-table -1))

(defun %general-dft-cached-cis! (f sign)
  (declare (optimize (speed 3) (safety 0) (debug 0))
           ((simple-array (complex fft-float) (#.+dp-size+)) f)
           ((integer -1 1) sign))
  (prog1 f
    (let ((table (if (= 1 sign) +cis+table+ +cis-table+)))
      (declare ((simple-array (complex fft-float) (#.(ash +dp-size+ -1))) table))
      (do ((m +dp-size+ (ash m -1))
           (shift 0 (+ shift 1)))
          ((= m 1))
        (declare ((integer 0 #.most-positive-fixnum) m shift))
        (let ((mh (ash m -1)))
          (dotimes (i mh)
            (do ((j i (+ j m)))
                ((>= j +dp-size+))
              (declare ((integer 0 #.most-positive-fixnum) j))
              (let* ((k (+ j mh))
                     (xt (- (aref f j) (aref f k)))
                     (cis-index (ash i shift)))
                (declare ((integer 0 #.most-positive-fixnum) k cis-index))
                (incf (aref f j) (aref f k))
                (setf (aref f k) (* (aref table cis-index) xt)))))))
      (let ((i 0))
        (declare ((integer 0 #.most-positive-fixnum) i))
        (loop for j from 1 below #.(- +dp-size+ 1)
              do (loop for k of-type (integer 0 #.most-positive-fixnum)
                          = #.(ash +dp-size+ -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 #.OPT
           ((simple-array (complex fft-float) (#.+dp-size+)) f))
  (if (zerop (length f))
      f
      (%general-dft-cached-cis! f 1)))

(declaim (inline inverse-dft!))
(defun inverse-dft! (f)
  (declare #.OPT
           ((simple-array (complex fft-float) (#.+dp-size+)) f))
  (prog1 f
    (let ((n (length f)))
      (unless (zerop n)
        (let ((/n (/ (coerce n 'fft-float))))
          (%general-dft-cached-cis! f -1)
          (dotimes (i n)
            (setf (aref f i) (* (aref f i) /n))))))))

(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

(defun main ()
  (declare #.OPT)
  (let* ((k (read))
         (dp (make-array +dp-size+ :element-type '(complex fft-float) :initial-element #c(0f0 0f0)))
         (multiplier (make-array +dp-size+ :element-type '(complex fft-float) :initial-element #c(0f0 0f0))))
    (declare ((integer 0 100000) k)
             ((simple-array (complex fft-float) (#.+dp-size+)) dp multiplier))
    (loop for i = 1 then (rem (* i 10) k)
          while (= #c(0f0 0f0) (aref dp i))
          do (setf (aref dp i) #c(1f0 0f0)
                   (aref multiplier i) #c(1f0 0f0)))
    (dft! multiplier)
    (loop for i from 1 to 45
          do (unless (= #c(0f0 0f0) (aref dp 0))
               (println i)
               (return-from main))
             (when (>= i 35)
               (println (let ((res 0))
                          (declare (fixnum res))
                          (loop (when (zerop k) (return res))
                                (multiple-value-bind (quot rem) (floor k 10)
                                  (incf res rem)
                                  (setq k quot)))))
               (return-from main))
             (dft! dp)
             (dotimes (i +dp-size+)
               (setf (aref dp i) (* (aref dp i) (aref multiplier i))))
             (inverse-dft! dp)
             (dotimes (i k)
               (if (or (> (realpart (aref dp i)) 0.5f0)
                       (> (realpart (aref dp (+ i k))) 0.5f0))
                   (setf (aref dp i) #c(1f0 0f0))
                   (setf (aref dp i) #c(0f0 0f0)))
               (setf (aref dp (+ i k)) #c(0f0 0f0)))
             (loop for i from (* 2 k) below +dp-size+
                   do (setf (aref dp i) #c(0f0 0f0))))))

#-swank(main)


Measured memory_bytes: 31328