杜教筛

标签：code   etc   直接   splay   prim   for   lse   getch   psi   

还是发出来提醒下自己,免得又给忘了...

杜教筛

• 一种可以在 \(O(n^{\frac 2 3})\) 时间内解决积性函数前缀和的操作.
• 求 \(S(n)=\sum_{i=1}^n f(i),f\) 为积性函数.
• 构造两个积性函数 \(g,h\) ,使得 \(f*g=h\).

\[ \begin{align*} \sum_{i=1}^n h(i)&= \sum_{i=1}^n \sum_{d|n}f(d)*g(\frac i d)\\&= \sum_{d=1}^n g(d) \cdot \sum_{i=1}^{\lfloor \frac n d\rfloor} f(i)\\&= g(1)\cdot S(n)+\sum_{d=2}^n g(d)\cdot S(\lfloor \frac n d \rfloor)\ \therefore\ g(1)\cdot S(n)&=\sum_{i=1}^n h(i)-\sum_{d=2}^n g(d)\cdot S(\lfloor \frac n d \rfloor). \end{align*} \]

• 这样,只要我们选取的 \(g,h\) 的前缀和很好求,计算 \(S(n)\) 就只有 \(O(n^{\frac 2 3})\) 的复杂度.

• 实现时可以暴力算前 \(\sqrt n\) 项的 \(S\) ,按照上面的式子递归计算,并用 \(hash\) 或 \(map\) 记忆化.

• 常用的几个 \(f*g=h\) :
• \(1.\ \mu\ *\ I=\epsilon\).
• \(2.\ \varphi\ *\ I=id\).
• \(3.\ id\ *\ \mu=\varphi\).

• \(4.\ f\ *\epsilon=f\).

• 另外,若 \(f(i)=\mu(i) \cdot i^k,\varphi(i) \cdot i^k\) , 也十分好构造.直接用原来对应的 \(g\) ,这样得到的 \(h\) 会多乘一个 \(n^k\) ,而这个东西的前缀和可以背公式(\(k\) 小)或插值搞出来.(\(k\) 大).

• 例:

简单的数学题

• 题意:求\(\sum_{i=1}^n \sum_{j=1}^n ij*gcd(i,j)\) , \(n\leq 1e10\) ,答案对一个质数 \(p\) 取模.

• 这道题用 \(\varphi\) 反演会简单一些.

\[ \begin{align*} \sum_{i=1}^n \sum_{j=1}^n ij*gcd(i,j)&= \sum_{i=1}^n \sum_{j=1}^n ij \sum_{k|gcd(i,j)} \varphi(k)\\&= \sum_{k=1}^n \varphi(k) \sum_{k|i}\sum_{k|j}ij\\&= \sum_{k=1}^n \varphi(k) \cdot k^2 \cdot (\sum_{i=1}^{\lfloor \frac n k \rfloor} i)^2. \end{align*} \]

• 现在外层整除分块是 \(\sqrt n\) 的,考虑到 \(n\) 的大小,需要快速计算出\(\sum_{k=l}^r \varphi(k) \cdot k^2\) ,需杜教筛解决.

#include"bits/stdc++.h"
using namespace std;
typedef long long LoveLive;
inline LoveLive read()
{
    LoveLive out=0,fh=1;
    char jp=getchar();
    while ((jp>'9'||jp<'0')&&jp!='-')
        jp=getchar();
    if (jp=='-')
        {
            fh=-1;
            jp=getchar();
        }
    while (jp>='0'&&jp<='9')
        {
            out=out*10+jp-'0';
            jp=getchar();
        }
    return out*fh;
}
const int MAXN=4600010;
int p;
LoveLive n;
LoveLive add(LoveLive a,LoveLive b)
{
    return (a + b) % p;
}
LoveLive mul(LoveLive a,LoveLive b)
{
    return (a * b) % p;
}
LoveLive fpow(LoveLive a,LoveLive b)
{
    a%=p;
    b%=p-1;
    int res=1;
    while(b)
        {
            if(b&1)
                res=mul(res,a);
            a=mul(a,a);
            b>>=1;
        }
    return res;
}
int prime[MAXN],cnt=0,ism[MAXN];
LoveLive phi[MAXN],sumphi[MAXN];
LoveLive inv6,inv2;
void init(int N)
{
    inv6=fpow(6,p-2);
    inv2=fpow(2,p-2);
    ism[1]=phi[1]=1;
    for(int i=2;i<=N;++i)
        {
            if(!ism[i])
                prime[++cnt]=i,phi[i]=i-1;
            for(int j=1;j<=cnt && i*prime[j]<=N;++j)
                {
                    ism[i*prime[j]]=1;
                    if(i%prime[j]==0)
                        {
                            phi[i*prime[j]]=phi[i]*prime[j];
                            break;
                        }
                    else
                        phi[i*prime[j]]=phi[i]*(prime[j]-1);
                }
        }
    for(int i=1;i<=N;++i)
        sumphi[i]=add(sumphi[i-1],mul(mul(i,i),phi[i]));
}
LoveLive s2(LoveLive x)
{
    x%=p;
    LoveLive res=x+1;
    res=mul(res,x);
    res=mul(res,2*x+1);
    res=mul(res,inv6);
    return res;
}
LoveLive s3(LoveLive x)
{
    x%=p;
    LoveLive tmp=mul(x,x+1);
    tmp=mul(tmp,inv2);
    return mul(tmp,tmp);
}
map<LoveLive,LoveLive> mp;
LoveLive calc(LoveLive x)
{
    if(x<=MAXN-10)
        return sumphi[x];
    if(mp[x])
        return mp[x];
    LoveLive res=s3(x);
    for(LoveLive l=2,r;l<=x;l=r+1)
        {
            r=x/(x/l);
            LoveLive tmp=add(s2(r),p-s2(l-1));
            tmp=mul(tmp,calc(x/l));
            res=add(res,p-tmp);
        }
    return mp[x]=res;
}
int main()
{
    p=read();
    n=read();
    init(MAXN-10);
    LoveLive ans=0;
    for(LoveLive l=1,r;l<=n;l=r+1)
        {
            r=n/(n/l);
            LoveLive tmp=add(calc(r),p-calc(l-1));
            tmp=mul(tmp,s3(n/l));
            ans=add(ans,tmp);
        }
    printf("%lld\n",ans);
    return 0;
}

杜教筛

标签：code   etc   直接   splay   prim   for   lse   getch   psi   

原文地址：https://www.cnblogs.com/jklover/p/10659279.html