(in-package :cl-user)
(eval-when (:compile-toplevel :load-toplevel :execute)
  (defparameter *opt*
    #+swank '(optimize (speed 3) (safety 2))
    #-swank '(optimize (speed 3) (safety 0) (debug 0)))
  #+swank (ql:quickload '(:cl-debug-print :fiveam :cp/util) :silent t)
  #+swank (use-package :cp/util :cl-user)
  #-swank (set-dispatch-macro-character
           #\# #\> (lambda (s c p) (declare (ignore c p)) `(values ,(read s nil nil t))))
  #+sbcl (dolist (f '(:popcnt :sse4)) (pushnew f sb-c:*backend-subfeatures*))
  #+sbcl (setq *random-state* (seed-random-state (nth-value 1 (get-time-of-day)))))
#-swank (eval-when (:compile-toplevel)
          (setq *break-on-signals* '(and warning (not style-warning))))
#+swank (set-dispatch-macro-character #\# #\> #'cl-debug-print:debug-print-reader)

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

(defconstant +mod+ 1000000007)

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

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

;; BEGIN_INSERTED_CONTENTS
(defpackage :cp/two-phase-simplex
  (:use :cl)
  (:export #:dual-primal!)
  (:documentation "Provides two phase simplex method. (dual-primal)

Reference:
Robert J. Vanderbei. Linear Programming: Foundations and Extensions. 5th edition."))
(in-package :cp/two-phase-simplex)

;; NOTE: not optimized at all

(defconstant +eps+ 1d-8)

(defconstant +neg-inf+ most-negative-double-float)
(defconstant +pos-inf+ most-positive-double-float)

;; dict: col(0), col(1), ...., col(n-1), row(0), row(1), ..., row(m-1)
;;      |---------- non-basic ---------| |---------- basic ----------|

(declaim (ftype (function * (values double-float &optional)) %pivot))
(defun %pivot (row col a b c arow acol dict)
  (declare (optimize (speed 3))
           ((mod #.array-dimension-limit) row col)
           ((simple-array double-float (* *)) a)
           ((simple-array double-float (*)) b c arow acol)
           ((simple-array fixnum (*)) dict))
  (destructuring-bind (m n) (array-dimensions a)
    (declare ((mod #.array-dimension-limit) m n))
    (rotatef (aref dict col) (aref dict (+ n row)))
    (let* ((apivot (aref a row col))
           (/apivot (/ apivot)))
      (dotimes (i m)
        (dotimes (j n)
          (decf (aref a i j) (* (aref acol i) (aref arow j) /apivot))))
      (dotimes (j n)
        (setf (aref a row j) (* (aref arow j) /apivot)))
      (dotimes (i m)
        (setf (aref a i col) (- (* (aref acol i) /apivot))))
      (setf (aref a row col) /apivot)
      (let ((brow (aref b row)))
        (dotimes (i m)
          (decf (aref b i) (* brow (aref acol i) /apivot)))
        (setf (aref b row) (* brow /apivot))
        (let ((ccol (aref c col)))
          (dotimes (j n)
            (decf (aref c j) (* ccol (aref arow j) /apivot)))
          (setf (aref c col) (- (* ccol /apivot)))
          (* ccol brow /apivot))))))

(defun %restore (b c dict)
  (declare (optimize (speed 3))
           ((simple-array double-float (*)) b c)
           ((simple-array fixnum (*)) dict))
  (let* ((m (length b))
         (n (length c))
         (res-primal (make-array n :element-type 'double-float :initial-element 0d0))
         (res-dual (make-array m :element-type 'double-float :initial-element 0d0)))
    (dotimes (i m)
      (let ((index (aref dict (+ n i))))
        (when (< index n)
          (setf (aref res-primal index) (aref b i)))))
    (dotimes (j n)
      (let ((index (aref dict j)))
        (when (>= index n)
          (setf (aref res-dual (- index n)) (- (aref c j))))))
    (values res-primal res-dual)))

(declaim (ftype (function * (values (simple-array fixnum (*)) &optional))
                %iota))
(defun %iota (n)
  (let ((res (make-array n :element-type 'fixnum :initial-element 0)))
    (dotimes (i n)
      (setf (aref res i) i))
    res))

(defun %primal-simplex! (a b c &optional dict)
  "Assumes b >= 0."
  (declare (optimize (speed 3))
           ((simple-array double-float (* *)) a)
           ((simple-array double-float (*)) b c)
           ((or null (simple-array fixnum (*))) dict))
  (destructuring-bind (m n) (array-dimensions a)
    (declare ((mod #.array-dimension-limit) m n))
    (let ((acol (make-array m :element-type 'double-float))
          (arow (make-array n :element-type 'double-float))
          (dict (or dict (%iota (+ m n))))
          (obj 0d0))
      (declare (double-float obj))
      (loop
        ;; pick largest coefficient
        (let (col
              (colmax +neg-inf+))
          (dotimes (j n)
            (when (> (aref c j) colmax)
              (setq colmax (aref c j)
                    col j)))
          (unless col
            (error "All coefficient ~A are too small." c))
          (when (< colmax +eps+)
            (return))
          (dotimes (i m)
            (setf (aref acol i) (aref a i col)))
          ;; select leaving variable
          (unless (find-if (lambda (x) (> x +eps+)) acol)
            (return-from %primal-simplex! (values :unbounded nil nil nil)))
          (let (row
                (rowmin +pos-inf+))
            (dotimes (i m)
              (when (> (aref acol i) +eps+)
                (let ((rate (/ (aref b i) (aref acol i))))
                  (when (< rate rowmin)
                    (setq row i
                          rowmin rate)))))
            (unless row
              (error "Pivot not found in column ~A." acol))
            (dotimes (j n)
              (setf (aref arow j) (aref a row j)))
            ;; pivot
            (incf obj (%pivot row col a b c arow acol dict)))))
      ;; restore primal & dual solutions
      (multiple-value-bind (res-primal res-dual) (%restore b c dict)
        (values obj res-primal res-dual dict)))))

(defun %dual-simplex! (a b c &optional dict)
  "Assumes c <= 0."
  (declare (optimize (speed 3))
           ((simple-array double-float (* *)) a)
           ((simple-array double-float (*)) b c)
           ((or null (simple-array fixnum (*))) dict))
  (destructuring-bind (m n) (array-dimensions a)
    (declare ((mod #.array-dimension-limit) m n))
    (let ((acol (make-array m :element-type 'double-float))
          (arow (make-array n :element-type 'double-float))
          (dict (or dict (%iota (+ m n))))
          (obj 0d0))
      (declare (double-float obj))
      (loop
        ;; pick least intercept
        (let (row
              (rowmin +pos-inf+))
          (dotimes (i m)
            (when (< (aref b i) rowmin)
              (setq rowmin (aref b i)
                    row i)))
          (unless row
            (error "All intercepts ~A are too large." b))
          (when (> rowmin +eps+)
            (return))
          (dotimes (j n)
            (setf (aref arow j) (aref a row j)))
          ;; select leaving variable
          (unless (find-if (lambda (x) (< x (- +eps+))) arow)
            (return-from %dual-simplex! (values :infeasible nil nil nil)))
          (let (col
                (colmin +pos-inf+))
            (dotimes (j n)
              (when (< (aref arow j) (- +eps+))
                (let ((rate (/ (aref c j) (aref arow j))))
                  (when (< rate colmin)
                    (setq col j
                          colmin rate)))))
            (unless col
              (error "Pivot not found in row ~A." arow))
            (dotimes (i m)
              (setf (aref acol i) (aref a i col)))
            (incf obj (%pivot row col a b c arow acol dict)))))
      ;; restore primal & dual solutions
      (multiple-value-bind (res-primal res-dual) (%restore b c dict)
        (values obj res-primal res-dual dict)))))

(declaim (ftype (function * (values (or double-float (member :unbounded :infeasible))
                                    (or null (simple-array double-float (*)))
                                    (or null (simple-array double-float (*)))
                                    (or null (simple-array fixnum (*)))
                                    &optional))
                dual-primal!))
(defun dual-primal! (a b c)
  "Maximizes cx subject to Ax <= b and x >= 0. Returns four values:

Optimal case:
- optimal objective value
- solutions for primal problem
- solutions for dual problem: min. by s.t. (A^t)y >= c and y >= 0
- current dictionary

Unbounded case:
- (values :unbounded nil nil nil)

Infeasible case:
- (values :infeasible nil nil nil)"
  (declare (optimize (speed 3))
           ((simple-array double-float (* *)) a)
           ((simple-array double-float (*)) b c))
  (destructuring-bind (m n) (array-dimensions a)
    (declare ((mod #.array-dimension-limit) m n))
    ;; Phase I: solve modified problem with objective = -x1-x2- ... -xn
    (let* ((c* (make-array n :element-type 'double-float :initial-element -1d0))
           (dict (%iota (+ n m)))
           (result1 (%dual-simplex! a b c* dict)))
      (when (eql result1 :infeasible)
        (return-from dual-primal! (values :infeasible nil nil nil)))
      ;; restore original objective function
      (let ((poses (make-array (+ n m) :element-type 'fixnum)))
        (dotimes (i (+ m n))
          (setf (aref poses (aref dict i)) i))
        (replace c* c)
        (fill c 0d0)
        (let ((obj 0d0))
          (declare (double-float obj))
          (dotimes (j n)
            (let* ((coef (aref c* j))
                   (pos (aref poses j)))
              (if (< pos n)
                  ;; xj is non-basic
                  (let ((col pos))
                    (incf (aref c col) coef))
                  ;; xj is basic
                  (let ((row (- pos n)))
                    (dotimes (j n)
                      (decf (aref c j) (* coef (aref a row j))))
                    (incf obj (* coef (aref b row)))))))
          ;; Phase II: solve original problem
          (multiple-value-bind (result2 primal dual) (%primal-simplex! a b c dict)
            (if (eql result2 :unbounded)
                (values result2 nil nil nil)
                (values (+ obj (the double-float result2)) primal dual dict))))))))

;; BEGIN_USE_PACKAGE
(eval-when (:compile-toplevel :load-toplevel :execute)
  (use-package :cp/two-phase-simplex :cl-user))
(in-package :cl-user)

;;;
;;; Body
;;;

(defun main ()
  (let* ((c (float (read) 1d0))
         (d (float (read) 1d0))
         (obj (make-array 2 :element-type 'double-float :initial-contents '(1000d0 2000d0)))
         (bs (coerce (vector c d) '(simple-array double-float (*))))
         (as (make-array '(2 2)
                         :element-type 'double-float
                         :initial-contents '((0.75d0 #.(float 2/7 1d0))
                                             (0.25d0 #.(float 5/7 1d0))))))
    (let ((res (dual-primal! as bs obj)))
      (println res))))

#-swank (main)

;;;
;;; Test
;;;

#+swank
(progn
  (defparameter *lisp-file-pathname* (uiop:current-lisp-file-pathname))
  (setq *default-pathname-defaults* (uiop:pathname-directory-pathname *lisp-file-pathname*))
  (uiop:chdir *default-pathname-defaults*)
  (defparameter *dat-pathname* (uiop:merge-pathnames* "test.dat" *lisp-file-pathname*))
  (defparameter *problem-url* "https://yukicoder.me/problems/771"))

#+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)))

#+(and sbcl (not swank))
(eval-when (:compile-toplevel)
  (when sb-c::*undefined-warnings*
    (error "undefined warnings: ~{~A~^ ~}" sb-c::*undefined-warnings*)))

;; To run: (5am:run! :sample)
#+swank
(5am:test :sample
  (5am:is
   (equal "18000
"
          (run "5 6
" nil)))
  (5am:is
   (equal "19692.3076923077
"
          (run "3 7
" nil)))
  (5am:is
   (equal "40000
"
          (run "32 10
" nil)))
  (5am:is
   (equal "3076923.0769230775
"
          (run "1000 1000
" nil))))