ZOJ 3614 Choir 二维前缀和 二维rmq
http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3614
题意:
给定n*m的矩阵,多次询问,对于给定a和b,找一个a*b的子矩阵,使得去掉该矩阵的一个最大值后,剩余数方差最小。
分析:
方差可以展开,通过前缀和预处理,之后枚举子矩阵左上角就可以o(1)计算方差。而去掉最大值就是rmq的工作了。代码很挫,好久没写二维rmq,开始写跪了,然后前缀和写得也不好,所以偷偷地再贴一份别人的。。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<iostream> 3 #include<cstring> 4 #include<cmath> 5 #include<algorithm> 6 #include<cstdlib> 7 #include<set> 8 #include<map> 9 #include<queue> 10 #include<ctime> 11 #include<string> 12 using namespace std; 13 14 const double eps = 1e-6; 15 int a, b, n, m, q, cas; 16 long long sum[400][400], sum2[400][400], v[400][400], v2[400][400]; 17 int dp[310][310][10][10]; 18 long long max(long long a, long long b) 19 { 20 if (a > b) return a; 21 else return b; 22 } 23 void init() 24 { 25 memset(sum, 0, sizeof(sum)); 26 memset(sum2, 0, sizeof(sum2)); 27 memset(dp, 0, sizeof(dp)); 28 for (int i = 1; i <= n; i++){ 29 sum[i][0] = 0; 30 sum2[i][0] = 0; 31 for (int j = 1; j <= m; j++){ 32 sum[i][j] = sum[i][j-1] + v[i][j]; 33 sum2[i][j] = sum2[i][j-1] + v2[i][j]; 34 } 35 } 36 for (int j = 1; j <= m; j++) 37 for (int i = 1; i <= n; i++){ 38 sum[i][j] = sum[i-1][j] + sum[i][j]; 39 sum2[i][j] = sum2[i-1][j] + sum2[i][j]; 40 } 41 for (int i = 1; i <= n; i++) 42 for (int j = 1; j <= m; j++) 43 dp[i][j][0][0] = v[i][j]; 44 for (int k = 0; (1 << k) <= n; k++){ 45 for (int l = 0; (1 << l) <= m; l++){ 46 if (k == 0 && l == 0) continue; 47 int t1 = 1 << (k-1), 48 t2 = 1 << (l-1); 49 t1 = max(0, t1), t2 = max(0, t2); 50 for (int i = 1; i + t1 <= n; i++) 51 for (int j = 1; j + t2 <= m; j++){ 52 if (k == 0 && l != 0){ 53 dp[i][j][k][l] = max(dp[i][j][k][l-1], dp[i][j+t2][k][l-1]); 54 } 55 else if (k != 0){ 56 dp[i][j][k][l] = max(dp[i][j][k-1][l], dp[i+t1][j][k-1][l]); 57 } 58 //printf("dp[%d][%d][%d][%d] %d\n", i, j, k, l, dp[i][j][k][l]); 59 } 60 } 61 } 62 } 63 int rmq(int x, int y, int xx, int yy) 64 { 65 int l = int(log(double(xx - x + 1))/log(2.0)), ll = int(log(double(yy - y + 1))/log(2.0)); 66 int maxx = dp[x][y][l][ll]; 67 maxx = max(maxx, dp[xx+1-(1<<l)][yy+1-(1<<ll)][l][ll]); 68 maxx = max(maxx, dp[xx+1-(1<<l)][y][l][ll]); 69 maxx = max(maxx, dp[x][yy+1-(1<<ll)][l][ll]); 70 return maxx; 71 } 72 void solve(int a, int b) 73 { 74 double ans = 1e20; 75 int ansx, ansy; 76 for (int i = 1; i+a-1 <= n; i++) 77 for (int j = 1; j+b-1 <= m; j++){ 78 long long high = rmq(i, j, i+a-1, j+b-1); 79 long long sumx = sum[i+a-1][j+b-1] - sum[i+a-1][j-1] - sum[i-1][j+b-1] + sum[i-1][j-1] - high; 80 //printf("%lld %lld %lld %lld\n", sum[i+a-1][j+b-1], sum[i+a-1][j-1], sum[i-1][j+b-1], sum[i-1][j-1]); 81 long long sumx2 = sum2[i+a-1][j+b-1] - sum2[i+a-1][j-1] - sum2[i-1][j+b-1] + sum2[i-1][j-1] - high * high; 82 83 //printf("%d %d %d %d %lld %lld %lld\n", i, j, i+a-1, j+b-1, high, sumx, sumx2); 84 85 double ave = (double)sumx / (double)(a*b-1); 86 double tmp = (double)sumx2 - ave * sumx * 2 + ave * ave * (a * b - 1); 87 //printf("%.4lf\n", tmp); 88 if (ans - tmp > eps){ 89 ans = tmp; 90 ansx = i; 91 ansy = j; 92 } 93 } 94 printf("(%d, %d), %.2lf\n", ansx, ansy, ans/(double)(a*b-1)); 95 } 96 int main() 97 { 98 freopen("1.in", "r", stdin); 99 freopen("2.out", "w", stdout); 100 cas = 0; 101 while(scanf("%d %d", &n, &m) != EOF) 102 { 103 for (int i = 1; i <= n; i++) 104 for (int j = 1; j <= m; j++){ 105 scanf("%lld", &v[i][j]); 106 v2[i][j] = v[i][j] * v[i][j]; 107 } 108 init(); 109 printf("Case %d:\n", ++cas); 110 scanf("%d", &q); 111 for (int i = 0; i < q; i++){ 112 scanf("%d %d", &a, &b); 113 solve(a, b); 114 } 115 } 116 return 0; 117 }
别人的:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 /* 2 *二维RMQ ST算法 3 *构造RMQ数组 makermq(int n,int m,int b[][]) O(n*m*log(n)*log(m))算法复杂度 4 *st[i][j][row][col] 表示 行从row ->row +2^i-1 列从col ->col +2^j-1 二维区间里最大值 5 *st[row][col][i][j] = 下行 6 *max{st[i][j-1][row][col],st[i-1][j][row][col],st[i][j-1][row][col+2^(j-1)],st[col][i-1][j][row+2^(i-1)]} 7 *查询RMQ rmq(int sx,int ex,int sy,int ey) 8 *同一维的将sx->ex 分为两个2^kx区间 将 sy->ey分为两个2^ky的区间 9 *kx=(int)log2(ex-sx+1) ky=(int)log2(ey-sy+1) 10 *log2()可能会CE 11 *查询结果为 12 *max{st[kx][ky][sx][sy],st[kx][ky][sx][ey-2^ky+1], 13 *st[kx][ky][ex-2^kx+1][sy],st[kx][ky][ex-2^kx+1][ey-2^ky+1]} 14 */ 15 #include <stdio.h> 16 #include <string.h> 17 #include <iostream> 18 #include <math.h> 19 using namespace std; 20 typedef long long ll; 21 const int N = 310; 22 const int M = 10; 23 int map[N][N]; 24 int st[M][M][N][N]; 25 ll A[N][N],B[N][N]; 26 int POW[] = { 1, 1 << 1, 1 << 2, 1 << 3, 1 << 4, 1 << 5, 1 << 6, 1 << 7, 27 1 << 8, 1 << 9, 1 << 10 }; 28 void makermq(int n, int m, int b[][N]) { 29 int row, col, i, j; 30 for (row = 1; row <= n; row++) { 31 for (col = 1; col <= m; col++) { 32 st[0][0][row][col] = b[row][col]; 33 } 34 } 35 for (i = 0; POW[i] <= n; i++) { 36 for (j = 0; POW[j] <= m; j++) { 37 if (i == 0 && j == 0) 38 continue; 39 for (row = 1; row + POW[i] - 1 <= n; row++) { 40 for (col = 1; col + POW[j] - 1 <= m; col++) { 41 if (i == 0) { 42 st[i][j][row][col] = max(st[i][j - 1][row][col], 43 st[i][j - 1][row][col + POW[j - 1]]); 44 } else { 45 st[i][j][row][col] = max(st[i - 1][j][row][col], st[i 46 - 1][j][row + POW[i - 1]][col]); 47 } 48 } 49 } 50 } 51 } 52 } 53 54 inline int get_k(int a, int b) { 55 return (int) (log((b - a + 1) * 1.0)/log(2.0)); 56 } 57 58 inline int RMQ(int x1, int x2, int y1, int y2) { 59 60 int kx = get_k(x1, x2), ky = get_k(y1, y2); 61 int max1, max2;// min1, min2; 62 max1 = max(st[kx][ky][x1][y1], st[kx][ky][x1][y2 - POW[ky] + 1]); 63 max2 = max(st[kx][ky][x2 - POW[kx] + 1][y1], 64 st[kx][ky][x2 - POW[kx] + 1][y2 - POW[ky] + 1]); 65 return max(max1, max2) ;//- min(min1, min2); 66 } 67 68 double cal(int x1, int y1, int x2, int y2 ,int n) { 69 ll re = (long long)RMQ(x1, x2, y1, y2); 70 ll aa = (long long)A[x2][y2]+A[x1-1][y1-1]-A[x2][y1-1]-A[x1-1][y2]; 71 ll bb = (long long)B[x2][y2]+B[x1-1][y1-1]-B[x2][y1-1]-B[x1-1][y2]; 72 aa-=re, bb-=re*re; 73 double ans = 1.0*(bb-1.0*aa*aa/n)/n; 74 return ans; 75 } 76 77 int main(){ 78 freopen("1.in", "r", stdin); 79 freopen("1.out", "w", stdout); 80 int q,i,j, ca = 1; 81 int n, m; 82 while(~scanf("%d%d",&n,&m)){ 83 for (i=1;i<=n;i++) { 84 for (j = 1;j <= m;j++) scanf("%d", &map[i][j]); 85 } 86 for (i = 0;i <= n;i++) A[i][0] = B[i][0] = 0; 87 for (i = 0;i <= m;i++) B[0][i] = A[0][i] = 0; 88 for (i = 1;i <= n;i++) { 89 for (j = 1;j <= m;j++) { 90 A[i][j] = A[i-1][j]+A[i][j-1]-A[i-1][j-1]+1ll*map[i][j]; 91 B[i][j] = B[i-1][j]+B[i][j-1]-B[i-1][j-1]+1ll*map[i][j]*map[i][j]; 92 } 93 } 94 makermq(n, m, map); 95 scanf("%d", &q); 96 int tx, ty; 97 printf("Case %d:\n", ca++); 98 while (q--) { 99 scanf("%d%d", &tx, &ty); 100 double ans = -1; 101 int ax, ay; 102 for (i = 1;i <= n-tx+1;i++) { 103 for (j = 1;j <= m-ty+1;j++) { 104 double tm = cal(i, j, tx+i-1, ty+j-1,(tx*ty-1)); 105 if (ans == -1 || tm < ans) { 106 ans = tm, ax = i, ay = j; 107 } 108 } 109 } 110 printf("(%d, %d), %.2lf\n", ax, ay, ans); 111 } 112 } 113 }