杜教筛

前置知识

积性函数

定义：满足对于任意互质的 \(x,y\)，有：\(f(x\cdot y)=f(x)\cdot f(y)\)。则称 \(f\) 为积性函数，例如 \(\mu,\varphi\)。

完全积性函数，则定义中没有 互质 这一前提，比如：

• \(\epsilon(n)=[n=1]\)

• \(I(n)=1\)

• \(id(n)=n\)

---

欧拉函数

\(\to\)

---

莫比乌斯函数

对于莫比乌斯函数，有：

\[\sum_{d|n}\mu(d)=\epsilon(n) \]

证明：

由定义可知，将 \(n\) 进行质因数拆分成 \(\prod_{i=1}^s {p_i}^{a_i}\) 的形式，所有对该式子有贡献的项 \(d\)，必定满足：\(d=\prod_{i=1}^s {p_i}^{k},k=0\ or \ 1\)。

那么就是每个质因子取或不取，答案就是 取偶数个 的方案数减去 取奇数个的方案数。

很显然当 \(s\) 为奇数时，奇偶的方案可以达成一一对应的关系，即每一个取奇数个的方案，剩下的没取的就是一个取偶数个的方案。

当 \(s\) 为偶数的时候，其贡献为：

\[\sum_{i=0}^n\begin{pmatrix} n\\i \end{pmatrix}\cdot (-1)^i \]

也就是：

\[\sum_{i=0}^n\begin{pmatrix} n\\i \end{pmatrix}\cdot (-1)^{i} \]

当 \(n\ge 1\) 时，根据二项式定理，转化为：

\[(1-1)^n \]

那么其值为 \(0\)。

至此，证毕。

---

狄利克雷卷积

对于两个积性函数 \(f,g\)，其卷积：

\[[x^n](f \ast g)=\sum_{d|n} [x^d]f(x)\cdot [x^{\frac{n}{d}}]g(x) \]

（这里因为最近多项式式子写多了，写成了这个鬼东西，但我相信并不影响理解）

狄利克雷卷积满足交换，结合律。

易得一些特殊的关系：

• \(f \ast \epsilon=f\)

• \(\mu \ast I= \epsilon\)

• \(\varphi \ast I= id\)

• \(\varphi=id\ast \mu\)

---

莫比乌斯反演

若有：

\[g(n)=\sum_{d|n} f(d) \]

将其写成：

\[g=f\ast I \]

由前面可得：

\[g\ast \mu=f \]

展开得：

\[f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d}) \]

（原来莫反是这个东西，我一直以来的莫反题都是怎么做的？可能是暴力拆式子然后合并，反正结果是一样的）

---

杜教筛

那么如何非线性求积性函数的前缀和呢？

比如对于积性函数 \(f\)，求 \(\sum_{i=1}^n\limits f(i)\)。

我们随便拿一个积性函数和它卷一下，得到：

\[\sum_{i=1}^n \sum_{d|i}g(d)\times f(\frac{i}{d}) \]

\[\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} f(i) \]

当 \(d=1\) 时，右边的 \(\sum\) 的结果为卷之前的结果。

那么就得到了：

\[\sum_{i=1}^n f(i)=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} f(i)-\sum_{d=2}^n g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} f(i) \]

得到这么复杂的式子有用吗？它的复杂度降下去了吗？

事实上已经降下去了，因为对于左边那个式子，\(\mu\) 和 \(\varphi\) 都可以通过卷 \(I\) 转化成可以 \(\operatorname O(1)\) 求解的东西。右边的话，看着就是一个非常可以递归的东西。

比方说对于 \(f=\mu\)，我们整理一下，得到：

\[\sum_{i=1}^n \mu(i)=1-\sum_{d=2}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} \mu(i) \]

同理，\(\varphi\) 可以得到：

\[\sum_{i=1}^n \varphi(i)=\begin{pmatrix} n\\ 2\end{pmatrix}-\sum_{d=2}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} \varphi(i) \]

