杜教筛

莫比乌斯函数前缀和

51nod - 1244
令$S(n)=\sum_{i=1}^n\mu(i)$，求$S(n),1 \le n \le 10^{10}$

做法

$$[n=1]=\sum_{d|n}\mu(d)$$
$$\sum_{i=1}\sum_{d|i}\mu(d)=\sum_{i=1}\sum_{j=1}^{\lfloor \frac ni \rfloor}\mu(j)=\sum _{ i=1}^n S(\frac ni)=1$$
$$S(1)=1-\sum _{i=2}^nS(\frac ni)$$
线性筛出S(n)的前 $n^{\frac 23}$项，然后利用分块优化求解上式，总时间复杂度 $O(n^\frac 23)$

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long LL;
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int p[5000010],pcnt,mu[5000000+10],sum[5000000+10];
bool f[5000010],vis[5000010];
LL g[5000000+10],l,r;
void prepare(){
    int i,j;
    mu[1]=sum[1]=1;
    for(i=2;i<=5000000;i++){
        if(!f[i]){
            p[++pcnt]=i;
            mu[i]=-1;
        }
        for(j=1;i*p[j]<=5000000;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
}
void read(){
    Read(l),Read(r);
}
LL solve(LL n,LL now){
    if(n<=5000000)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=1;
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret-=(j-i+1)*solve(n/i,now*i);
    }
    return ret;
}
int main()
{
    prepare();
    read();
    LL t1=solve(r,1);
    memset(vis,0,sizeof vis);
    LL t2=solve(l-1,1);
    printf("%lld\n",t1-t2);
}

这种利用狄利克雷卷积筛前缀和的思路被称为杜教筛。

欧拉函数前缀和

51nod - 1239
令$S(n)=\sum_{i=1}^n\varphi(i)$，求$S(n),1 \le n \le 10^{10}$

做法

$$n=\sum_{d|n}\varphi(d)$$
$$\sum _{i=1}^n \sum_{d|i}\varphi(d)=\sum _{i=1}^n\sum_{j=1}^{\lfloor \frac ni \rfloor}\varphi(j)=\sum_{i=1}^n S(\frac ni)=\sum _{i=1}^ni$$
$$S(n)=\sum _{i=1}^ni-\sum_{i=2}S(\frac ni)$$
线性筛出S(n)的前 $n^{\frac 23}$项，然后利用分块优化求解上式，总时间复杂度 $O(n^\frac 23)$

代码

#include<cstdio>
#define MOD 1000000007
#include<algorithm>
using namespace std;
typedef long long LL;
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int p[5000010],pcnt,phi[5000010];
bool f[5000010],vis[5000010];
LL g[5000010],sum[5000010],n;
void prepare(){
    int i,j;
    sum[1]=phi[1]=1;
    for(i=2;i<=5000010;i++){
        if(!f[i]){
            p[++pcnt]=i;
            phi[i]=i-1;
        }
        for(j=1;i*p[j]<=5000010;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
        sum[i]=(sum[i-1]+phi[i])%MOD;
    }
}
LL solve(LL n,LL now){
    if(n<=5000010)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=0;
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        (ret-=(j-i+1)*solve(n/i,now*i)%MOD)%=MOD;
    }
    ret+=(n%MOD*((n+1)%MOD)%MOD*500000004%MOD)%MOD;
    return ret%MOD;
}
int main()
{
    prepare();
    Read(n);
    printf("%lld\n",(solve(n,1)+MOD)%MOD);
}

BZOJ3944 - Sum

上面两题二合一

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long LL;
template<class T>
void Read(T &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int p[1000010],pcnt,mu[1000000+10],sum[1000000+10],phi[1000000+10],T;
bool f[1000010],vis[1000010];
LL n,sum2[1000000+10];
typedef pair<LL,LL>pll;
pll g[1000000+10];
void prepare(){
    int i,j;
    mu[1]=sum[1]=1;
    phi[1]=sum2[1]=1;
    for(i=2;i<=1000000;i++){
        if(!f[i]){
            p[++pcnt]=i;
            mu[i]=-1;
            phi[i]=i-1;
        }
        for(j=1;i*p[j]<=1000000;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                phi[p[j]*i]=p[j]*phi[i];
                break;
            }
            mu[i*p[j]]=-mu[i];
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
        sum[i]=sum[i-1]+mu[i];
        sum2[i]=sum2[i-1]+phi[i];
    }
}
void read(){
    Read(n);
}
pll solve(LL n,LL now){
    if(n<=1000000)
        return pll(sum[n],sum2[n]);
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j;
    pll &ret=g[now]=pll(1,(n+1)*n/2);
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        pll t=solve(n/i,now*i);
        ret=pll(ret.first-(j-i+1)*t.first,ret.second-(j-i+1)*t.second);
    }
    return ret;
}
int main()
{
    prepare();
    Read(T);
    while(T--){
        memset(vis,0,sizeof vis);
        read();
        pll t=solve(n,1);
        printf("%lld %lld\n",t.second,t.first);
    }
}

HDU5608 - function

已知$N^2-3N+2=\sum_{d|N} f(d)$，求$\sum_{i=1}^{N} f(i)~mod~10^9+7$
多组数据，$T \leq 500, N \leq 10^9$，其中只有5组$N > 10^6$

做法

$$F(N)=N^2-3N+2=\sum_{d|N} f(d)$$
$$\begin{align}\sum _{i=1}^n\sum _{d|i}f(d)&=\sum _{ i=1}^ni^2-3i+2\\&=\sum _{i=1}^n (i-1)(i-2)\\&=2\times C(n,3)\\&=\sum_{i=1}^n\sum_{j=1}^{\lfloor \frac ni\rfloor}f(j)\end{align}$$
$$\sum _{i=1} ^n f(i)=2\times C(n,3)-\sum_{i=2}^n\sum_{j=1}^{\lfloor \frac ni\rfloor}f(j)$$
然而f(j)的前缀和并不好线性筛，所以时间复杂度为 $O(n^\frac 34)$
有没有其他方法呢？
我们莫比乌斯反演一下。
$$f(n)=\sum_{d|n} \mu(d)F(\frac nd)$$
$$\sum_{i=1}^nf(n)=\sum_{i=1}^n\sum_{d|i}\mu(d) F(\frac nd)=\sum_{i=1}^nF(i)\sum_{j=1}^{\lfloor \frac ni \rfloor}\mu(d)$$
我们就可以像求莫比乌斯函数的前缀和那样做，时间复杂度 $O(n^\frac 23)$

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#define MAXL 1000000
#define MOD 1000000007
using namespace std;
template<class T>
void Read(T &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int n,p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10],g[MAXL+10],ans,T;
bool f[MAXL+10];
bool vis[MAXL+10];
typedef long long LL;
void prepare(){
    int i,j;
    mu[1]=sum[1]=1;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            p[++pcnt]=i;
            mu[i]=-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
        sum[i]=mu[i]+sum[i-1];
    }
}
inline LL Sum(LL x){
    return x*(x-1)%MOD*(x-2)%MOD*((MOD+1)/3)%MOD;
}
int solve(int n,int now){
    if(n<=MAXL)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    int i,j,&ret=g[now]=1;
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret-1ll*(j-i+1)*solve(n/i,now*i))%MOD;
    }
    return ret;
}
int main()
{
    prepare();
    Read(T);
    while(T--){
        Read(n);
        if(n>=MAXL)
            memset(vis,0,sizeof vis);
        int i,j;
        ans=0;
        for(i=1;i<=n;i=j+1){
            j=n/(n/i);
            ans=(ans+((Sum(j)-Sum(i-1))%MOD*solve(n/i,i))%MOD)%MOD;
        }
        printf("%d\n",(ans+MOD)%MOD);
    }
}

