[fzu2016]How many tuples 解题报告

这题tm什么鬼啊。。$10^{10}$的做法竟然500ms ac。
而且更奇怪的是。。似乎别人都打了一个表，他们都是怎么做的呢？

直接裸反演：

$$\sum_{i=1}\mu (i) \prod_{j=1}^m\lfloor {a_j \over i} \rfloor$$
这样如果不考虑求 $s(n)=\sum_{i=1}^n \mu (i)$的话，是 $O(m\sqrt A \log m)(A=10^8)$的。
求 $s(n)$的话，可以用杜教筛。先线筛出 $\mu(i),i\le P,P=10^7$。然后有如下式子： $$\sum_{i=1}^n\sum_{j|i}\mu(j)=1$$ ,枚举 $i\over j$，会得到 $$\sum_{i=1}^n\sum_{j=1}^{\lfloor {n\over i}\rfloor}\mu (j)=\sum_{i=1}^ns(\lfloor {n\over i} \rfloor )=1$$ ，即 $$s(n)=1-\sum_{i=2}^ns(\lfloor{n\over i}\rfloor)$$
按照这个递推式处理n>P的情况，注意到对于每一个a，n>P的情况只会有 $A\over P$种，所以拿map存就好了。这部分的时间复杂度是 $O(\sum_{i=1}^{A\over P}\sqrt {n\over i}+({A\over P})^2\log Tm{A\over P})=O(\sqrt {n{A\over P}}+({A\over P})^2\log Tm{A\over P})$
所以总时间复杂度是 $O(T(m\sqrt A \log m+\sqrt {n{A\over P}}+({A\over P})^2\log Tm{A\over P}))≈10^{10}$。

复杂度太高了。
然后我大为困惑。。想了一整天该怎么做，把式子各种乱化，或者试图不用反演做。。然后感觉实在不会做，各种搜题解，搜不到。。然后发现出题人叫AekdyCoin，似乎数学非常神，人称“福大核武，景润后人”，怪不得这题这么神。。。然后我发现vjudge a了这个题，就去vjudge注册了个账号看看能不能看代码，发现a的人都没共享代码。。。然后我绝望之中不知道怎么回事搜到了拉马努金，听说他没接受过系统教育，他的学习就是用自己的方法证明一本子的定理，深受启发，决定用自己的方法搞定这道题。于是又想了一上午，却毫无进展。。这时我突然发现ac status里面平均码长有50k，最短的也有10k，所以这肯定是道打表题！所以我试图各种打表，搞了一下午，并没有什么结果。。绝望之中我感到我比拉马努金差了几万个陈景润，是不可能像他一样了，于是想找个大神问一问，就去知乎上随便翻一翻，发现叉姐好像是愿意帮新人解答问题的那种，就想去找下怎么联系叉姐。先去uoj群里看了看，没有。于是我打算去叉姐icpc-camp的那个网站发个帖什么的，但是注册需要邮箱验证，但是我的邮箱并没有收到邮件。直接登录显示没有此用户，再注册显示邮箱已经被使用。。然后我发现重设密码也会发邮件，于是点了重设。毫无卵用。。搞了半天之后我决定换个邮箱试试，也许是qq太弱，网易说不定就可以。。结果网易还是不行，并没有邮件，然后点了重设密码，本来已经不抱希望了。。结果竟然收到邮件了！！妈蛋终于注册上了！！！然后跟叉姐问题，她竟然秒回了！瞬间好感动！然而她说她觉得这种比赛不可能全是极限数据。。我感觉非常不满意。。但是也没办法了，既然叉姐都这么说了，我就写写试试看吧。。没想到还真a了！。。

代码：

#include<cstdio>
#include<iostream>
using namespace std;
#include<algorithm>
#include<map>
const int T=1e4+5,M=20+5,A=1e8+5;
struct HS{
    int a,r;
}heap[M];
int m;
void down(int node){
    for(int next=node<<1;next<=m;node=next,next<<=1){
        if(next<m&&heap[next+1].r<heap[next].r)++next;
        if(heap[node].r<heap[next].r)return;
        swap(heap[node],heap[next]);
    }
}

const int B=2e4;
//const int B=2;
typedef long long LL;
const int Mod=1e9+7;
LL power(LL prod,int x){
    LL ans=1;
    for(;x;x>>=1){
        if(x&1)ans=ans*prod%Mod;
        prod=prod*prod%Mod;
    }
    return ans;
}
LL fac[100000+5],rev[100000+5];
const int P=1e7;
int prime[1000005];
bool p[10000005];
int smu[10000005];
map<int,int> Map;
int query(int n){
    if(n<=P)return smu[n];
    else if(Map[n])return Map[n];
    else{
        int ans=1;
        for(int i=2,j;i<=n;i=j+1){
            j=n/(n/i);
            ans-=(LL)(j-i+1)*query(n/i);
        }
        return Map[n]=ans;
    }
}

int a[M];
int main(){
    freopen("fzu_2106.in","r",stdin);
    freopen("fzu_2106.out","w",stdout);

    fac[0]=1;
    for(int i=1;i<=100000;++i)fac[i]=fac[i-1]*i%Mod;
    LL rfac=power(fac[100000],Mod-2);
    for(int i=99999;i>=0;--i){
        rev[i+1]=fac[i]*rfac%Mod;
        rfac=rfac*(i+1)%Mod;
    }
    smu[1]=1;
    for(int i=2;i<=P;++i){
        if(!p[i]){
            prime[++prime[0]]=i;
            smu[i]=-1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=P;++j){
            p[i*prime[j]]=1;
            if(i%prime[j])smu[i*prime[j]]=-smu[i];
            else break;
        }
    }
    for(int i=2;i<=P;++i)smu[i]+=smu[i-1];

    int t;
    scanf("%d",&t);
    while(t--){
        scanf("%d",&m);
        for(int i=1;i<=m;++i)scanf("%d",a+i);
        sort(a+1,a+m+1);
        LL ans=0,prod;
        if(a[1]<=B){
            for(int i=a[1];i;--i)
                if(smu[i]!=smu[i-1]){
                    prod=1;
                    for(int j=m;j;--j)prod=prod*(a[j]/i)%Mod;
                    ans=(ans+(smu[i]-smu[i-1])*prod)%Mod;
                }
        }
        else{
            for(int i=B;i;--i)
                if(smu[i]!=smu[i-1]){
                    prod=1;
                    for(int j=m;j;--j)prod=prod*(a[j]/i)%Mod;
                    ans=(ans+(smu[i]-smu[i-1])*prod)%Mod;
                }
            prod=1;
            for(int j=m;j;--j)prod=prod*(a[j]/(B+1))%Mod;
            for(int j=m;j;--j){
                heap[j]=(HS){a[j],a[j]/(a[j]/(B+1))};
                down(j);
            }
            bool flag=1;
            for(int i,tmp=B;flag;tmp=i){
                i=heap[1].r;
                ans=(ans+(query(i)-query(tmp))*prod)%Mod;
                for(;heap[1].r==i;down(1)){
                    if(i==heap[1].a){
                        flag=0;
                        break;
                    }
                    prod=prod*rev[heap[1].a/i]%Mod*(heap[1].a/(i+1))%Mod;
                    heap[1].r=heap[1].a/(heap[1].a/(i+1));
                }
            }
        }
        printf("%I64d\n",(ans+Mod)%Mod);
    }
}

总结：
①一定要检查数组大小！！