BZOJ.4514.[SDOI2016]数字配对(费用流SPFA 二分图)

时间：2023-03-08 18:21:37

BZOJ

洛谷

\(Solution\)

很显然的建二分图后跑最大费用流，但有个问题是一个数是只能用一次的，这样二分图两部分都有这个数。

那么就用两倍的。如果\(i\)可以向\(j'\)连边，\(j\)也向\(i'\)连边，如果上一次走了\(i->j'\)，那么这一次一定走\(j->i'\)。

每次跑最大费用流，直至有一次费用变成负，然后加上当前正权值能抵消它的流量，最后总流量除以2就可以了。

\(Another Solution\)

两个数能匹配首先要能整除，其次它们所有质因子的次数和一定只相差1.

于是可以按这个次数的奇偶性建二分图，分别连源点汇点，一遍最大流就可以了。

ps：sb的我硬是把\(i\to j\)的费用设成\(0\)，然后源点汇点连\(\log b[i]\)的边，然后再取幂。。最后写完发现不对，不会处理上面的问题。。贪心也调不对。

数字匹配\(O(sqrt(n))\)直接枚举就可以，闲的写什么MillerRabin+set。。

考试的时候一个题写了三百多行 还不算颓

时隔一年后的考试写出来了2333（虽然写过还是道水题）

第一种方法：

//3212kb	112ms

#include <set>

#include <queue>

#include <cmath>

#include <cstdio>

#include <cctype>

#include <algorithm>

#define gc() getchar()

typedef long long LL;

const int N=410,M=100005,INF=1e9;

const int P[9]={2,3,5,7,13,17,31};

int n,src,des,A[N],B[N],C[N],Enum,H[N],nxt[M],fr[M],to[M],cap[M],pre[N];

bool inq[N];

LL dis[N],cost[M];

std::queue<int> q;

std::set<int> is_p[2];

inline int read()

{

	int now=0,f=1;register char c=gc();

	for(;!isdigit(c);c=gc()) if(c=='-') f=-1;

	for(;isdigit(c);now=now*10+c-'0',c=gc());

	return now*f;

}

inline void AddEdge(int u,int v,int w,LL c)

{

	to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, cap[Enum]=w, cost[Enum]=c;

	to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, cap[Enum]=0, cost[Enum]=-c;

}

int FP(LL x,int k,int p)

{

	LL t=1;

	for(; k; k>>=1,x=x*x%p)

		if(k&1) t=t*x%p;

	return (int)t;

}

bool Miller_Rabin(int x)

{

	if(x==2) return 1;

	if(!(x&1)||x==1) return 0;

	for(int i=0; i<7; ++i)

		if(x==P[i]) return 1;

		else if(!(x%P[i])) return 0;

	int u=x-1,t=0;

	while(!(u&1)) u>>=1,++t;

	for(int i=0; i<7; ++i)

	{

		LL now=FP(P[i],u,x),las;

		for(int j=1; j<=t; ++j)

		{

			las=now, (now*=now)%=x;

			if(now==1&&las!=1&&las!=x-1) return 0;

		}

		if(now!=1) return 0;

	}

	return 1;

}

bool SPFA()

{

	for(int i=1; i<=des; ++i) dis[i]=-1e12;

	dis[src]=0, q.push(src);

	while(!q.empty())

	{

		int x=q.front(); q.pop();

		inq[x]=0;

		for(int v,i=H[x]; i; i=nxt[i])

			if(cap[i]&&dis[to[i]]<dis[x]+cost[i])

			{

				dis[v=to[i]]=dis[x]+cost[i], pre[v]=i;

				if(!inq[v]) q.push(v),inq[v]=1;

			}

	}

	return dis[des]>-1e12;

}

LL Solve()

{

	int flow=0,mn;

	LL sum_cost=0;

	while(SPFA())

	{//注意这只是一个二分图。。

		mn=1e9;

		for(int i=des; i!=src; i=fr[pre[i]]) mn=std::min(cap[pre[i]],mn);

		if(sum_cost+1ll*dis[des]*mn>=0){//dis[des]就是费用。

			for(int i=des; i!=src; i=fr[pre[i]]) cap[pre[i]]-=mn,cap[pre[i]^1]+=mn;

			sum_cost+=1ll*dis[des]*mn, flow+=mn;

		}

		else{

			flow+=(int)(sum_cost/std::abs(dis[des]));//加上剩余正价值能补多少负值。

			break;

		}

	}

	return flow>>1;

}