51nod1237 - 最大公约数之和 V3

G=0;
for(i=1;i<N;i++)
for(j=1;j<=N;j++)
{
    G = (G + gcd(i,j)) % 1000000007;
}

求G

做法

令$A(n)=\sum _{i=1}^ngcd(i,n),F(n)=\sum _{i=1}^n A(i)$

$$A(n)=\sum _{d|n}d\times\varphi(\frac nd)$$
$$F(n)=\sum_{i=1}^n\sum_{d|i}d\times\varphi (\frac id)=\sum _{i=1}^ni\sum_{j=1}\varphi(j)$$
$$ans=F(n)\times 2-\sum_{i=1}^n i$$

计算欧拉函数的前缀和即可，时间复杂度$O(n^\frac 23)$

代码

#include<cstdio>
#define MOD 1000000007
#include<algorithm>
using namespace std;
typedef long long LL;
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int p[5000010],pcnt,phi[5000010];
bool f[5000010],vis[5000010];
LL g[5000010],sum[5000010],n,ans;
void prepare(){
    int i,j;
    sum[1]=phi[1]=1;
    for(i=2;i<=5000010;i++){
        if(!f[i]){
            p[++pcnt]=i;
            phi[i]=i-1;
        }
        for(j=1;i*p[j]<=5000010;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
        sum[i]=(sum[i-1]+phi[i])%MOD;
    }
}
LL solve(LL n,LL now){
    if(n<=5000010)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=0;
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        (ret-=(j-i+1)*solve(n/i,now*i)%MOD)%=MOD;
    }
    ret+=(n%MOD*((n+1)%MOD)%MOD*500000004%MOD)%MOD;
    return ret%MOD;
}
inline LL Sum(LL x){
    return x%MOD*((x+1)%MOD)%MOD*500000004%MOD;
}
int main()
{
    prepare();
    Read(n);
    LL i,j;
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ans=(ans+(Sum(j)-Sum(i-1))*solve(n/i,i))%MOD;
    }
    printf("%lld\n",((((ans+MOD)%MOD)*2-Sum(n))%MOD+MOD)%MOD);
}

