round15

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

分析：考虑固定B，如何去确定A;
首先，0在A中，那么所有B中的元素就都在笛卡尔和中；考虑第一个B中没有出现的数x，那么x-1必然在A中（不然1和这个数作用就会重复），然后我们把x-1和B中所有元素再做一遍笛卡尔和；
接着我们再找到第一个和中没有出现过的数x，x-1必然在A中，依此类推，直到穷尽1~$n^2$
从中可以看出，一个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) =\sum_{d|n,d>1}g(\frac{n}{d},m)$$
$$g(n,m) =\sum_{d|m,d>1}f(n,\frac{m}{d})$$
一个惊人的关系是f(n,m)=g(m,n),可以用数学归纳法简单证明
由此可以导出
$$f(n,m) =\sum_{d|n,d>1}f(m,\frac{n}{d})$$
而我们要算的答案即 $h(n,n)=f(n,n)+g(n,n)=2f(n,n)$
鉴于数据范围巨大，当前的式子还不足以解决这个问题；另一个观察是f(n,m)只与n,m的因子构成有关；因此因子构成相同的可以只算一次；用dfs爆搜所有不同的因子构成，发现只有<5000个状态；因而可以把这5000个答案本地打表出来，这样复杂度就只剩下因式分解的复杂度；为了本地打表我利用了推导出的另一个公式
$$h(n,m)=\sum_{d|n,d>1}-\mu(d)h(\frac{n}{d},m)+\sum_{d|m,d>1}-\mu(d)h(n,\frac{m}{d})$$
由于 $\mu$函数的特殊性，我们可以预处理出一张转移表来加速；另外h(n,m)=h(m,n)，因而我们可以只算一半；这样使得总计算次数在 $10^9$左右，可以跑出结果；
为了节省时间，我采用的是<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]]);
    }
}