杜教筛学习笔记

刀尖舞蹈

如果数据很容易爆 l o n g   l o n g long\ long long long的话,你有几种处理方式:

1:在爆炸的边缘试探

正确写法: a n s = ( U L L ) n ∗ ( n + 1 L L ) / 2 L L ; ans = (ULL) n * (n + 1LL) / 2LL; ans=(ULL)n(n+1LL)/2LL;

错误写法: a n s = ( U L L ) n ∗ ( n + 1 ) / 2 ; ans = (ULL) n * (n + 1) / 2; ans=(ULL)n(n+1)/2;

1 L L 1LL 1LL表示加一个 l o n g   l o n g long\ long long long类型的 1 1 1,尽量不要跨类型加。

2: i n t _ 128 int\_128 int_128禁术

3:利用语言特性(大数模乘)

r e t u r n ( L L ) a ∗ b − ( L L ) ( ( L D B ) a ∗ b / p ) ∗ p ) return (LL) a * b - (LL)((LDB)a * b / p) * p) return(LL)ab(LL)((LDB)ab/p)p)

4:高精度(最不幸的)

5: 改python

杜教筛

用于在快于线性的时间内求一个函数的前缀和。

你的目标是求 f f f的前缀和 F ( n ) F(n) F(n)

我们寻找一个好求前缀和的函数 g g g,(一般是 1 1 1函数)

构造 h = f ∗ g h=f*g h=fg,这里 ∗ * 是狄利克雷卷积。记 H H H h h h的前缀和。

那么有:

F ( n ) = H ( n ) − ∑ i = 2 n g ( i ) F ( n i ) g ( 1 ) F(n)={\frac{H(n)-\sum_{i=2}^ng(i)F(\frac{n}{i})}{g(1)}} F(n)=g(1)H(n)i=2ng(i)F(in)

我们对 n i \frac{n}{i} in数论分块来做,这就要求 H , g H,g H,g是可以快速求出前缀和的。

首先预处理一部分可以线性得到的 F F F,一般预处理 n 2 3 n^{\frac{2}{3}} n32为宜。

对于每个询问,跑杜教筛。可以用 c + + 11 c++11 c++11特性下的 u n o r d e r e d _ m a p unordered\_map unordered_map来存储求过的 F F F

ϕ \phi ϕ的前缀和。

ϕ ∗ 1 = I = i d 1 \phi*1=I=id_1 ϕ1=I=id1

预处理就是一般的 ϕ \phi ϕ筛。

inline LL dlsphi(int n) {
    if (n < N) return F[n];
    if (q1[n]) return q1[n];
    LL res = 0;
    for (register int L = 2, R; R < INF && L <= n; L = R + 1) {
        R = n / (n / L);
        res += (ULL) (R - L + 1) * dlsphi(n / L);
    }
    return q1[n] = (ULL)n * (n + 1LL) / 2LL - res;
}

这里的 n n n i n t int int级别的,所以要特殊处理下防止越界。

μ \mu μ的前缀和

μ ∗ 1 = ϵ \mu*1=\epsilon μ1=ϵ

预处理也就是一般的 μ \mu μ筛。

inline int dlsmu(int n) {
    if (n < N) return F_[n];
    if (q2[n]) return q2[n];
    int res = (n >= 1);
    for (register int L = 2, R; R < INF && L <= n; L = R + 1) {
        R = n / (n / L);
        res -= (ULL) (R - L + 1) * dlsmu(n / L);
    }
    return q2[n] = res;
}
F ( n ) = ∑ i = 1 n ϕ ( i ) ∗ i 2 F(n)=\sum_{i=1}^n\phi(i)*i^2 F(n)=i=1nϕ(i)i2

f ( i ) = ϕ ( i ) ∗ i 2 f(i)=\phi(i)*i^2 f(i)=ϕ(i)i2
看到了 i 2 i^2 i2,我们考虑让构造的 g g g里面也有 i 2 i^2 i2达到消掉的目的。

g ( i ) = i 2 g(i)=i^2 g(i)=i2 h = f ∗ g h=f*g h=fg