51nod1227 - 平均最小公倍数

Lcm(a,b)表示a和b的最小公倍数，A(n)表示Lcm(n,i)的平均数(1 <= i <= n)，
例如：A(4) = (Lcm(1,4) + Lcm(2,4) + Lcm(3,4) + Lcm(4,4)) / 4 = (4 + 4 + 12 + 4) / 4 = 6。
F(a, b) = A(a) + A(a + 1) + …… A(b)。(F(a,b) = ∑A(k), a <= k <= b)
例如：F(2, 4) = A(2) + A(3) + A(4) = 2 + 4 + 6 = 12。
给出a,b，计算F(a, b)，由于结果可能很大，输出F(a, b) % 1000000007的结果即可。

做法

令$A(n)=\frac{\sum_{i=1}^nlcm(i,n)}n,F(n)=\sum_{i=1}^n A(n)$

$$\begin{align}A(n)&=\frac{n\sum_{d|n}\sum_{i=1}^{\frac nd}[gcd(i,\frac nd)=1]i}n\\&=\sum_{d|n}\sum_{i=1}^{\frac nd}[gcd(i,\frac nd)=1]i\\&=\sum_{d|n}\frac nd\times \frac{\varphi(\frac nd)+[\frac nd=1]}2\\&=1+\sum_{d|n}\frac{d\times \varphi (d)}2\end{align}$$
令 $G(n)=2\times F(n)-n$
$$G(n)=\sum_{i=1}^n\sum_{d|i}d\times\varphi(d)=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac ni\rfloor}j\times \varphi(j)$$
怎样快速求 $\sum_{i=1}^{n}i\times\varphi(i)$呢？
我们构造一下，发现 $\sum_{d|n}d\times \varphi(d)\times \frac nd=n^2$
$$\sum_{i=1}^n\sum_{d|i}d\times \varphi(d)\times \frac id=\sum_{i=1}^ni\sum_{j=1}^{\lfloor\frac ni\rfloor}j\times \varphi(j)$$
线性筛 $\sum_{i=1}^{n}i\times\varphi(i)$的qian $n^\frac23$项，总时间复杂度 $O(n^\frac 23)$

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
#define MOD 1000000007
#define MAXL 1000000
typedef long long LL;
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
LL l,r;
int p[MAXL+10],sum[MAXL+10],pcnt,phi[MAXL+10],ans;
LL g[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
void prepare(){
    int i,j;
    phi[1]=sum[1]=1;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            p[++pcnt]=i;
            phi[i]=i-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
        sum[i]=(sum[i-1]+1ll*i*phi[i])%MOD;
    }
}
LL Sum2(LL n){
    n%=MOD;
    return n%MOD*((n+1)%MOD)%MOD*(2*n+1)%MOD*((MOD+1)/6)%MOD;
}
LL Sum(LL n){
    n%=MOD;
    return n*(n+1)%MOD*((MOD+1)/2)%MOD;
}
LL Sum3(LL n){
    return 1ll*Sum(n)*Sum(n)%MOD;
}
LL solve(LL n,LL now){
    if(n<=MAXL)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=Sum2(n);
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret-(Sum(j)-Sum(i-1))*solve(n/i,now*i))%MOD;
    }
    return ret;
}
LL cal(LL n){
    memset(vis,0,sizeof vis);
    LL i,j,ret(0);
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret+(j-i+1)*solve(n/i,i))%MOD;
    }
    return (ret+n)*((MOD+1)>>1)%MOD;
}
int main()
{
    prepare();
    Read(l),Read(r);
    printf("%lld\n",((cal(r)-cal(l-1))%MOD+MOD)%MOD);
}

