結果
| 問題 |
No.453 製薬会社
|
| コンテスト | |
| ユーザー |
 sansaqua
|
| 提出日時 | 2021-05-19 23:33:01 |
| 言語 | Common Lisp (sbcl 2.5.0) |
| 結果 |
AC
|
| 実行時間 | 12 ms / 2,000 ms |
| コード長 | 12,459 bytes |
| コンパイル時間 | 1,568 ms |
| コンパイル使用メモリ | 56,756 KB |
| 実行使用メモリ | 31,548 KB |
| 最終ジャッジ日時 | 2024-10-10 05:39:10 |
| 合計ジャッジ時間 | 2,404 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 9 |
コンパイルメッセージ
; compiling file "/home/judge/data/code/Main.lisp" (written 10 OCT 2024 05:39:07 AM): ; file: /home/judge/data/code/Main.lisp ; in: DEFUN %PIVOT ; (DEFUN CP/TWO-PHASE-SIMPLEX::%PIVOT ; (CP/TWO-PHASE-SIMPLEX::ROW CP/TWO-PHASE-SIMPLEX::COL ; CP/TWO-PHASE-SIMPLEX::A CP/TWO-PHASE-SIMPLEX::B ; CP/TWO-PHASE-SIMPLEX::C CP/TWO-PHASE-SIMPLEX::AROW ; CP/TWO-PHASE-SIMPLEX::ACOL CP/TWO-PHASE-SIMPLEX::DICT) ; (DECLARE (OPTIMIZE (SPEED 3)) ; ((MOD 17592186044416) CP/TWO-PHASE-SIMPLEX::ROW ; CP/TWO-PHASE-SIMPLEX::COL) ; ((SIMPLE-ARRAY DOUBLE-FLOAT (* *)) CP/TWO-PHASE-SIMPLEX::A) ; ((SIMPLE-ARRAY DOUBLE-FLOAT (*)) CP/TWO-PHASE-SIMPLEX::B ; CP/TWO-PHASE-SIMPLEX::C CP/TWO-PHASE-SIMPLEX::AROW ; CP/TWO-PHASE-SIMPLEX::ACOL) ; ((SIMPLE-ARRAY FIXNUM (*)) CP/TWO-PHASE-SIMPLEX::DICT)) ; (DESTRUCTURING-BIND ; (CP/TWO-PHASE-SIMPLEX::M CP/TWO-PHASE-SIMPLEX::N) ; (ARRAY-DIMENSIONS CP/TWO-PHASE-SIMPLEX::A) ; (DECLARE ; ((MOD 17592186044416) CP/TWO-PHASE-SIMPLEX::M CP/TWO-PHASE-SIMPLEX::N)) ; (ROTATEF (AREF CP/TWO-PHASE-SIMPLEX::DICT CP/TWO-PHASE-SIMPLEX::COL) ; (AREF CP/TWO-PHASE-SIMPLEX::DICT ; (+ CP/TWO-PHASE-SIMPLEX::N CP/TWO-PHASE-SIMPLEX::ROW))) ; (LET* ((CP/TWO-PHASE-SIMPLEX::APIVOT #) (CP/TWO-PHASE-SIMPLEX::/APIVOT #)) ; (DOTIMES (CP/TWO-PHASE-SIMPLEX::I CP/TWO-PHASE-SIMPLEX::M) ; (DOTIMES # #)) ; (DOTIMES (CP/TWO-PHASE-SIMPLEX::J CP/TWO-PHASE-SIMPLEX::N) (SETF #)) ; (DOTIMES (CP/TWO-PHASE-SIMPLEX::I CP/TWO-PHASE-SIMPLEX::M) (SETF #)) ; (SETF # CP/TWO-PHASE-SIMPLEX::/APIVOT) ; (LET (#) ; (DOTIMES # #) ; (SETF #) ; (LET # ; # ; # ; #))))) ; --> SB-IMPL::%DEFUN SB-IMPL::%DEFUN SB-INT:NAMED-LAM
ソースコード
(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))))