int main()

{

	n=read(),Enum=1,src=0,des=n<<1|1;

	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();

	for(int i=1; i<=n; ++i) AddEdge(src,i,B[i],0),AddEdge(i+n,des,B[i],0);

	for(int i=1; i<=n; ++i)

		for(int j=1; j<=n; ++j)

			if(A[i]>A[j]&&!(A[i]%A[j]))

			{

				int t=A[i]/A[j];

				if(is_p[0].count(t)) continue;

				if(is_p[1].count(t)) AddEdge(i,j+n,INF,1ll*C[i]*C[j]),AddEdge(j,i+n,INF,1ll*C[i]*C[j]);

				else{

					bool f=Miller_Rabin(t);

					if(f) is_p[1].insert(t),AddEdge(i,j+n,INF,1ll*C[i]*C[j]),AddEdge(j,i+n,INF,1ll*C[i]*C[j]);

					else is_p[0].insert(t);

				}

			}

	printf("%lld",Solve());

	return 0;

}

第二种方法（19.2.17）：

/*

2840kb	28ms

常规做建二分图跑费用流很容易，但因为拆了点每对数之间可能选两次。但是既然选了一次那就一定会选第二次。最后费用即将变成负数时判一下还能选多少个，然后流量除以二即可。

因为匹配条件是$i$可整除$j$且结果是质数，那么$i,j$之间质因子次数一定相差一，可以直接判，同时根据质因子次数奇偶性可以分别连源点汇点，就没有之前的问题了。

另外这张二分图里费用就是$dis[T]$。

*/

#include <queue>

#include <cstdio>

#include <cctype>

#include <algorithm>

#define gc() getchar()

typedef long long LL;

const int N=205,M=(N*N+N)*2,INF=1<<30;

const LL INFll=1ll<<60;

int T,Enum,H[N],nxt[M],fr[M],to[M],cap[M],pre[N];

LL dis[N],cost[M];

inline int read()

{

	int now=0,f=1;register char c=gc();

	for(;!isdigit(c);c=='-'&&(f=-1),c=gc());

	for(;isdigit(c);now=now*10+c-48,c=gc());

	return now*f;

}

inline void AE(int u,int v,int w,LL c)

{

	to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, cap[Enum]=w, cost[Enum]=c;

	to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, cap[Enum]=0, cost[Enum]=-c;

}

int Calc(int x)

{

	int t=0;

	for(int i=2; 1ll*i*i<=x; ++i)

		if(!(x%i))

		{

			x/=i, ++t;

			while(!(x%i)) x/=i, ++t;

		}

	if(x!=1) ++t;

	return t;

}

bool SPFA()

{

	static bool inq[N];

	static std::queue<int> q;

	for(int i=1; i<=T; ++i) dis[i]=-INFll;

	q.push(0), dis[0]=0;

	while(!q.empty())

	{

		int x=q.front(); q.pop();

		inq[x]=0;

		for(int i=H[x],v; i; i=nxt[i])

			if(dis[v=to[i]]<dis[x]+cost[i] && cap[i])

				pre[v]=i, dis[v]=dis[x]+cost[i], !inq[v]&&(q.push(v),inq[v]=1);

	}

	return dis[T]>-INFll;

}

int Solve()

{

	int res=0; LL now=0;

	while(SPFA())

	{

		int mn=INF;

		for(int i=T; i; i=fr[pre[i]]) mn=std::min(mn,cap[pre[i]]);

		if(now+dis[T]*mn>=0)

		{

			now+=dis[T]*mn, res+=mn;

			for(int i=T; i; i=fr[pre[i]]) cap[pre[i]]-=mn, cap[pre[i]^1]+=mn;

		}

		else

		{

			res+=now/std::abs(dis[T]);//abs!!

			break;

		}

	}

	return res;

}

int main()