51nod1238 - 最小公倍数之和 V3

G=0;
for(i=1;i<N;i++)
for(j=1;j<=N;j++)
{
    G = (G + lcm(i,j)) % 1000000007;
}

求G

做法

令$A(n)=\sum _{i=1}^nlcm(i,n),F(n)=\sum _{i=1}^n A(i)$

$$\begin{align}A(n)&=n\sum_{d|n}\sum_{i=1}^{\frac nd}[gcd(i,\frac nd)=1]i\\&=n\sum_{d|n}\frac nd\times \frac{\varphi(\frac nd)+[\frac nd=1]}{2}\\&=n+n\sum_{d|n}\frac{d\times\varphi(d)}{2}\end{align}$$
令 $G(n)=2\times F(n)-\sum _{i=1}^ni$
$$G(n)=\sum_{i=1}^ni\sum_{d|i}d\times \varphi(d)=\sum_{i=1}^ni\sum_{j=1}^{\lfloor\frac ni\rfloor}j^2\times\varphi(j)$$
怎样快速求 $\sum_{i=1}^{n}i^2\times\varphi(i)$呢？
类似上一题，我们构造一下
$$\sum_{d|n}d^2\times\varphi(d)\times (\frac nd) ^2=n^3$$
$$\sum _{i=1}^n\sum_{d|i}d^2\times\varphi(d)\times (\frac id) ^2=\sum_{i=1}^ni^2\sum_{j=1}^{\lfloor \frac ni \rfloor}j^2\times \varphi(j)=\sum_{i=1}^n i^3$$
$$\sum_{i=1}^{n}i^2\times\varphi(i)=\sum _{i=1}^n i^3-\sum_{i=2}^ni^2\sum_{j=1}^{\lfloor \frac ni \rfloor}j^2\times \varphi(j)$$
线性筛 $\sum_{i=1}^{n}i^2\times\varphi(i)$的前 $n^\frac 23$项，总的时间复杂度 $O(n^\frac 23)$，答案就是 $G(n)$

代码

