Project Euler Problem 675

ORZ foreverlasting聚聚，QQ上问了他好久才会了这题（所以我们又聊了好久的Gal）

我们先来尝试推导一下\(S\)的性质，我们利用狄利克雷卷积来推：

\[2^\omega=I\ast|\mu|\]

这个很好理解吧，考虑一下它的组合意义即可

然后两边同卷上\(I\)有：

\[2^\omega \ast I=I\ast I\ast |\mu|=d\ast |\mu|\]

后面还是同样，考虑\(d\ast |\mu|\)的组合意义，一正一反的情况下其实就是\(d(n^2)\)

因此我们有了\(2^\omega\ast I=d(n^2)\)，即\(S(n)=d(n^2)\)

那么显然\(S\)现在是个积性函数了，答案又是阶乘的形式，因此可以从\(n-1\)的答案推到\(n\)来

考虑一个非常暴力的过程，每次暴力分解质因数，复杂度大概是\(O(n\sqrt n)\)的

然后你只需要一台好一点的电脑我仿佛已经闻到了CPU的香气

然后考虑怎么优化这个过程，我们发现类似于某个套路，这种方法之所以慢是因为会出现不必要的枚举，因此我们只需要记录一下每个数的最小质因数，然后每次直接除去即可，顺带把贡献算一下

这样的复杂度很迷啊，加藤聚聚说是一个\(\log\)的，我感觉还要再少点，毕竟向下除至少去掉一个\(2\)

那么我们就可以很快的做掉这道题了（用自己的笔记本跑了2s就出来了）

#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=10000000,mod=1000000087;
int prime[N+5],cnt,mnp[N+5],bkt[N+5],inv[(N<<1)+5],ret=1,ans;
#define Pi prime[j]
inline void init(void)
{
    RI i,j; for (mnp[1]=1,i=2;i<=N;++i)
    {
        if (!mnp[i]) mnp[i]=i,prime[++cnt]=i;
        for (j=1;j<=cnt&&1LL*i*Pi<=N;++j)
        {
            mnp[i*Pi]=Pi; if (i%Pi==0) break;
        }
    }
    for (inv[0]=inv[1]=1,i=2;i<=(N<<1)+1;++i)
    inv[i]=1LL*inv[mod%i]*(mod-mod/i)%mod;
}
#undef Pi
inline void inc(int& x,CI y)
{
    if ((x+=y)>=mod) x-=mod;
}
inline int sum(CI x,CI y)
{
    int t=x+y; return t>=mod?t-mod:t;
}
int main()
{
    freopen("ans.txt","w",stdout);
    init(); for (RI i=2;i<=N;++i)
    {
        ret=1LL*ret*inv[bkt[i]+1]%mod; inc(bkt[i],2); ret=1LL*ret*(bkt[i]+1)%mod;
        for (int x=i;x!=mnp[x];x/=mnp[x])
        {
            if (x/mnp[x]==mnp[x])
            {
                ret=1LL*ret*inv[bkt[x]+1]%mod*inv[bkt[mnp[x]]+1]%mod;
                inc(bkt[mnp[x]],sum(bkt[x],bkt[x]));
                ret=1LL*ret*(bkt[mnp[x]]+1)%mod; bkt[x]=0;
            } else
            {
                ret=1LL*ret*inv[bkt[x]+1]%mod*inv[bkt[mnp[x]]+1]%mod*inv[bkt[x/mnp[x]]+1]%mod;
                inc(bkt[mnp[x]],bkt[x]); inc(bkt[x/mnp[x]],bkt[x]);
                ret=1LL*ret*(bkt[mnp[x]]+1)%mod*(bkt[x/mnp[x]]+1)%mod; bkt[x]=0;
            }
        }
        inc(ans,ret);
    }
    return printf("%d",ans),0;
}

原文地址：https://www.cnblogs.com/cjjsb/p/11518890.html