{

	static int A[N],B[N],C[N],tm[N];

	freopen("pair.in","r",stdin);

	freopen("pair.out","w",stdout);

	int n=read(); Enum=1, T=n+1;

	for(int i=1; i<=n; ++i) tm[i]=Calc(A[i]=read());

	for(int i=1; i<=n; ++i) B[i]=read();

	for(int i=1; i<=n; ++i) C[i]=read();

	for(int i=1; i<=n; ++i)

	{

		if(tm[i]&1)

		{

			AE(0,i,B[i],0);

			for(int j=1; j<=n; ++j)

				if(!(tm[j]&1) && tm[i]-tm[j]==1 && !(A[i]%A[j])) AE(i,j,INF,1ll*C[i]*C[j]);

		}

		else

		{

			AE(i,T,B[i],0);

			for(int j=1; j<=n; ++j)

				if(tm[j]&1 && tm[i]-tm[j]==1 && !(A[i]%A[j])) AE(j,i,INF,1ll*C[i]*C[j]);

		}

	}

	printf("%d\n",Solve());

	return 0;

}

弃疗的暴力+最大流（18.4.1）：

#include <set>

#include <queue>

#include <cstdio>

#include <cctype>

#include <algorithm>

#define gc() getchar()

typedef long long LL;

const int N=410,M=100005,INF=1e9;

const int P[9]={2,3,5,7,13,17,31};

int n,range,src,des,A[N],B[N],C[N],Enum,cur[N],H[N],nxt[M],fr[M],to[M],cap[M],pre[N],lev[N],num[N],que[N];

std::queue<int> q;

std::set<int> is_p[2];

inline int read()

{

	int now=0,f=1;register char c=gc();

	for(;!isdigit(c);c=gc()) if(c=='-') f=-1;

	for(;isdigit(c);now=now*10+c-'0',c=gc());

	return now*f;

}

inline void AddEdge(int u,int v,int w)

{

	printf("%d->%d %d\n",u,v,w);

	to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, cap[Enum]=w;

	to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, cap[Enum]=0;

}

int FP(LL x,int k,int p)

{

	LL t=1;

	for(; k; k>>=1,x=x*x%p)

		if(k&1) t=t*x%p;

	return (int)t;

}

bool Miller_Rabin(int x)

{

	if(x==2) return 1;

	if(!(x&1)||x==1) return 0;

	for(int i=0; i<7; ++i)

		if(x==P[i]) return 1;

		else if(!(x%P[i])) return 0;

	int u=x-1,t=0;

	while(!(u&1)) u>>=1,++t;

	for(int i=0; i<7; ++i)

	{

		LL now=FP(P[i],u,x),las;

		for(int j=1; j<=t; ++j)

		{

			las=now, (now*=now)%=x;

			if(now==1&&las!=1&&las!=x-1) return 0;

		}

		if(now!=1) return 0;

	}

	return 1;

}

namespace Subtask1

{

	int Ans;

	bool can[13][13],vis[13];

	void DFS(int x,int tm,LL sum)

	{//O(n!)

		if(x>n)

			if(sum>=0 && tm>Ans) Ans=tm;

			else ;

		else{

			if(vis[x]) DFS(x+1,tm,sum);

			else{

				vis[x]=1;

				for(int i=1; i<=n; ++i)

					if(can[x][i]&&!vis[i])

						vis[i]=1, DFS(x+1,tm+1,sum+1ll*C[x]*C[i]), vis[i]=0;

				vis[x]=0;

				DFS(x+1,tm,sum);

			}

		}

	}

	void Solve()

	{

		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();

		for(int i=1; i<=n; ++i)

			for(int j=1; j<=n; ++j)

				if(!(A[i]%A[j])&&i!=j)

				{

					int t=A[i]/A[j];

					if(is_p[0].count(t)) continue;

					if(is_p[1].count(t)) can[i][j]=1;

					else{

						bool f=Miller_Rabin(t);

						if(f) is_p[1].insert(t),can[i][j]=1;

						else is_p[0].insert(t);

					}

				}

		DFS(1,0,0);

		printf("%d",Ans);

	}

}

bool BFS()

{

	for(int i=src; i<des; ++i) lev[i]=des+1;

	lev[des]=0, que[0]=des; int h=0,t=1;

	while(h<t)

	{

		int x=que[h++];

		for(int i=H[x]; i; i=nxt[i])

			if(lev[to[i]]==des+1 && cap[i^1])

				lev[to[i]]=lev[x]+1, que[t++]=to[i];

	}

	return lev[src]<=des;

}

void Augment(){

	for(int i=des; i!=src; i=fr[pre[i]])

		--cap[pre[i]], ++cap[pre[i]^1];

}