#include<cstdio>
#include<algorithm>
using namespace std;
#define MOD 1000000007
#define MAXL 10000000
typedef long long LL;
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
LL n;
int p[MAXL+10],sum[MAXL+10],pcnt,phi[MAXL+10],ans;
LL g[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
void prepare(){
    int i,j;
    phi[1]=sum[1]=1;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            p[++pcnt]=i;
            phi[i]=i-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
        sum[i]=(sum[i-1]+1ll*i*i%MOD*phi[i])%MOD;
    }
}
LL Sum2(LL n){
    n%=MOD;
    return n%MOD*((n+1)%MOD)%MOD*(2*n+1)%MOD*((MOD+1)/6)%MOD;
}
LL Sum(LL n){
    n%=MOD;
    return n*(n+1)%MOD*((MOD+1)/2)%MOD;
}
LL Sum3(LL n){
    return 1ll*Sum(n)*Sum(n)%MOD;
}
LL solve(LL n,LL now){
    if(n<=MAXL)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=Sum3(n);
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret-(Sum2(j)-Sum2(i-1))*solve(n/i,now*i))%MOD;
    }
    return ret;
}
int main()
{
    prepare();
    Read(n);
    LL i,j;
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ans=(ans+(Sum(j)-Sum(i-1))*solve(n/i,i))%MOD;
    }
    printf("%d\n",(ans+MOD)%MOD);
}

SPOJ DIVCNT2 - Counting Divisors(square)

$$S_2(n) = \sum _{i=1}^n \sigma_0(i^2).$$
求 $S_2(N)$.

做法

$$\sigma_0(n^2)=\sum_{d|n}2^{\omega(d)}$$
$\omega(n)$表示n的不同质因子的数量，可以看做在 $d^2$中减少 $d$的一个质因数集合。
$2^{\omega(n)}$可以看做n的无平方因子的数的个数。
$$2^{\omega(n)}=\sum_{d|n}\mu^2(d)$$
$$\begin{align}S_2(n)&amp;=\sum_{i=1}^n\sum_{j|i}\sum_{d|j}\mu^2(d)\\&=1*1*\mu^2\\&=\mu^2*\sigma\\&=\sum_{i=1}^n\sum_{d|i}\mu^2(\frac id)\sigma(d)\\&=\sum_{i=1}^n\mu^2(i)\sum_{j=1}^{\lfloor\frac ni\rfloor}\sigma(d)\end{align}$$
时间复杂度据说 $O(n^\frac 23)$

代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
template<class T>
void Read(T &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
typedef long long LL;
int p[10000010],mu[100000010],cf[100000010],L=1e7,n,pcnt,T,smu[100000010];
bool f[100000010];
char cs[100000010];
LL a[10000+10],mxa;
void read(){
    Read(n);
}
void prepare(int n){
    int i,j;
    smu[1]=mu[1]=1;
    cf[1]=1;
    for(i=2;i<=n;i++){
        if(!f[i]){
            mu[i]=-1;
            p[++pcnt]=i;
            cs[i]=1,cf[i]=2;
        }
        for(j=1;i*p[j]<=n;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                cs[i*p[j]]=cs[i]+1;
                cf[i*p[j]]=cf[i]/(cs[i]+1)*(cs[i]+2);
                break;
            }
            mu[i*p[j]]=-mu[i];
            cs[i*p[j]]=1;
            cf[i*p[j]]=cf[i]<<1;
        }
        cf[i]+=cf[i-1];
        smu[i]=abs(mu[i])+smu[i-1];
    }
}
LL cal(LL x){
    if(x<=L)
        return cf[x];
    LL ret=0,i,j;
    for(i=1;i<=x;i=j+1){
        j=x/(x/i);
        ret+=(j-i+1)*(x/i);
    }
    return ret;
}
LL Sum(LL x){
    if(x<=L)
        return smu[x];
    LL ret(0);
    for(LL i=1;i*i<=x;i++)
        if(mu[i])
            ret+=mu[i]*(x/(i*i));
    return ret;
}
LL solve(LL n){
    int t=sqrt(n+0.5);
    LL i,j;
    LL ret(0),pre=smu[t];
    for(i=1;i<=t;i++)
        if(mu[i])
            ret+=cal(n/i);
    for(i=t+1;i<=n;i=j+1){
        j=n/(n/i);
        LL tt=Sum(j);
        ret+=(tt-pre)*cal(n/i);
        pre=tt;
    }
    return ret;
}
int main()
{
    Read(T);
    for(int i=1;i<=T;i++)
        Read(a[i]),mxa=max(mxa,a[i]);
    L=max(1.0*L,pow(mxa,2.0/3));
    prepare(L);
    for(int i=1;i<=T;i++)
        printf("%lld\n",solve(a[i]));
}

