ICPC 2014 Taichung F

少し前に練習会で2014の台湾の問題をやった時に，まったく解けなかった．
しかし，CodechefのLTIME29に参加した後に解法が思いついて実装したら通ってしまった．
解いているチーム数も少ないので，解説を書きます．

問題概要

n行m列の行列が与えられる．この行列の中から要素数k以上となるように部分行列を選択する．
選択できる部分行列の中で，部分行列の要素の平均値の最大値を四捨五入して，小数点以下第三位まで求めよ．

制約



行列の各要素は5000以下の非負整数　

解法
まず，二次元累積和を計算する．
次に，選択する部分行列の行の範囲を固定して考える．
そうすると，列方向における1次元の累積和sumが計算できる．
ここで，部分行列の行数をrowと定義する．また，行数が固定された時にkより大きくなるのに必要な列数の下限はlen=ceil(k/row)と決まる．

そして，平均最大化を考える
二分探索を使って，平均値X以上になるように列方向の要素を選択できるかを考える．
この時に，　　選択した部分行列の要素の総和/選択した部分行列の要素数>=x　という判定式を考える．
この判定式を変形すると，1個の要素を選択した時にかかるコストは，行列の要素-row*Xであり，このコストの総和が0以上になれば選択できると考えられる．
そこで，sum[i]に -row*i*x のペナルティを追加する．
２つの累積和を選択する間隔がlen以上になるように前からiについてのループを回して次の式の最大値を求める．sum[i+len]-min(sum[0],sum[1],...,sum[i])
この時に，最大値0以上になれば選択できると判定できる．

あとは，二分探索を下限を更新しながら行の範囲を全探索すれば良い．

1テストケースあたりの計算量は，Tを二分探索の回数とすると， 求める平均値の最大値は制約より5000で，出力の精度から二分探索の区間が0.0005以下になれば良いので，T=25くらい．
上限値をあてはめると，600*(1000000)*25=1.5*10^10となる．通らなさそう...
しかし，考察をしてみると，調べるべき部分行列の大きさSとした時に，k<=S<=2kの範囲を調べれば良いことが分かる．(参考 JAG 2013 Day 4 C FoxObservation)
これにより探索範囲は大幅に枝刈りできるので，ICPC Live Archiveで実行時間10sくらいで通る．(ICPC Live ArchiveのTLEが60sなので，枝刈りしなくても通るかも)

しかし，大会本番のTLEは3sということでさらなる高速化を施す．
毎回の二分探索を開始する前に，現在の下限値以上になるように選択できるかを調べておく．そうすると，選択できない場合には二分探索をしなく済む．
ここまで，やると1sくらい通ります．

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

#define range(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,n)  range(i,0,n)

const double eps = 1e-10;

ll mat[610][10010];

int n,m,k;
double tmp[10010];

bool check(double x,int row,int len){
	double cmin=0.0;
	rep(i,m-len+1){
		cmin=min(cmin,tmp[i]-1.0*i*row*x);
		double cur=tmp[i+len]-1.0*(i+len)*row*x-cmin;
		if(cur>eps) return true;
	}
	return false;
}

int main(void){
	int t;
	scanf("%d",&t);
	rep(loop,t){
		scanf("%d %d %d",&n,&m,&k);
		rep(i,n)rep(j,m) scanf("%lld",&(mat[i+1][j+1]));
		rep(i,n)rep(j,m) mat[i+1][j+1]+=mat[i+1][j]+mat[i][j+1]-mat[i][j];

		double left=0.0,right=5001.0;
		rep(ui,n+1)rep(li,ui){
			int row=ui-li;
			if(2*k<row) continue;
			if(row*m<k) continue;
			right=5001.0;
			int len=(k+row-1)/row;
			rep(j,m+1) tmp[j]=1.0*mat[ui][j]-mat[li][j];

			if(check(left,row,len)==false) continue; // important points
			while(right-left>0.0005){
				double mid=(left+right)/2.0;
				if(check(mid,row,len))
					left=mid;
				else
					right=mid;
			}
		}
		printf("%.3f\n",left);
	}
	return 0;
}