那么对于当前的 \(n\)，整除分块一下，然后递归调用。

在整除分块的时候，证明过向下取整和除运算先后顺序改一下不会影响结果。

所以其实整个过程中调用的 \(n'\) 其实就是 \(n\) 整除得到的值之一，只有 \(\sqrt n\) 种，然后递归是 \(\log\) 的，记忆化以后每一层都是 \(\sqrt {n'}\) 的复杂度。复杂度我不会严谨的计算，但是我感性理解它是对的。

然后对于比较小的 \(n'\)，递归调用比较愚蠢，可以在之前线性筛一波。

（会不会快我就不知道了）

可以证明其复杂度为 \(\operatorname O(n^{\frac{2}{3}})\)。

“【证明显然，由读者自行完成】”

---

代码

模板题的 \(n\) 是 \(\texttt{maxint}\) 范围内的，那么就会出现一些比较恶心的东西。

比如说过程量炸 long long。

#include <stdio.h>
#include <map>
#define LL long long
#define ULL long long
using namespace std;
const int N=1e7+3;
const int P=8e5+3;
inline int rin()
{
    int s=0;
    bool bj=false;
    char c=getchar();
    for(;(c>'9'||c<'0')&&c!='-';c=getchar());
    if(c=='-')bj=true,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())s=(s<<1)+(s<<3)+(c^'0');
    if(bj)s=-s;
    return s;
}

LL mu[N];
LL var[N];
bool pri[N];
int prime[P];
int cutt;

inline void init()
{
    mu[1]=var[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!pri[i])
        {
            prime[++cutt]=i;
            var[i]=i-1;
            mu[i]=-1;
        }
        for(int j=1;j<=cutt;j++)
        {
            int now=i*prime[j];
            if(now>=N)break;
            pri[now]=true;
            if(i%prime[j]==0)
            {
                var[now]=var[i]*prime[j];
                break;
            }
            var[now]=var[i]*(prime[j]-1);
            mu[now]=-mu[i];
        }
        var[i]+=var[i-1];
        mu[i]+=mu[i-1];
    }
    return;
}

map<int,LL>q_mu;
map<int,LL>q_var;
inline LL work_mu(int n)
{
    if(n<N)return mu[n];
    LL sum=q_mu[n];
    if(sum)return sum;
    
    sum=1;
    for(int l=2,r;l<=n&&r!=2147483647;l=r+1)
    {
        int now=(n/l);
        r=n/now;
        sum-=work_mu(now)*(r-l+1);
    }
    q_mu[n]=sum;
    return sum;
}
inline LL work_var(int n)
{
    if(n<N)return var[n];
    LL sum=q_var[n];
    if(sum)return sum;
    
    ULL ans;
    ans=(n&1)?(1ULL*(1LL+n)>>1)*n:(1ULL*n>>1)*(1LL+n);
    for(int l=2,r;l<=n&&r!=2147483647;l=r+1)
    {
        int now=(n/l);
        r=n/now;
        ans-=1ULL*work_var(now)*(r-l+1);
    }
    sum=ans;
    q_var[n]=sum;
    return sum;
}

inline void work()
{
    int n=rin();
    printf("%lld ",work_var(n));
    printf("%lld\n",work_mu(n));
}
int main()
{
    int i,j;
    init();
    for(int T=rin();T;T--)work();
    return 0;
}

练习题

简单的数学题

\[\sum_{d=1}^{n}\sum_{i=1}^n\sum_{j=1}^ni\times j\times d\times [\gcd(i,j)=d] \]

\[\sum_{d=1}^{n}d^3\times \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}i\times j\times [\gcd(i,j)=1] \]

\[\sum_{d=1}^{n}d^3\times \sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\times p^2\times \sum_{i=1}^{\lfloor\frac{n}{T}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{T}\rfloor}i\times j \]

