round15

F:问有多少对n个元素的集合A,B满足它们的笛卡尔和恰由1 ~ n2构成，其中A包含0，B包含1
n<=1012,T<=5000

分析：考虑固定B，如何去确定A;
首先，0在A中，那么所有B中的元素就都在笛卡尔和中；考虑第一个B中没有出现的数x，那么x-1必然在A中（不然1和这个数作用就会重复），然后我们把x-1和B中所有元素再做一遍笛卡尔和；
接着我们再找到第一个和中没有出现过的数x，x-1必然在A中，依此类推，直到穷尽1~n2
从中可以看出，一个B最多对应一个A；于是我们现在的目标就是找到有多少个合法的B；
再观察B有什么性质
假设1~x都在B中，而x+1不在B中，那么x必然在A中，因而可以推断出x+1~2x都不在B中；进一步k*x+1~(k+1)*x要么都在B中，要么都不在
另一方面，假设2～x都不在B中，而x+1在B中，与上面同理可以得到只有kx+1可能在B中
按照这个思路，我们可以令f(n,m)代表n个数，填满1~n*m，且2必须在B中的方案数；g(n,m)代表n个数，填满1~n*m，且2必须不在B中的方案数；
那么有

f(n,m)=∑d|n,d>1g(nd,m)

g(n,m)=∑d|m,d>1f(n,md)

一个惊人的关系是f(n,m)=g(m,n),可以用数学归纳法简单证明
由此可以导出

f(n,m)=∑d|n,d>1f(m,nd)

而我们要算的答案即h(n,n)=f(n,n)+g(n,n)=2f(n,n)
鉴于数据范围巨大，当前的式子还不足以解决这个问题；另一个观察是f(n,m)只与n,m的因子构成有关；因此因子构成相同的可以只算一次；用dfs爆搜所有不同的因子构成，发现只有<5000个状态；因而可以把这5000个答案本地打表出来，这样复杂度就只剩下因式分解的复杂度；为了本地打表我利用了推导出的另一个公式

h(n,m)=∑d|n,d>1−μ(d)h(nd,m)+∑d|m,d>1−μ(d)h(n,md)

由于μ函数的特殊性，我们可以预处理出一张转移表来加速；另外h(n,m)=h(m,n)，因而我们可以只算一半；这样使得总计算次数在109左右，可以跑出结果；
为了节省时间，我采用的是<100w预处理，>100wPollard-rho分解的方法；但我发现仍然很慢TAT

总结:这题我想了很久，就是最后没有想到因子同构QAQ

#include<stdio.h>#include<algorithm>#include<iostream>#include<map>#include<vector>#include<string.h>#include<math.h>#include<assert.h>using namespace std;typedef long long LL;const int Maxn=1000002,Maxstate=5000;const LL M=10000000061LL,totlim=1e12+2;typedef pair<int,int>pi;LL a[Maxn];map<pi,LL>mpF,mpG;LL calF(int n,int m);LL tot=0;map<vector<int>,int >Id;vector<int>All[5000];vector<pi>G[5000];LL dp[100][100];bool isp[Maxn];int least[Maxn];vector<int>pri;void getp(){    for(int i=2;i<Maxn;i++){        if(!isp[i])pri.push_back(i),least[i]=i;        for(int j=0;j<pri.size()&&pri[j]*i<Maxn;j++){            isp[pri[j]*i]=1;            least[pri[j]*i]=pri[j];            if(i%pri[j]==0){break;}        }    }}vector<LL>yinzi;LL mul(LL a,LL b,LL n){return (a*b-(LL)(a/(long double)n*b+1e-3)*n+n)%n;}LL powmod(LL x,LL y,LL mod){    LL ret=1;    x%=mod;    while(y){        if(y&1)ret=mul(ret,x,mod);        y>>=1;        x=mul(x,x,mod);    }    return ret;}bool check(LL a,LL n){    LL m=n-1;    int s=0;    while(~m&1)m>>=1,s++;    LL x=powmod(a,m,n);    for(int i=1;i<=s;i++){        LL y=mul(x,x,n);        if(y==1&&x!=1&&x!=n-1)return 0;        x=y;    }    return x==1;}LL Miller_Rabin(LL x){    if(x==1)return 0;    if(x==2)return 1;    if(~x&1)return 0;    int times=3;//Pro=1/4^times    while(times--)if(!check(rand()%(x-2)+2,x))return 0;    return 1;}LL Pollard_rho(LL n,int c){    LL factor=1;    int cirsize=1;    LL x=rand()%n,xfixed=x;    while(factor==1){        for(int i=0;i<cirsize;i++){            x=(mul(x,x,n)+c)%n;            if(x==xfixed)return 1;            LL d=__gcd(n,(x-xfixed+n)%n);//            if(d>1&&d<n)return d;        }        cirsize<<=1;        xfixed=x;    }    assert(0);}void findfac(LL x,int c){    if(x==1)return;    if(x<Maxn){        yinzi.push_back(least[x]);        findfac(x/least[x],c);        return;    }    if(Miller_Rabin(x)){//ispri        yinzi.push_back(x);        return;    }    LL tmp;    while((tmp=Pollard_rho(x,c--))==1);    findfac(x/tmp,c);findfac(tmp,c);}void dfs(LL n,vector<int>&V){//n<m    All[Id.size()]=V;    int tmp=Id.size();    Id[V]=tmp;    int bef=V.size()?V.back():60;    for(int i=1;i<=bef;i++){        if(n*pri[V.size()]>=totlim)break;        n*=pri[V.size()];        V.push_back(i);        dfs(n,V);        V.pop_back();    }}void precal(){    for(int i=0;i<Id.size();i++){        for(int mask=1;mask<1<<All[i].size();mask++){            vector<int>t=All[i];            for(int k=0;k<t.size();k++){                if(mask>>k&1){                    t[k]--;                 }            }            sort(t.begin(),t.end(),greater<int>());            while(t.size()&&!t.back())t.pop_back();            G[i].push_back(pi(Id[t],__builtin_popcount(mask)&1?1:-1));        }    }    LL tot=0;    for(int i=0;i<Id.size();i++){        for(int j=i;j<Id.size();j++){            if(i==0&&j==0)dp[i][j]=1;            else dp[i][j]=0;            for(int k=0;k<G[i].size();k++){                int nxt=G[i][k].first,flag=G[i][k].second;                dp[i][j]+=dp[nxt][j]*flag;                dp[i][j]%=M;                tot++;            }            for(int k=0;k<G[j].size();k++){                int nxt=G[j][k].first,flag=G[j][k].second;                dp[i][j]+=dp[min(nxt,i)][max(nxt,i)]*flag;                dp[i][j]%=M;                tot++;            }        }    }    for(int i=0;i<Id.size();i++)printf("%lldLL,",(dp[i][i]+M)%M);}map<LL,vector<int> >have;int main(){    //freopen("mysol.txt","w",stdout);    getp();    vector<int>State;    dfs(1,State);    //precal();    int cas=1;    //return 0;    int _;scanf("%d",&_);    while(_--){        LL x;scanf("%lld",&x);        vector<int>V;        if(have.find(x)!=have.end())V=have[x];        else{            yinzi.clear();            findfac(x,2730);                    sort(yinzi.begin(),yinzi.end());            for(int i=0,j;i<yinzi.size();i=j){                for(j=i+1;j<yinzi.size()&&yinzi[i]==yinzi[j];j++);                V.push_back(j-i);            }            sort(V.begin(),V.end(),greater<int>());            have[x]=V;        }        printf("Case #%d: %lld\n",cas++,sol[Id[V]]);    }}