int ISAP()

{

	if(!BFS()) return 0;

	for(int i=src; i<=des; ++i) ++num[lev[i]],cur[i]=H[i];

	int x=src,res=0; bool f=0;

	while(lev[src]<=des)

	{

		if(x==des) x=src,f=0,++res,Augment();

		bool can=0;

		for(int i=cur[x]; i; i=nxt[i])

			if((to[i]==des&&f)||(to[i]!=des &&  cap[i]))

			{

				printf("%d->%d lev:%d f:%d ",x,to[i],lev[x],f);

				if(x>n&&to[i]<=n) f=1;

				can=1, cur[x]=i, pre[x=to[i]]=i;

				printf("now:%d\n",f);

				break;

			}

		if(!can)

		{

			int mn=des;

			for(int i=H[x]; i; i=nxt[i])

				if((to[i]==des&&f)||(to[i]!=des&&cap[i])) mn=std::min(mn,lev[to[i]]);

			if(!--num[lev[x]]) break;

			printf("Change:%d %d->",x,lev[x]);

			++num[lev[x]=mn+1];

			printf("%d\n",lev[x]);

			cur[x]=H[x];

			if(x!=src)

			{

				if(x<=n&&fr[pre[x]]>n) f=0;

				x=fr[pre[x]];

			}

		}

	}

	return res;

}

int main()

{

	freopen("pair.in","r",stdin);

	freopen("pair.out","w",stdout);

	n=read();

	if(n<=10) {Subtask1::Solve(); return 0;}

	Enum=1,src=0,des=n<<1|1;

	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();

	for(int i=1; i<=n; ++i) AddEdge(src,i,INF),AddEdge(i+n,des,INF),AddEdge(i,i+n,B[i]);

	for(int i=1; i<=n; ++i)

		for(int j=1; j<=n; ++j)

			if(!(A[i]%A[j])&&i!=j)

			{

				int t=A[i]/A[j];

				if(is_p[0].count(t)) continue;

				if(is_p[1].count(t)) AddEdge(i+n,j,INF);

				else{

					bool f=Miller_Rabin(t);

					if(f) is_p[1].insert(t),AddEdge(i+n,j,INF);

					else is_p[0].insert(t);

				}

			}

	printf("%d",ISAP());

	return 0;

}

丧心病狂的费用流(其实就是sb了)（18.4.1）

#include <set>

#include <queue>

#include <cmath>

#include <cstdio>

#include <cctype>

#include <algorithm>

#define gc() getchar()

typedef long long LL;

const int N=410,M=100005,INF=1e9;

const int P[9]={2,3,5,7,13,17,31};

int n,range,src,des,A[N],B[N],C[N],Enum,H[N],nxt[M],fr[M],to[M],cap[M],Cap[M],pre[N];

double dis[N],cost[M],neg[M];

bool inq[N],ignore[M];

std::queue<int> q;

std::set<int> is_p[2];

inline int read()

{

	int now=0,f=1;register char c=gc();

	for(;!isdigit(c);c=gc()) if(c=='-') f=-1;

	for(;isdigit(c);now=now*10+c-'0',c=gc());

	return now*f;

}

inline void AddEdge(int u,int v,int w,int c)

{

	if(!c){

		to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, Cap[Enum]=cap[Enum]=w, cost[Enum]=0, ignore[Enum]=1;

		to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, Cap[Enum]=cap[Enum]=0, cost[Enum]=0, ignore[Enum]=1;

	}

	else if(c<0){

		to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, Cap[Enum]=cap[Enum]=w, cost[Enum]=-log2((double)(-c)), neg[Enum]=1;

		to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, Cap[Enum]=cap[Enum]=0, cost[Enum]=log2((double)(-c)), neg[Enum]=1;

	}

	else{

		to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, Cap[Enum]=cap[Enum]=w, cost[Enum]=log2((double)c);

		to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, Cap[Enum]=cap[Enum]=0, cost[Enum]=-log2((double)c);

	}

	printf("%d->%d %d %d %.3lf\n",u,v,w,c,cost[Enum-1]);

}

int FP(LL x,int k,int p)

{

	LL t=1;

	for(; k; k>>=1,x=x*x%p)

		if(k&1) t=t*x%p;

	return (int)t;

}

bool Miller_Rabin(int x)