51nod1222 - 最小公倍数计数

定义F(n)表示最小公倍数为n的二元组的数量。
即：如果存在两个数（二元组）X，Y(X <= Y)，它们的最小公倍数为N，则F(n)的计数加1。
例如：F(6) = 5，因为[2,3] [1,6] [2,6] [3,6] [6,6]的最小公倍数等于6。

给出一个区间[a,b]，求最小公倍数在这个区间的不同二元组的数量。
例如：a = 4，b = 6。符合条件的二元组包括：
[1,4] [2,4] [4,4] [1,5] [5,5] [2,3] [1,6] [2,6] [3,6] [6,6]，共10组不同的组合。
Input
输入数据包括2个数：a, b，中间用空格分隔（1 <= a <= b <= 10^11)。
Output
输出最小公倍数在这个区间的不同二元组的数量。

做法

令$S(n)=\sum_{i=1}^n F(i)$，$ans[a,b]=S(b)-S(a-1)$

$$S(n)=\sum_{d=1}^n\sum_{i=1}^{n}\sum_{j=1}^n[gcd(a,b)=1][a\le b][abd<=n]$$
莫比乌斯反演一下, $$S(n)=\sum_{r=1}^{\sqrt n}\mu (r)\sum_{d=1}^n\sum_{i=1}^{n}\sum_{j=1}^n[a\le b][abd<=\frac n{r^2}]$$
我们发现，其实 $\sum_{r=1}^{\sqrt n}\mu (r)\sum_{d=1}^n\sum_{i=1}^{n}\sum_{j=1}^n[abd<=\frac n{r^2}]$还是很好算的，我们枚d,i就可以了，然后分别统计一个相同，两个相同，三个相同的数量，然后就可以轻易算出答案了，据说复杂度 $O(n^\frac 23)$

代码

#include<cstdio>
#include<algorithm>
#include<cmath>
#define MAXL 1000000
using namespace std;
typedef long long LL;
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
LL l,r;
int p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10];
bool f[MAXL+10];
void read(){
    Read(l),Read(r);
}
void prepare(){
    mu[1]=sum[1]=1;
    int i,j;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            p[++pcnt]=i;
            mu[i]=-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
}
LL cal(LL n){
    LL i,cnt1=0,cnt2=0,cnt3=0,j,k;
    for(i=1;i*i*i<=n;i++){
        cnt3++;
        for(j=i+1;i*j*j<=n;j++){
            cnt2++;
        //  for(k=j+1;i*j*k<=n;k++)
        //      cnt1++;
            cnt1+=max(0ll,n/(i*j)-j);
        }
    }
    for(i=1;i*i*i<=n;i++)
        cnt2+=max(n/(i*i)-i,0ll);
    return cnt1*6+cnt2*3+cnt3;
}
LL solve(LL n){
    LL i,j,ans=0;
    for(i=1;i*i<=n;i=j+1){
        j=n/(n/(i*i));
        j=sqrt(j+0.5);
        ans+=(sum[j]-sum[i-1])*cal(n/(i*i));
    }
    return (ans+n)/2;
}
int main()
{
    prepare();
    read();
    printf("%lld\n",solve(r)-solve(l-1));
}

BZOJ4176 - Lucas的数论

去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。

在整理以前的试题时，发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”，其中 表示i的约数个数。他现在长大了，题目也变难了。
求如下表达式的值：