h ( i ) = ∑ d ∣ i g ( i d ) f ( d ) = i 2 ∑ d ∣ i ϕ ( d ) = i 3 h(i)=\sum_{d|i}g(\frac{i}{d})f(d)=i^2\sum_{d|i}\phi(d)=i^3 h(i)=dig(di)f(d)=i2diϕ(d)=i3

loj6229

目前遇到的最难的推式子题。。。

写一下推导的过程吧。。

要求 ∑ i = 1 n ∑ j = 1 i l c m ( i , j ) g c d ( i , j ) \sum_{i=1}^n\sum_{j=1}^i\frac{lcm(i,j)}{gcd(i,j)} i=1nj=1igcd(i,j)lcm(i,j)

考虑求 ∑ i = 1 n ∑ j = 1 n l c m ( i , j ) g c d ( i , j ) \sum_{i=1}^n\sum_{j=1}^n\frac{lcm(i,j)}{gcd(i,j)} i=1nj=1ngcd(i,j)lcm(i,j) 然后 + n +n +n再除以 2 2 2就是答案。

= ∑ i = 1 n ∑ j = 1 n i j g c d 2 ( i , j ) =\sum_{i=1}^n\sum_{j=1}^n\frac{ij}{gcd^2(i,j)} =i=1nj=1ngcd2(i,j)ij

= ∑ d = 1 n 1 d 2 ∑ i = 1 n ∑ j = 1 n i j [ g c d ( i , j ) = d ] =\sum_{d=1}^n\frac{1}{d^2}\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=d] =d=1nd21i=1nj=1nij[gcd(i,j)=d]

= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ i j [ g c d ( i , j ) = 1 ] =\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij[gcd(i,j)=1] =d=1ni=1dnj=1dnij[gcd(i,j)=1]

= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ i j ∑ d ′ ∣ g c d ( i , j ) μ ( d ′ ) =\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum_{d&#x27;|gcd(i,j)}\mu(d&#x27;) =d=1ni=1dnj=1dnijdgcd(i,j)μ(d)

= ∑ d = 1 n ∑ d ′ = 1 ⌊ n d ⌋ μ ( d ′ ) ( d ′ ) 2 ∑ i = 1 ⌊ n d d ′ ⌋ ∑ j = 1 ⌊ n d d ′ ⌋ i j =\sum_{d=1}^n\sum_{d&#x27;=1}^{\lfloor\frac{n}{d}\rfloor}\mu(d&#x27;)(d&#x27;)^2\sum_{i=1}^{\lfloor\frac{n}{dd&#x27;}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dd&#x27;}\rfloor}ij =d=1nd=1dnμ(d)(d)2i=1ddnj=1ddnij

q = d ′ d q=d&#x27;d q=dd

∑ q = 1 n ∑ d ′ ∣ q μ ( d ′ ) ( d ′ ) 2 ( ∑ i = 1 ⌊ n q ⌋ i ) 2 \sum_{q=1}^n\sum_{d&#x27;|q}\mu(d&#x27;)(d&#x27;)^2(\sum_{i=1}^{\lfloor\frac{n}{q}\rfloor}i)^2 q=1ndqμ(d)(d)2(i=1qni)2

f ( i ) = ∑ d ∣ i μ ( d ) d 2 f(i)=\sum_{d|i}\mu(d)d^2 f(i)=diμ(d)d2,发现后面那坨是个数论分块的东西。
所以任务就是求 f ( i ) f(i) f(i)的前缀和。

构造 g ( i ) = i 2 , t ( i ) = μ ( i ) i 2 g(i)=i^2,t(i)=\mu(i)i^2 g(i)=i2,t(i)=μ(i)i2,发现直接求 h = f ∗ g h=f*g h=fg不太好搞

接下来就是见证奇迹的时刻。

f = t ∗ 1 h = f ∗ g = t ∗ g ∗ 1 f=t*1\quad\quad h=f*g=t*g*1 f=t1h=fg=tg1

发现 t ∗ g = ϵ t*g=\epsilon tg=ϵ

证明:假设 l = t ∗ g l=t*g l=tg

l ( i ) = ∑ d ∣ i t ( d ) ∗ g ( i d ) = ∑ d ∣ i μ ( d ) i 2 = i 2 [ i = = 1 ] = ϵ l(i)=\sum_{d|i}t(d)*g(\frac{i}{d})=\sum_{d|i}\mu(d)i^2=i^2[i==1]=\epsilon l(i)=dit(d)g(di)=diμ(d)i2=i2[i==1]=ϵ

所以: h = t ∗ g ∗ 1 = 1 h=t*g*1=1 h=tg1=1

所以我们可以在 O ( 1 ) O(1) O(1)的时间内算出 H ( n ) = ∑ i = 1 n h ( i ) = n H(n)=\sum_{i=1}^nh(i)=n H(n)=i=1nh(i)=n

根据杜教筛的式子:

F ( n ) = H ( n ) − ∑ i = 2 n g ( i ) F ( n i ) g ( 1 ) F(n)=\frac{H(n)-\sum_{i=2}^ng(i)F(\frac{n}{i})}{g(1)} F(n)=g(1)H(n)i=2ng(i)F(in)

F ( n ) = H ( n ) − ∑ i = 2 n i 2 F ( n i ) F(n)=H(n)-\sum_{i=2}^ni^2F(\frac{n}{i}) F(n)=H(n)i=2ni2F(in),平方和的区间和可以通过两个通项公式倒腾出来。

复杂度的话,简单分析来就外层是一个大数论分块 O ( n 0.5 ) O(n^{0.5}) O(n0.5),里面杜教筛的 n n n依赖于外层数论分块的 L , R L,R L,R

这个东西可以归结成一类的复杂度分析:

∑ i = 1 n f ( n i ) × g ( i ) \sum_{i=1}^{n} f(\frac{n}{i}) \times g(i) i=1nf(in)×g(i),其中 g g g可以使用杜教筛求前缀和。

考虑在求 g g g的值的时候,杜教筛里面是有记忆化的。也就是一个数论分块,所以这次恰好会访问到那些上次外层调用杜教筛的地方。复杂度还是 O ( n 2 3 ) O(n^\frac{2}{3}) O(n32)

c o d e s : codes: codes:

#include <cstdio>
#include <map>
using namespace std;
typedef long long LL;
map<int, LL> q;
const int N = 1e6 + 5, P = 1e9 + 7;
int prime[N], vis[N], tot;
LL F[N], mu[N], inv2, inv6;
int n;
inline LL q_p(LL a, LL b, LL p) {
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}
inline void eluer() {
    vis[1] = 1; mu[1] = 1;
    inv2 = q_p(2, P - 2, P);
    inv6 = q_p(6, P - 2, P);
    for (int i = 2; i < N; ++i) {
        if (!vis[i]) {
            prime[++tot] = i;
            mu[i] = ((1 - 1LL * i * i) % P + P) % P;
        }
        for (int j = 1; j <= tot; ++j) {
            LL cur = 1LL * i * prime[j];
            if (cur >= N) break;
            vis[cur] = 1;
            if (i % prime[j] == 0) {
                mu[cur] = mu[i];
                break;
            }
            else mu[cur] = (mu[i] * mu[prime[j]]) % P;
        }
    }
    for (int i = 1; i < N; ++i) F[i] = (F[i - 1] + mu[i]) % P;
}
inline LL calc(int n) {
    return inv6 * n % P * (n + 1) % P * (2 * n + 1) % P;
}
inline LL dls(int n) {
    if (n < N) return F[n];
    if (q[n]) return q[n];
    LL res = n;
    for (int L = 2, R; L <= n; L = R + 1) {
        R = n / (n / L);
        res = ((res - ((calc(R) - calc(L - 1)) % P + P) % P * dls(n / L) % P) % P + P) % P;
    }
    return q[n] = res;
}
int main() {
    scanf("%d", &n);
    eluer();
    LL ans = 0;
    for (int L = 1, R; L <= n; L = R + 1) {
        R = n / (n / L);
        LL res = ((dls(R) - dls(L - 1)) % P + P) % P;
        res = res * (1 + (n / L)) % P * (1 + (n / L)) % P;
        res = res * (n / L) % P * (n / L) % P;
        res = res * inv2 % P * inv2 % P;
        ans = (ans + res) % P;
    }
    printf("%lld\n", (ans + n) % P * inv2 % P);
    return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值