{

	if(x==2) return 1;

	if(!(x&1)||x==1) return 0;

	for(int i=0; i<7; ++i)

		if(x==P[i]) return 1;

		else if(!(x%P[i])) return 0;

	int u=x-1,t=0;

	while(!(u&1)) u>>=1,++t;

	for(int i=0; i<7; ++i)

	{

		LL now=FP(P[i],u,x),las;

		for(int j=1; j<=t; ++j)

		{

			las=now, (now*=now)%=x;

			if(now==1&&las!=1&&las!=x-1) return 0;

		}

		if(now!=1) return 0;

	}

	return 1;

}

bool SPFA()

{

	for(int i=1; i<=des; ++i) dis[i]=1e14;

	dis[src]=0, q.push(src);

	while(!q.empty())

	{

		int x=q.front(); q.pop();

		inq[x]=0;

		for(int v,i=H[x]; i; i=nxt[i])

			if(cap[i]&&(neg[i]||dis[to[i]]>dis[x]+cost[i]))

			{

				dis[v=to[i]]=dis[x]+cost[i], pre[v]=i;

				if(!inq[v]) q.push(v);

			}

	}

	return dis[des]<1e14;

}

bool debug=1;

double MCMF(LL &flow)

{

	int mn=1e9; double res=0,tmp1=0;

	bool f=0;

	for(int i=des; i!=src; i=fr[pre[i]])

	{

		mn=std::min(cap[pre[i]],mn);

		if(ignore[i]&&i<=range) f=1;

	}

//	printf("flow:%I64d ",flow);

	if(mn>=flow) mn=flow,flow=0;

	else flow-=mn;

	if(f){

		for(int i=des; i!=src; i=fr[pre[i]])

		{

			if(debug) printf("%d<-",i),

			cap[pre[i]]-=mn, cap[pre[i]^1]+=mn;

			if(neg[i]) ;

		}

		if(debug) printf("%d\n",src);

		if(debug) printf("then:%I64d %d\n",flow,mn);

		return 0;

	}

	for(int i=des; i!=src; i=fr[pre[i]])

	{

		if(debug) printf("%d<-",i),

		cap[pre[i]]-=mn, cap[pre[i]^1]+=mn, res+=cost[pre[i]];

		if(neg[pre[i]]&&!cost[pre[i]]) tmp1+=1.0;

	}

	if(debug) printf("%d\n",src);

	if(debug) printf("then:%I64d res:%.3lf 2:%.3lf\n",flow,res,res<0 ? -1.0*mn*pow(2.0,-res) : 1.0*mn*pow(2.0,res));

	return res<0 ? -1.0*mn*pow(2.0,-res)-tmp1 : 1.0*mn*pow(2.0,res)-tmp1;

}

bool Check(LL flow)

{

	double sigma=0;

	LL tmp=flow;

	for(int i=2; i<=Enum; ++i) cap[i]=Cap[i];

	while(SPFA())

	{

		sigma+=MCMF(flow);

		if(!flow||sigma>0) break;

	}

	if(debug) printf("End:flow:%I64d bef:%I64d sigma:%.3lf\n\n",flow,tmp,sigma);

	return !flow&&sigma<=0;

}

int main()

{

	freopen("pair.in","r",stdin);

//	freopen("pair.out","w",stdout);

	n=read(),Enum=1,src=0,des=n<<1|1;

	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();

	for(int i=1; i<=n; ++i) AddEdge(src,i,INF,0),AddEdge(i+n,des,INF,0),AddEdge(i,i+n,B[i],-C[i]);

//	range=Enum;

	for(int i=1; i<=n; ++i)

		for(int j=1; j<=n; ++j)

			if(!(A[i]%A[j])&&i!=j)

			{

				int t=A[i]/A[j];

				if(is_p[0].count(t)) continue;

				if(is_p[1].count(t)) AddEdge(i+n,j,INF,0);

				else{

					bool f=Miller_Rabin(t);

					if(f) is_p[1].insert(t),AddEdge(i+n,j,INF,0);

					else is_p[0].insert(t);

				}

			}

//	long long l=0,r=1e10,mid,ans=0;

	long long l=1,r=15,mid,ans=0;

	while(l<=r)

		if(Check(mid=l+r>>1)) ans=mid,l=mid+1;

		else r=mid-1;

	printf("%I64d",ans);

	return 0;

}

相关文章