其中 表示ij的约数个数。
他发现答案有点大，只需要输出模1000000007的值。

做法

$d|ij$等价于$\frac d{gcd(i,d)}|j$

$$\begin{align}\sum_{i=1}^n\sum_{d=1}^{n^2}\lfloor\frac{n\times gcd(i,d)}{d}\rfloor&=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd \rfloor}\sum_{j=1}^{\lfloor\frac{ n^2}{d}\rfloor}\lfloor\frac nj\rfloor[gcd(i,j)=1]\\&=\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac nd\rfloor}\sum_{j=1}^n\lfloor \frac nj\rfloor\sum_{k|i,k|j}\mu(k)\\&=\sum_{k=1}^n\mu(k)\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{j=1}^{\lfloor\frac nk \rfloor}\lfloor\frac n{kj}\rfloor\\&=\sum_{k=1}^n\mu(k)\sum_{d=1}^n\lfloor\frac{n}{dk}\rfloor\sum_{j=1}^{\lfloor\frac nk \rfloor}\lfloor\frac n{kj}\rfloor\\&=\sum_{k=1}^n\mu(k)\sum_{d=1}^{\lfloor\frac nk \rfloor}\lfloor\frac{n}{dk}\rfloor\sum_{j=1}^{\lfloor\frac nk \rfloor}\lfloor\frac n{kj}\rfloor\\&=\sum_{k=1}^n\mu(k)(\sum_{d=1}^{\lfloor\frac nk \rfloor}\lfloor\frac{n}{dk}\rfloor)^2\end{align}$$

据说时间复杂度$O(n^\frac 34)$

代码

#include<cstdio>
#include<algorithm>
#define MOD 1000000007
#define MAXL 1000000
using namespace std;
typedef long long LL;
LL n;
int p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
LL g[MAXL+10],ans;
void prepare(){
    mu[1]=sum[1]=1;
    int i,j;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            p[++pcnt]=i;
            mu[i]=-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
}
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
    }
}
inline LL cal(LL n){
    LL i,j,ret(0);
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret+(j-i+1)*(n/i))%MOD;
    }
    return ret;
}
LL solve(LL n,LL now){
    if(n<=MAXL)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=1;
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret-(j-i+1)*solve(n/i,now*i))%MOD;
    }
    return ret;
}
int main()
{
    prepare();
    Read(n);
    LL i,j;
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        LL t=cal(n/i);
        ans=(ans+(solve(j,n/i)-solve(i-1,n/((i-1)?(i-1):1)))*t%MOD*t)%MOD;
    }
    printf("%lld\n",(ans+MOD)%MOD);
}

51nod1220 - 约数之和

d(k)表示k的所有约数的和。d(6) = 1 + 2 + 3 + 6 = 10。
定义S(N) = ∑1<=i<=N ∑1<=j<=N d(i*j)。
例如：S(3) = d(1) + d(2) + d(3) + d(2) + d(4) + d(6) + d(3) + d(6) + d(9) = 59，S(1000) = 563576517282。
给出正整数N，求S(N)，由于结果可能会很大，输出Mod 1000000007(10^9 + 7)的结果。

做法

和上一题类似，做法还是看http://www.51nod.com/question/index.html#!questionId=1027

代码

