sdoi最优图像image（网络流）

题意：有一个N*M的01矩阵a，已知每一行有多少个1以及每一列有多少1。现在这个矩阵搞丢了，但是会告诉你a[i,j]是1的概率p[i,j]（一个[0,100)的整数表示百分率）。让你还原出一个概率最大的符合条件的01矩阵，任意输出一个。

据说这题只能用zkw费用流过。。这个太奇怪了。

这题调了我一个下午，很有启发意义。

首先这是经典的矩阵还原模型，就是行做X部，列做Y部，之间的连边的流量代表a[i,j]的值，之前做过两道类似的题，因此这里还是想到了的。

然后就是最精华的地方了，概率应该是所有的p[i,j]的乘积，但是费用流的模型是和，我开始试着改一下费用流的写法，就是把中间找最短路的地方改成乘法，这样理论上可行但是乘不到几个数就会严重暴long long了。。然后我就上网搜了，我惊讶地发现竟然没几个人写这道题的题解，并且没人贴程序。这题有这么偏门吗。。找了半天才找到，要把费用设为概率的对数！这个真是太巧妙了，因为两个正数的大小关系完全等同于他们取对数的大小，然后两个数相乘就对应了他们的对数相加！！！然后要注意因为是最大费用，将他们费用全部取反，最后再取反回来即可。

然后我就无脑地写了一发，然后就交了，结果只过了前四个点，后面全T了。。然后我就各种调，我开始以为是连边有问题产生了负权环，然后下了组数据写了个SPFA判定，结果没负环。然后我就翻来覆去看了几遍费用流的模板，觉得没问题，，然后就抓狂了，在各种细节上试，，最后终于试出来了，就是update里面判是否是最短路径上的点那个dis的判定我直接用的等号，但事实上由于精度问题转移几次之后该相等的已经不相等了，于是我写了个eps和自定义了个等号。。好消息，没T了，但是都WA了，，，然后我又找了半天，觉得好像取对数的时候精度不好，于是把概率全部乘了很大一个数之后再取的对数，然后终于过了。以前openjudge上面有道题就是，要先乘一个数再做除法，不然精度要爆。

这题没看到别人写的代码，也完全不知道精度会有这种问题，全是自己摸索的，写出来还是很愉快。

这道题的启发：

1、比较乘积但是太大了，可以换成取对数之后的加法。

2、像什么取对数，开根号之类容易挂精度的，可以先乘一个很大的数。

3、系统的对数函数好像并不是O(1)的（或者是常数有点大）。系统自带的有log，log2，log10，分别是以e,2,10为底，其中以2位底的会比另外两个快三倍，如果只是用来比大小，最后用log2。

#include<cstdio>#include<cstring>#include<cassert>#include<cmath>#include<algorithm>using namespace std;#define DB double#define clr(a) memset(a,0,sizeof a)inline DB min(const DB&a, const DB&b){return a < b ? a : b; }#define rep(a,b,c) for (int a=b;a<=c;++a)const int inf = 0x3f3f3f3f;const int MAXN = 205;const int MAXM = 100000;int N, M;const DB eps = 1e-10;bool cmp(DB a, DB b) //自己定义在精度误差范围内的等号{if (a < b) swap(a, b);return a-b <= eps;}struct Ed {int to, cap;DB cost;Ed*nxt, *back;};struct FlowNet{Ed Edge[MAXM], *ecnt, *adj[MAXN];FlowNet () { ecnt=Edge; }DB dis[MAXN];bool vis[MAXN];int vn, S, T, flow;DB tot;inline void adde(int a, int b, int c, DB d){(++ecnt)->to = b;ecnt->cap = c;ecnt->cost = d;ecnt->nxt = adj[a];ecnt->back = ecnt+1;adj[a] = ecnt;(++ecnt)->to = a;ecnt->cap = 0;ecnt->cost = -d;ecnt->nxt = adj[b];ecnt->back = ecnt-1;adj[b] = ecnt;}void init(int n, int s, int t){tot = 0; flow = 0;vn = n; S = s; T = t;clr(dis), clr(adj);ecnt = Edge;}bool update(){DB tmp = 1e10;rep(i, 1, vn) if (vis[i])for (Ed *p = adj[i]; p; p=p->nxt)if (p->cap > 0 && !vis[p->to])tmp = min(tmp, dis[p->to]-dis[i] + p->cost);if (tmp == 1e10) return 0;for (int i = 1; i<=vn; ++i)if (vis[i]) dis[i] += tmp;return 1;}int aug(int u, int augco){if (u == T){tot += dis[S] * augco;flow += augco;return augco;}vis[u] = 1;int delta, augc = augco;for (Ed*p = adj[u]; p && augc; p=p->nxt){int&v = p->to;if (!vis[v] && p->cap && cmp(dis[u],dis[v]+p->cost)) //写等号要出问题{delta = min(p->cap, augc);delta = aug(v, delta);p->cap -= delta, p->back->cap += delta;augc -= delta;}}return augco - augc;}DB mcmf(){do {do clr(vis); while (aug(S, inf));} while (update());return tot;}void drawMap(int img[105][105]){rep(i, 1, N)for (Ed*p = adj[i]; p; p=p->nxt)if (p->to > i && p->to<=N+M)img[i][p->to - N] = p->back->cap;}} G;int img[105][105];int main(){scanf("%d%d", &N, &M);G.init(N+M+2, N+M+1, N+M+2);DB tmp;rep(i, 1, N) rep(j, 1, M){scanf("%lf", &tmp);G.adde(i, j+N, 1, -log2(tmp*1e6)); //调整精度}int xx;rep(i, 1, N){scanf("%d", &xx);G.adde(G.S, i, xx, 0);}rep(j, 1, M){scanf("%d", &xx);G.adde(N+j, G.T, xx, 0);}G.mcmf();G.drawMap(img);rep(i, 1, N){rep(j, 1, M) printf("%d", img[i][j]);puts("");}return 0;}