\[\sum_{d=1}^{n}d^3\times \sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\times p^2\times \begin{pmatrix}\lfloor\frac{n}{T}\rfloor+1\\2\end{pmatrix}^2 \]

\[\sum_{T=1}^{n}\sum_{p|T}\mu(p)\times p^2 \times (\frac{T}{p})^3\times \begin{pmatrix}\lfloor\frac{n}{T}\rfloor+1\\2\end{pmatrix}^2 \]

\[\sum_{T=1}^{n}(T)^2 \cdot \begin{pmatrix}\lfloor\frac{n}{T}\rfloor+1\\2\end{pmatrix}^2\cdot \sum_{p|T}\mu(p)\times \frac{(T)}{p} \]

\[\sum_{T=1}^{n}\varphi(T)\times (T)^2 \times \begin{pmatrix}\lfloor\frac{n}{T}\rfloor+1\\2\end{pmatrix}^2 \]

整除分块一下，然后杜教筛即可，将 \(\varphi\) 和 \(id^2\) 卷一下。

#include <stdio.h>
#include <unordered_map>
#define LL long long
using namespace std;
const int N=1e7+3;
const int P=7e5+3;
inline LL rin()
{
    LL s=0;
    bool bj=false;
    char c=getchar();
    for(;(c>'9'||c<'0')&&c!='-';c=getchar());
    if(c=='-')bj=true,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())s=(s<<1)+(s<<3)+(c^'0');
    if(bj)s=-s;
    return s;
}
int p;
LL n;

inline int prpr(int x,int y){return 1LL*x*y%p;}
inline LL prpr(LL x,LL y){return (x%p)*(y%p)%p;}
inline LL p2(LL x){x%=p;return x*x%p;}
inline LL n2(LL now)
{
    LL a=now,b=now+1,c=(now<<1)+1;
    for(;true;)
    {
        if(a%2==0){a>>=1;break;}
        if(b%2==0){b>>=1;break;}
        if(c%2==0){c>>=1;break;}
    }
    for(;true;)
    {
        if(a%3==0){a/=3;break;}
        if(b%3==0){b/=3;break;}
        if(c%3==0){c/=3;break;}
    }
    return prpr(prpr(a%p,b%p),c%p);
}
inline LL n3(LL now){return (now&1)?(prpr(p2(now+1>>1),p2(now))):(prpr(p2(now>>1),p2(now+1)));}

LL var[N];
bool pri[N];
int prime[P];
int cutt;

inline void init()
{
    var[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!pri[i])prime[++cutt]=i,var[i]=i-1;
        for(int j=1;j<=cutt;j++)
        {
            int now=i*prime[j];
            if(now>=N)break;
            pri[now]=true;
            if(i%prime[j]==0)
            {
                var[now]=var[i]*prime[j];
                break;
            }
            var[now]=var[i]*(prime[j]-1);
        }
        var[i]=var[i]*i%p*i%p;
        var[i]=(var[i-1]+var[i])%p;
    }
    return;
}

unordered_map<LL,LL>q;
inline LL work_var(LL n)
{
    if(n<N)return var[n];
    if(q.count(n))return q[n];
    LL sum=n3(n);
    for(LL l=2,r;l<=n;l=r+1)
    {
        LL now=n/l;
        r=n/now;
        sum-=prpr(n2(r)-n2(l-1),work_var(now));
        sum%=p;
    }
    sum=(sum%p+p)%p;
    return q[n]=sum;
}

inline void work(LL n)
{
    LL ans=0;
    for(LL l=1,r;l<=n;l=r+1)
    {
        LL now=n/l;
        r=n/now;
        ans+=prpr(p2((now&1)?(prpr(now+1>>1,now)):(prpr(now+1,now>>1))),work_var(r)-work_var(l-1));
    }
    ans=(ans%p+p)%p;
    printf("%lld\n",ans);
    return;
}
int main()
{
    int i,j;
    p=rin();n=rin();
    init();
    work(n);
    return 0;
}