[BZOJ4514][Sdoi2016]数字配对（数论+费用流）

题目描述

传送门

题解

首先判断两个数是否能配对
线筛$10^7$范围内的质数，然后超过这个范围的用Miller_Rabbin——直接暴力也可以
拆点xiyi，对于能配对的两个点ij，连边xi->yj,inf,ci*cj，xj->yi,inf,ci*cj
然后对于每一个点i，连边s->xi,bi,0，yi->t,bi,0
这样建图之后跑最大费用费用流，因为是每一次找一条最大的费用增广，所以如果当前的总费用不满足非负就直接推出。注意这里最后一次增广不一定满流，所以要计算一下在非负的条件下最多增广多大的流
得出的最大流/2就是答案，如果为奇数需要下取整，因为只有跑出来对称的两次才能算是完整的配对。

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;
#define LL long long
#define N 410
#define inf 2100000000
#define Abs(x) ((x>0)?x:-x)

int n,s,t,maxflow;LL disc,mincost;
int a[N],b[N];LL c[N];
int p[10000005],prime[5000005];
int tot,point[N],nxt[N*N*2],v[N*N*2],remain[N*N*2];LL u[N*N*2];
LL dis[N];int last[N];bool vis[N];
queue <int> q;

void get(int n)
{
    for (int i=2;i<=n;++i)
    {
        if (!p[i]) prime[++prime[0]]=i;
        for (int j=1;j<=prime[0]&&i*prime[j]<=n;++j)
        {
            p[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
}
LL fast_pow(LL a,int p,LL Mod)
{
    LL ans=1;
    for (;p;p>>=1,a=a*a%Mod)
        if (p&1)
            ans=ans*a%Mod;
    return ans;
}
bool miller_rabbin(int n)
{
    if (n==2) return 1;
    if (n<=1||!(n&1)) return 0;
    int t=0,a,x,y,u=n-1;
    while (!(u&1)) ++t,u>>=1;
    for (int i=0;i<=10;++i)
    {
        a=rand()%(n-1)+1;
        x=fast_pow(a,u,n);
        for (int j=0;j<t;++j)
        {
            y=(LL)x*x%n;
            if (y==1&&x!=1&&x!=n-1) return 0;
            x=y;
        }
        if (x!=1) return 0;
    }
    return 1;
}
bool check(int x)
{
    if (x==1) return 0;
    if (x<=10000000&&!p[x]) return 1;
    return miller_rabbin(x);
}
void addedge(int x,int y,int cap,LL z)
{
    ++tot; nxt[tot]=point[x]; point[x]=tot; v[tot]=y; remain[tot]=cap; u[tot]=z;
    ++tot; nxt[tot]=point[y]; point[y]=tot; v[tot]=x; remain[tot]=0; u[tot]=-z;
}
int addflow(int s,int t)
{
    int now=t,ans=inf;
    while (now!=s)
    {
        ans=min(ans,remain[last[now]]);
        now=v[last[now]^1];
    }
    now=t;
    while (now!=s)
    {
        remain[last[now]]-=ans;
        remain[last[now]^1]+=ans;
        now=v[last[now]^1];
    }
    return ans;
}
bool spfa(int s,int t)
{
    memset(dis,128,sizeof(dis));disc=dis[s];dis[s]=0;
    memset(vis,0,sizeof(vis));vis[s]=1;
    q.push(s);
    while (!q.empty())
    {
        int now=q.front();q.pop();
        vis[now]=0;
        for (int i=point[now];i!=-1;i=nxt[i])
            if (dis[v[i]]<dis[now]+u[i]&&remain[i])
            {
                dis[v[i]]=dis[now]+u[i];
                last[v[i]]=i;
                if (!vis[v[i]]) vis[v[i]]=1,q.push(v[i]);
            }
    }
    if (dis[t]==disc) return 0;
    int flow=addflow(s,t);
    if (mincost+(LL)flow*dis[t]<0)
    {
        maxflow+=mincost/Abs(dis[t]);
        return 0;
    }
    maxflow+=flow;
    mincost+=flow*dis[t];
    return 1;
}
int main()
{
    get(10000000);
    scanf("%d",&n);
    for (int i=1;i<=n;++i) scanf("%d",&a[i]);
    for (int i=1;i<=n;++i) scanf("%d",&b[i]);
    for (int i=1;i<=n;++i) scanf("%lld",&c[i]);
    tot=-1;memset(point,-1,sizeof(point));
    s=n+n+1,t=s+1;
    for (int i=1;i<=n;++i) addedge(s,i,b[i],0),addedge(n+i,t,b[i],0);
    for (int i=1;i<=n;++i)
        for (int j=i+1;j<=n;++j)
        {
            int x=a[i],y=a[j];
            if (x<y) swap(x,y);
            if (x%y==0&&check(x/y))
            {
                addedge(i,n+j,inf,c[i]*c[j]);
                addedge(j,n+i,inf,c[i]*c[j]);
            }
        }
    while (spfa(s,t));
    printf("%d\n",maxflow>>1);
}