#include<cstdio>
#include<algorithm>
#define MOD 1000000007
#define MAXL 1000000
using namespace std;
typedef long long LL;
LL n;
int p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
LL g[MAXL+10],ans;
void prepare(){
    mu[1]=sum[1]=1;
    int i,j;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            p[++pcnt]=i;
            mu[i]=-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
        sum[i]=(sum[i-1]+mu[i]*i)%MOD;
    }
}
void Read(LL &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
    }
}
inline LL Sum(LL n){
    return n*(n+1)/2%MOD;
}
inline LL cal(LL n){
    LL i,j,ret(0);
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret+(Sum(j)-Sum(i-1))*(n/i))%MOD;
    }
    return ret;
}
inline LL cal2(LL n){
    LL i,j,ret(0);
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret+(j-i+1)*Sum(n/i))%MOD;
    }
    return ret;
}
LL solve(LL n,LL now){
    if(n<=MAXL)
        return sum[n];
    if(vis[now])
        return g[now];
    vis[now]=1;
    LL i,j,&ret=g[now]=1;
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret-(Sum(j)-Sum(i-1))*solve(n/i,now*i))%MOD;
    }
    return ret;
}
int main()
{
    prepare();
    Read(n);
    LL i,j;
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        ans=(ans+(solve(j,n/i)-solve(i-1,n/((i-1)?(i-1):1)))*cal2(n/i)%MOD*cal(n/i))%MOD;
    }
    printf("%lld\n",(ans+MOD)%MOD);
}

BZOJ3512 - DZY Loves Math IV

给定n,m，求 模10^9+7的值。

做法

令$S(n,m)=\sum_{i=1}^m\varphi(i\times n),ans=\sum_{i=1}^nS(i,m)$

$$S(n,m)=\begin{cases}p^{a-1}S(\frac np,m) &存在一个n的质因数p的质数>1\\ \varphi(p)*sum(\frac{n}{p},m)+sum(n,\lfloor \frac{m}{p}\rfloor) \end{cases}$$
然后记忆化搜索+杜教筛求解即可。

代码

#include<cstdio>
#include<algorithm>
#include<map>
#define MAXL 1000000
#define MOD 1000000007
using namespace std;
using namespace __gnu_cxx;
typedef pair<int,int>pii;
map<pii,int>h;
void Read(int &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
#define MAXN 100000
typedef long long LL;
int n,m,ans[MAXN+10],sf[MAXL+10],p[MAXL+10],pcnt,sum[MAXL+10],phi[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
LL g[MAXL+10];
void read(){
    Read(n),Read(m);
}
void prepare(){
    int i,j;
    phi[1]=sum[1]=1;
    for(i=2;i<=MAXL;i++){
        if(!f[i]){
            sf[i]=i;
            p[++pcnt]=i;
            phi[i]=i-1;
        }
        for(j=1;i*p[j]<=MAXL;j++){
            f[i*p[j]]=1;
            sf[i*p[j]]=sf[i];
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
        sum[i]=(sum[i-1]+phi[i])%MOD;
    }
}
LL Get_fac(int n){
    int i,ret(1);
    for(i=1;p[i]*p[i]<=n;i++){
        if(n%p[i]==0){
            n/=p[i];
            while(n%p[i]==0)
                ret*=p[i],n/=p[i];
        }
    }
    return ret;
}
LL Sum(LL n){
    return n*(n+1)/2%MOD;
}
LL cal(LL n){
    if(n<=MAXL)
        return sum[n];
    if(vis[m/n])
        return g[m/n];
    vis[m/n]=1;
    LL i,j,&ret=g[m/n]=Sum(n);
    for(i=2;i<=n;i=j+1){
        j=n/(n/i);
        ret=(ret-(j-i+1)*cal(n/i))%MOD;
    }
    return ret;
}
LL solve(LL n,LL m){
    if(h.count(pii(n,m)))
        return h[pii(n,m)];
    if(!m)
        return 0;
    if(n==1)
        return cal(m);
    if(m==1)
        return phi[n];
    return h[pii(n,m)]=(phi[sf[n]]*solve(n/sf[n],m)+solve(n,m/sf[n]))%MOD;
}
LL solve(){
    int i,t,ret=0;
    for(i=1;i<=n;i++){
    //  printf("%d\n",i);
        t=Get_fac(i);
        if(t!=1)
            ret=(ret+1ll*t*ans[i/t])%MOD;
        else
            ret=(ret+(ans[i]=solve(i,m)))%MOD;
    }
    return ret;
}
int main()
{
    prepare();
    read();
    printf("%lld\n",(solve()+MOD)%MOD);
}