BZOJ 4514: [Sdoi2016]数字配对 [费用流 数论]

4514: [Sdoi2016]数字配对

题意：
有 n 种数字，第 i 种数字是 ai、有 bi 个，权值是 ci。
若两个数字 ai、aj 满足，ai 是 aj 的倍数，且 ai/aj 是一个质数，
那么这两个数字可以配对，并获得 ci×cj 的价值。
一个数字只能参与一次配对，可以不参与配对。
在获得的价值总和不小于 0 的前提下，求最多进行多少次配对。

---

显然可以配对的两点之间可以连费用为\(c_i \times c_j\)的边
一开始想拆开节点限制流量，但这样没法求配对次数啊
应该深入分析连边两点的性质
因为差一个质因子，只有可能是奇数个质因子向偶数个质因子连边
这样就变成二分图了，在与s和t的连边上限制流量就行了

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long ll;
#define fir first
#define sec second
const int N=1005, E=1e5+5, M=32000, INF=1e9;
inline int read() {
	char c=getchar(); int x=0, f=1;
	while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
	while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
	return x*f;
}

int n, a[N], b[N], s, t; ll c[N];
struct edge{int v, c, f, ne; ll w;}e[M];
int cnt=1, h[N];
inline void ins(int u, int v, int c, ll w) { //printf("ins %d %d %d  %lld\n",u,v,c,w);
	e[++cnt]=(edge){v, c, 0, h[u], w}; h[u]=cnt;
	e[++cnt]=(edge){u, 0, 0, h[v], -w}; h[v]=cnt;
}

int q[N], head, tail, inq[N]; ll d[N];
pair<int, int> pre[N];
inline void lop(int &x) {if(x==N) x=1;else if(x==0) x=N-1;}
bool spfa(int s, int t) {
	for(int i=s; i<=t; i++) inq[i]=0, d[i]=-1e15;
	head=tail=1;
	q[tail++]=s; inq[s]=1; d[s]=0;
	pre[t].fir = -1;
	while(head!=tail) {
		int u=q[head++]; inq[u]=0; lop(head);
		for(int i=h[u];i;i=e[i].ne) {
			int v=e[i].v;
			if(d[v]<d[u]+e[i].w && e[i].c>e[i].f) {
				d[v]=d[u]+e[i].w; 
				pre[v] = make_pair(u, i);
				if(!inq[v]) {
					inq[v]=1;
					if(d[v]>d[q[head]]) head--, lop(head), q[head]=v;
					else q[tail++]=v, lop(tail);
				}
			}
		}
	}
	return pre[t].fir != -1;
}
int ek(int s, int t) {
	int flow=0; ll cost=0;
	while(spfa(s, t)) {
		int f=INF, x;
		for(int i=t; i!=s; i=pre[i].fir) x=pre[i].sec, f=min(f, e[x].c-e[x].f);
		//printf("hi %d %d  %d\n",f,d[t],cost);
		if(d[t]<0 && d[t]*f+cost<0) {
			int x = cost/(-d[t]); //printf("x %d\n",x);
			flow+=x; 
			break;
		} else flow+=f, cost+=d[t]*f;
		for(int i=t; i!=s; i=pre[i].fir) x=pre[i].sec, e[x].f+=f, e[x^1].f-=f;
	}
	return flow;
}

int p[M], notp[M];
void sieve(int n) {
	for(int i=2; i<=n; i++) {
		if(!notp[i]) p[++p[0]]=i;
		for(int j=1; j<=p[0] && i*p[j]<=n; j++) {
			notp[i*p[j]]=1;
			if(i%p[j]==0) break;
		}
	}
}
int fact(int x) {
	int m=sqrt(x), ans=0;
	for(int i=1; p[i]<=m; i++) 
		while(x%p[i]==0) ans++, x/=p[i];
	if(x!=1) ans++;
	return ans;
}
inline bool isp(int x) {
	if(x==1) return false;
	if(x==2) return true;
	int m=sqrt(x);
	for(int i=1; p[i]<=m; i++) if(x%p[i]==0) return false;
	return true;
}
inline bool legal(int a, int b) {
	if(a<b) swap(a, b);
	if(a%b) return false;
	return isp(a/b);
}
int le[N], ri[N];
void build() {
	s=0; t=n+n+1;
	for(int i=1; i<=n; i++) {
		if(fact(a[i])&1) le[++le[0]]=i;
		else ri[++ri[0]]=i;
	}
	//for(int i=1; i<=le[0]; i++) printf("%d ",le[i]); puts(" le");
	//for(int i=1; i<=ri[0]; i++) printf("%d ",ri[i]); puts(" ri");

	for(int i=1; i<=le[0]; i++)
		for(int j=1; j<=ri[0]; j++) 
			if(legal(a[le[i]], a[ri[j]])) ins(le[i], n+ri[j], INF, c[le[i]]*c[ri[j]]);// printf("legal %d %d\n",le[i], ri[j]);;
	for(int i=1; i<=le[0]; i++) ins(s, le[i], b[le[i]], 0);
	for(int i=1; i<=ri[0]; i++) ins(n+ri[i], t, b[ri[i]], 0);
}
int main() {
	freopen("in","r",stdin);
	sieve(M-1); 
	n=read();
	for(int i=1; i<=n; i++) a[i]=read();
	for(int i=1; i<=n; i++) b[i]=read();
	for(int i=1; i<=n; i++) c[i]=read();
	build();
	int flow=ek(s, t);
	printf("%d", flow);
}