[BZOJ3512]DZY Loves Math IV 杜教筛+记忆化搜索

最新推荐文章于 2020-09-11 18:43:43 发布

原创最新推荐文章于 2020-09-11 18:43:43 发布 · 438 阅读

1 ·

CC 4.0 BY-SA版权

dp 同时被 3 个专栏收录

71 篇文章

订阅专栏

数论

30 篇文章

订阅专栏

杜教筛

7 篇文章

订阅专栏

本文探讨了在数论中如何高效地计算特定形式的求和问题，通过将整数分解为质因数的乘积并利用积性函数的性质，提出了一种优化算法。该算法能够有效减少计算复杂度，特别适用于求解形如S(n,m)的表达式。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

发现 $n$ 只有 $=10^5$ ，我们先考虑 $S(n,m)=\sum_{i=1}^m \varphi(ni)$ 怎么求。
先把 $n$ 分成两部分， $n_1=\prod p_i$ ，即所有质因子一次幂的乘积， $n_2=\prod p_i^{c_i-1}$ ，即剩下的。不难发现 $\varphi(n)=\varphi(n_1)*n_2$ 。

S (n, m) = n 2 \sum i = 1 m φ (n 1 i) = n 2 \sum i = 1 m φ ( n 1 ) φ ( i ) gcd ( n 1 , i ) φ ( gcd ( n 1 , i ) )

$S(n,m)=n_2\sum_{i=1}^m\varphi(n_1i)=n_2\sum_{i=1}^m\frac{\varphi(n_1)\varphi(i)\gcd(n_1,i)}{\varphi(\gcd(n_1,i))}$
然后因为

φφ $\varphi$ 是积性函数，所以当

gcd(b,a/b)=1gcd(b,a/b)=1 $\gcd(b,a/b)=1$ 时，

φ(a)φ(b)=φ(ab)φ(a)φ(b)=φ(ab) $\frac{\varphi(a)}{\varphi(b)}=\varphi(\frac{a}{b})$ ，所以

= n 2 \sum i = 1 m φ (n 1 gcd ( n 1 , i )) φ (i) gcd (n 1, i) = n 2 \sum i = 1 m φ (n 1 gcd ( n 1 , i )) φ (i) \sum d | gcd (n 1, i) φ (d)

$=n_2\sum_{i=1}^m\varphi(\frac{n_1}{\gcd(n_1,i)})\varphi(i)\gcd(n_1,i)=n_2\sum_{i=1}^m\varphi(\frac{n_1}{\gcd(n_1,i)})\varphi(i)\sum_{d|\gcd(n_1,i)}\varphi(d)$
这里用

∑d|nφ(d)∑d|nφ(d) $\sum_{d|n}\varphi(d)$ 代替

nn $n$ 是为了去掉

gcd

$\gcd$ 的限制，并且和前面的合并。

= n 2 \sum i = 1 m φ (i) \sum d | gcd (n 1, i) φ (n 1 d) = n 2 \sum d | n 1 φ (n 1 d) \sum i = 1 ⌊ m d ⌋ φ (d i) = n 2 \sum d | n 1 φ (n 1 d) S (d, m d)

$=n_2\sum_{i=1}^m\varphi(i)\sum_{d|\gcd(n_1,i)}\varphi(\frac{n_1}{d})=n_2\sum_{d|n_1}\varphi(\frac{n_1}{d})\sum_{i=1}^{\lfloor \frac{m}{d}\rfloor}\varphi(di)=n_2\sum_{d|n_1}\varphi(\frac{n_1}{d})S(d,\frac{m}{d})$
注意到

SS $S$ 中的参数最多只有

n \sqrt{m}

$n\sqrt{m}$ 种，可以记忆化，当

n=1n=1 $n=1$ 的时候就是杜教筛啦。
听说复杂度

O(nm−−√+m23)O(nm+m23) $O(n\sqrt{m}+m^{\frac{2}{3}})$ 。可是还是不太明白求了所有的

∑mdi=1φ(i)∑i=1mdφ(i) $\sum_{i=1}^{\frac{m}{d}}\varphi(i)$ 为什么还是

O(m23)O(m23) $O(m^{\frac{2}{3}})$ 。。。
代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<cmath>
#include<vector>
#define MP make_pair
#define PB push_back
#define fs first
#define sc second
#define ll long long
#define N 100010
#define U 6000000
#define P 1000000
#define ID(x,y) (x*1234567891+y)
#define up(x,y) x=(x+(y))%mod
using namespace std;
const int mod=1000000007;
bool flag[U+5];
int n,m,pri[U+5],minp[U+5],num,sp[U+5],phi[U+5];
struct hash
{
    vector<pair<ll,ll> > a[P];
    void ins(ll x,ll y)
    {
        int d=x%P;
        a[d].PB(MP(x,y));
    }
    ll qry(ll x)
    {
        int d=x%P;
        for(int i=a[d].size()-1;i>=0;i--)
            if(a[d][i].fs==x) return a[d][i].sc;
        return -1;   
    }
}h;
void getpri()
{
    memset(flag,1,sizeof(flag));
    flag[1]=0;phi[1]=1;
    for(int i=2;i<=U;i++)
    {
        if(flag[i]) pri[++num]=i,phi[i]=i-1,minp[i]=i;
        for(int j=1;j<=num&&i*pri[j]<=U;j++)
        {
            int v=i*pri[j];
            flag[v]=0;minp[v]=pri[j];
            if(i%pri[j]==0) {phi[v]=phi[i]*pri[j];break;}
            phi[v]=phi[i]*phi[pri[j]];          
        }
    }
}
void fj(int x,int *t)
{
    while(x>1)
    {
        t[++t[0]]=minp[x];
        while(x%t[t[0]]==0) x/=t[t[0]];
    }
}
ll qphi(int x)
{
    if(x<=U) return sp[x];
    ll re=((ll)x*(x+1)/2)%mod;
    for(int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        up(re,mod-qphi(x/l)*(r-l+1)%mod);
    }
    return re;
}
ll S(ll a,ll b)
{
    if(!b) return 0;
    if(h.qry(ID(a,b))!=-1) return h.qry(ID(a,b));
    if(a==1)
    {
        h.ins(ID(a,b),qphi(b));
        return h.qry(ID(a,b));
    }
    int t[12],w=1;
    ll re=0;
    t[0]=0;
    fj(a,t);
    for(int i=1;i<=t[0];i++)
        w*=t[i];
    for(int i=floor(sqrt(w));i;i--)
        if(w%i==0) up(re,S(i,b/i)*phi[w/i]),up(re,S(w/i,b/(w/i))*phi[i]);
    re=re*a/w%mod;
    h.ins(ID(a,b),re);
    return re;  
}
int main()
{
    getpri();
    for(int i=1;i<=U;i++)
        sp[i]=(phi[i]+sp[i-1])%mod; 
    scanf("%d%d",&n,&m);
    ll ans=0;
    for(int i=1;i<=n;i++)
        up(ans,S(i,m));
    printf("%lld",ans);
    return 0;
}