(in-package :cl-user)
(eval-when (:compile-toplevel :load-toplevel :execute)
  (sb-int:defconstant-eqx opt
    #+swank '(optimize (speed 3) (safety 2))
    #-swank '(optimize (speed 3) (safety 0) (debug 0))
    #'equal)
  #+swank (ql:quickload '(:cl-debug-print :fiveam) :silent t)
  #-swank (set-dispatch-macro-character
           #\# #\> (lambda (s c p) (declare (ignore c p)) `(values ,(read s nil nil t)))))
#+swank (cl-syntax:use-syntax cl-debug-print:debug-print-syntax)

(defmacro define-int-types (&rest bits)
  `(progn
     ,@(mapcar (lambda (b) `(deftype ,(intern (format nil "UINT~A" b)) () '(unsigned-byte ,b))) bits)
     ,@(mapcar (lambda (b) `(deftype ,(intern (format nil "INT~A" b)) () '(signed-byte ,b))) bits)))
(define-int-types 2 4 7 8 15 16 31 32 62 63 64)

(defconstant +mod+ 200003)
(defconstant +generator+ 2)

;; BEGIN_INSERTED_CONTENTS
(defpackage :cp/read-fixnum
  (:use :cl)
  (:export #:read-fixnum))
(in-package :cp/read-fixnum)

(declaim (ftype (function * (values fixnum &optional)) read-fixnum))
(defun read-fixnum (&optional (in *standard-input*))
  "NOTE: cannot read -2^62"
  (declare #.cl-user::opt)
  (macrolet ((%read-byte ()
               `(the (unsigned-byte 8)
                     #+swank (char-code (read-char in nil #\Nul))
                     #-swank (sb-impl::ansi-stream-read-byte in nil #.(char-code #\Nul) nil))))
    (let* ((minus nil)
           (result (loop (let ((byte (%read-byte)))
                           (cond ((<= 48 byte 57)
                                  (return (- byte 48)))
                                 ((zerop byte) ; #\Nul
                                  (error "Read EOF or #\Nul."))
                                 ((= byte #.(char-code #\-))
                                  (setq minus t)))))))
      (declare ((integer 0 #.most-positive-fixnum) result))
      (loop
        (let* ((byte (%read-byte)))
          (if (<= 48 byte 57)
              (setq result (+ (- byte 48)
                              (* 10 (the (integer 0 #.(floor most-positive-fixnum 10))
                                         result))))
              (return (if minus (- result) result))))))))

(defpackage :cp/mod-power-table
  (:use :cl)
  (:export #:make-mod-power-table))
(in-package :cp/mod-power-table)

(declaim (inline make-mod-power-table))
(defun make-mod-power-table (base length modulus &optional (element-type '(unsigned-byte 31)))
  "Returns a vector of the given length: VECTOR[x] := BASE^x mod MODULUS."
  (declare (fixnum base)
           ((integer 0 #.most-positive-fixnum) length)
           ((integer 1 #.most-positive-fixnum) modulus))
  (let ((res (make-array length :element-type element-type)))
    (unless (zerop length)
      (setf (aref res 0) 1)
      (loop for i from 1 below length
            do (setf (aref res i)
                     (mod (* base (aref res (- i 1))) modulus))))
    res))

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

(defpackage :cp/fft-real
  (:use :cl)
  (:export #:fft-float #:with-fixed-length-fft #:dft! #:inverse-dft! #:convolute!))
(in-package :cp/fft-real)

(deftype fft-float () 'double-float)

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

(defun %dft! (f)
  (declare (optimize (speed 3) (safety 0))
           ((simple-array fft-float (*)) f))
  (prog1 f
    (let ((n (length f))
          (theta (coerce (* -2 pi) 'fft-float)))
      (declare (fft-float theta))
      (assert (power2-p n))
      ;; 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 theta (/ theta 2))
          (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 ((wreal (cos (* theta i)))
                  (wimag (sin (* theta i))))
              (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! (f)
  (declare (optimize (speed 3) (safety 0))
           ((simple-array fft-float (*)) f))
  (prog1 f
    (let* ((n (length f))
           (theta (/ (coerce (* 2 pi) 'fft-float) n)))
      (declare (fft-float theta))
      (assert (power2-p n))
      (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 ((wreal (cos (* theta i)))
                  (wimag (sin (* theta i))))
              (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 theta (* theta 2)))
      ;; 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))))))))

;; 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) (safety 0)))
  (check-type n (integer 0 #.most-positive-fixnum))
  (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-length-fft (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 (optimize (speed 3) (safety 0))
           ((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 (optimize (speed 3) (safety 0))
           ((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
      (if *cos-table*
          (%dft-fixed-base! f)
          (%dft! 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)))))
          (if *cos-table*
              (%inverse-dft-fixed-base! f)
              (%inverse-dft! f))
          (dotimes (i n)
            (setf (aref f i) (* (aref f i) factor))))))))

(declaim (inline convolute!))
(defun convolute! (f1 f2 &optional result-vector)
  "Returns the convolution of two vectors F1 and F2. A new vector is created when
RESULT-VECTOR is null. This function destructively modifies F1 and F2. (They can
be restored by INVERSE-DFT!.)"
  (declare ((simple-array fft-float (*)) f1 f2)
           ((or null (simple-array fft-float (*))) result-vector))
  (let ((n (length f1)))
    (assert (and (power2-p n)
                 (= n (length f2))))
    (dft! f1)
    (dft! f2)
    (let ((f (or result-vector (make-array n :element-type 'fft-float))))
      (unless (zerop n)
        (setf (aref f 0)
              (* (aref f1 0) (aref f2 0)))
        (setf (aref f (ash n -1))
              (* (aref f1 (ash n -1)) (aref f2 (ash n -1)))))
      (loop for i from 1 below (ash n -1)
            for value1 of-type fft-float
               = (- (* (aref f1 i) (aref f2 i))
                    (* (aref f1 (- n i)) (aref f2 (- n i))))
            for value2 of-type fft-float
               = (+ (* (aref f1 i) (aref f2 (- n i)))
                    (* (aref f1 (- n i)) (aref f2 i)))
            do (setf (aref f i) value1)
               (setf (aref f (- n i)) value2))
      (inverse-dft! f))))

;; BEGIN_USE_PACKAGE
(eval-when (:compile-toplevel :load-toplevel :execute)
  (use-package :cp/fft-real :cl-user)
  (use-package :cp/mod-power-table :cl-user)
  (use-package :cp/read-fixnum :cl-user))
(in-package :cl-user)

(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
;;;

(declaim (ftype (function * (values (simple-array fft-float (*)) &optional)) fft-square!))
(defun fft-square! (vector)
  (declare #.opt
           ((simple-array fft-float (*)) vector))
  (let* ((n (length vector)))
    (dft! vector)
    (let ((res (make-array n :element-type 'fft-float)))
      (unless (zerop n)
        (setf (aref res 0)
              (expt (aref vector 0) 2))
        (setf (aref res (ash n -1))
              (expt (aref vector (ash n -1)) 2)))
      (loop for i from 1 below (ash n -1)
            for value1 of-type fft-float
               = (- (expt (aref vector i) 2)
                    (expt (aref vector (- n i)) 2))
            for value2 of-type fft-float
               = (+ (* (aref vector i) (aref vector (- n i)))
                    (* (aref vector (- n i)) (aref vector i)))
            do (setf (aref res i) value1)
               (setf (aref res (- n i)) value2))
      (inverse-dft! res))))

(defun main ()
  (let* ((n (read))
         (len (sb-int:power-of-two-ceiling (* +mod+ 2)))
         (counter (make-array +mod+ :element-type 'uint31 :initial-element 0))
         (powers (make-mod-power-table +generator+ (- +mod+ 1) +mod+))
         (where (make-array +mod+ :element-type 'uint31 :initial-element 0))
         (sub 0))
    (declare (uint62 sub n))
    (dotimes (i (length powers))
      (setf (aref where (aref powers i)) i))
    ;; #>where
    (dotimes (_ n)
      (let ((a (read-fixnum)))
        (declare (uint31 a))
        (unless (zerop a)
          (incf sub (mod (* a a) +mod+))
          (incf (aref counter (aref where a))))))
    (let ((counter (adjust-array (map '(simple-array double-float (*))
                                      (lambda (x) (float x 1d0))
                                      counter)
                                 len :initial-element 0d0)))
      (let ((f (fft-square! (copy-seq counter)))
            (res 0))
        (dotimes (i (length f))
          (let* ((coef (aref powers (mod i (- +mod+ 1))))
                 (count (round (aref f i))))
            (incf res (* coef count))))
        (println (ash (- res sub) -1))))))

#-swank (main)

;;;
;;; Test and benchmark
;;;

#+swank
(defun get-clipbrd ()
  (with-output-to-string (out)
    #+os-windows (run-program "powershell.exe" '("-Command" "Get-Clipboard") :output out :search t)
    #+os-unix (run-program "xsel" '("-b" "-o") :output out :search t)))

#+swank (defparameter *this-pathname* (uiop:current-lisp-file-pathname))
#+swank (defparameter *dat-pathname* (uiop:merge-pathnames* "test.dat" *this-pathname*))

#+swank
(defun run (&optional thing (out *standard-output*))
  "THING := null | string | symbol | pathname

null: run #'MAIN using the text on clipboard as input.
string: run #'MAIN using the string as input.
symbol: alias of FIVEAM:RUN!.
pathname: run #'MAIN using the text file as input."
  (let* ((*standard-output* (or out (make-string-output-stream)))
         (res (etypecase thing
                (null
                 (with-input-from-string (*standard-input* (delete #\Return (get-clipbrd)))
                   (main)))
                (string
                 (with-input-from-string (*standard-input* (delete #\Return thing))
                   (main)))
                (symbol (5am:run! thing))
                (pathname
                 (with-open-file (*standard-input* thing)
                   (main))))))
    (if out res (get-output-stream-string *standard-output*))))

#+swank
(defun gen-dat ()
  (uiop:with-output-file (out *dat-pathname* :if-exists :supersede)
    (format out "")))

#+swank
(defun bench (&optional (out (make-broadcast-stream)))
  (time (run *dat-pathname* out)))

;; To run: (5am:run! :sample)
#+swank
(it.bese.fiveam:test :sample
  (it.bese.fiveam:is
   (equal "474287
"
          (run "4
2019 0 2020 200002
" nil)))
  (it.bese.fiveam:is
   (equal "600013
"
          (run "5
1 1 2 2 100000
" nil))))


Measured memory_bytes: 50008