• min_25筛法亚线性时间内求解积性函数的前缀和


    min_25筛

    参考链接:

    https://www.luogu.com.cn/blog/wucstdio/solution-p5325

    https://www.cnblogs.com/zkyJuruo/p/13445596.html#处理整个函数的前缀和

    wucstdio同学讲的很清楚了,但是中间递推g的式子错了,应该是\(p_j^k(g(\lfloor\frac{n}{p_j}\rfloor, j-1)-g(p_{j-1}, j-1))\).

    min_25筛可以用来在亚线性时间内求解满足条件的积性函数的前缀和。

    条件为: 当\(p\)为质数时,\(f(p)\)是一个关于\(p\)的项数较少的多项式,或可以快速求值;\(f(p^c)\)可以快速求值

    \(S(n,j)=\sum_{i=2}^{n}f(i)\times[i的最小质因子比第j个质数大]\),那么我们要求的就是\(f(1)+S(n,0)\);其中\(f\)是满足条件的函数。

    \(S\)的求解,被分为了\(i\)为质数和\(i\)为合数俩部分。

    1.质数部分的前缀和

    假设完全积性函数\(f'\)在质数处取值与\(f\)相同。

    定义\(g(n,j)\)\([1,n]\)中用前\(j\)个质数执行埃拉斯托尼筛法后,还未被筛去的数对答案的贡献。

    \[g(n,j)=\sum_{i=1}^{n}f'(i)\times[i的最小质因子比第j个质数大或i为质数] \]

    那么质数部分对答案的贡献

    \[Sum\_Prime(n)=\sum_{i=1}^{n}f(i)\times[i\in P]=g(n,|P|) \]

    其中\(P\)表示\([1,n]\)中全体质数的集合。

    \(g(n,0)=\sum_{i=2}^{n}f'(i)\)

    \(g(n,|P|)\)\(f'\)在合数处的取值无关,故在中间递推时可以用\(f'\)代替\(f\)

    \(p_j\)为第\(j\)个质数。

    \(p_j^2>n\)时,显然有\(g(n,j)=g(n,j-1)\) (没有以\(p_j\)作为最小质因数的合数能够被筛去)

    \(p_j^2<=n\)时,在从\(g(n,j-1)\)\(g(n,j)\)的转移中新被筛去的数,是恰好以\(p_j\)作为最小质因数的合数。

    利用完全积性函数的性质把它提取出来,得到

    \[g(n,j)=g(n,j-1)-f'(p_j)(g(\lfloor\frac{n}{p_j}\rfloor,j-1)-g(p_{j-1},j-1)) \]

    提取了一个\(p_j\)后,剩下的数中合数的最小质因子都大于等于\(p_j\)

    \(g(\lfloor\frac{n}{p_j}\rfloor,j-1)\)还包括了前\(j-1\)个质数对答案的贡献,应该减去。

    最终得到了上述转移方程,能够用于计算\(g(n,|P|)\)

    我们还需要处理出所有的\(\lfloor\frac{i}{p_j}\rfloor\),这需要利用到以下性质:

    \[\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{n}{ab}\rfloor \]

    这意味着我们不需要求出转移过程中所有的\(n\),只需要根据整除分块的那套理论对\(n\)求出所有\(\lfloor\frac{n}{x}\rfloor\)即可。

    根据整除分块的那套理论,对于\(i\),满足\(\lfloor\frac{n}{i}\rfloor=k\)的所有的整数\(i\)所在区间的左端点记为\(l\)的话,右端点就是\(\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor=\lfloor\frac{n}{k}\rfloor\),记为\(r\)

    证明:

    先证当\(i\)\(\lfloor\frac{n}{k}\rfloor\)时满足\(\lfloor\frac{n}{i}\rfloor=k\)

    因为\(\lfloor\frac{n}{k}\rfloor\le\frac{n}{k}\)

    \(\lfloor\frac{n}{\lfloor\frac{n}{k}\rfloor}\rfloor\ge k\)

    再证\(r\ge i\):

    \[\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor\ge\lfloor\frac{n}{\frac{n}{i}}\rfloor=i \]

    再证\(r=\lfloor\frac{n}{k}\rfloor\)

    \[\lfloor\frac{n}{i}\rfloor=k\leftrightarrow k\le \frac{n}{i}< k+1\lrarr\frac{1}{k+1}<\frac{i}{n}\le\frac{1}{k}\lrarr\frac{n}{k+1}< i \le \frac{n}{k} \]

    又因为\(i\)是正整数,故\(r=\lfloor\frac{n}{k}\rfloor\)

    转移过程中的数都是\(\lfloor\frac{n}{x}\rfloor\)的形式,数量大概是\(O(\sqrt n)\)级别。故我们可以手动离散化一下,用\(idx1\)记录\(\lfloor\frac{n}{x}\rfloor\le \sqrt{n}\)时的\(\lfloor\frac{n}{x}\rfloor\)的下标,\(idx2\)记录另一种情况时的下标.

    2.合数部分的前缀和

    枚举每一个数的最小质因数及其次数,由于这里满足互质的条件,可以直接把这个质因数的幂提取出来。

    \[S(n,j)=\sum_{i=2}^{n}f(i)\times[i的最小质因子比第j个质数大] \]

    那么合数部分的前缀和为

    \[\sum_{k>j}\sum_{e=1}^{p_k^e \le n}f(p_k^e)(S(\lfloor\frac{n}{p_k^e}\rfloor, k)+[e>1]) \]

    \([e>1]\)表示\(p_k^e\)自身对答案的贡献(\(S\)没有计算f(1)的贡献,并且质数部分已经计算过了)

    那么

    \[S(n,j) = \sum_{i=1}^{n}[i\in P]\times f(i)-\sum_{i=1}^{j}f(p_i)+\sum_{k>j}\sum_{e=1}^{p_k^e \le n}f(p_k^e)(S(\lfloor\frac{n}{p_k^e}\rfloor, k)+[e>1]) \]

    那么最终答案是\(S(n,0)+f(1)\)

    \[g(n,|P|)+\sum_{k>0}\sum_{e=1}^{p_k^e \le n}f(p_k^e)(S(\lfloor\frac{n}{p_k^e}\rfloor, k)+[e>1]) \]

    实操:洛谷P5325 【模板】Min_25筛

    符合条件的函数为\(f(x)\),且\(f(p^k)=p^k(p^k-1)\)

    把多项式\(f(p^k)\)表示为若干个单项式的和,得

    \[f(p^k)=p^{2k}-p^k \]

    依次取\(f'(x)=x^2,f'(x)=x\).

    在求解\(g\)时,我们要计算\(f'\)在小于等于\(\sqrt n\)的质数处的前缀和\(g(p_{j-1},j-1)\),分别将一次幂和二次幂的前缀和记为sum1和sum2,这个可以在线性筛的时候顺便完成。

    然后用整除分块求出所有\(\lfloor\frac{n}{x}\rfloor = k\),并离散化存一下.

    对于所有的\(k\), 我们需要求出\([1,k]\)的前缀和。

    \(f'\)拆分成一次项和二次项,分别求前缀和即可。但这个前缀和不能递推,需要用等差数列求和公式以及二次幂求和公式快速计算出\(k\)处的前缀和。

    求解\(g\)的时候要注意,需要保证\(f'\)是完全积性函数,也就是说,必须分别求出一次幂作为\(f'\)和二次幂作为\(f'\)时的\(g\),不妨分别设为\(g_1,g_2\)

    注意\(l,r\)得为\(long long\)

    首先整除分块并求出\(g(i, 0)\)

        for (ll l = 1, r; l <= n; l = r + 1) {
            ll k = n / l;
            r = n / k;
            w[++totw] = k;
            //求g(i,0)
            //f'(i)=i^2-i
            //g(i,0) = \sum_{j=2}^{n}f'(j) = 一次幂的前缀和减去二次幂的前缀和减去f'(1)
            //这题f'(1)是0 不用减了。。。
            ll tk = k % md;
            g2[totw] = tk * (tk+1) % md * inv2 % md * (2*tk+1) % md * inv3 % md;
            g2[totw]--;
            g1[totw] = (tk + 1ll) * tk % md * inv2 % md; 
            g1[totw]--;    
            if (k <= sqrtn) {
                idx1[k] = totw;
            } else idx2[r] = totw;
        }
    

    然后\(dp\)求出\(g(i, |P|)\)

    for (int i = 1; i <= tot; i++) {
            for (int j = 1; j <= totw; j++) {
                if (pri[i] * pri[i] > w[j]) break;
                ll k = w[j] / pri[i];
                ll pre = k <= sqrtn ? idx1[k] : idx2[n / k];
                g2[j] = (g2[j] - pri[i] * pri[i] % md * (g2[pre] - sum2[i-1] + md) % md + md) % md;
                g1[j] = (g1[j] - pri[i] * (g1[pre]-sum1[i-1] + md) % md + md) % md;
            }
        }
    

    接下来求\(S\)

    ll S(ll i, ll j) {
        //printf("i == %lld j == %lld\n", i, j);
        //S(n,j)=\sum_{i=2}^{n}f(i)[i的最小质因子比第j个质数大]
        if (pri[j] >= i) 
            return 0;
        //质数部分的前缀和,pri[j]在sqrtn以内,一定是由idx1表示。
        int id = i > sqrtn ? idx2[n/i] : idx1[i];
        ll sp = ((g2[id] - g1[id] + sum1[j] - sum2[j]) % md + md) % md;
        for (ll k = j + 1; k <= tot && pri[k] * pri[k] <= i; k++) { //超过sqrti的对答案没贡献
            for (ll p = pri[k]; p <= i; p *= pri[k]) {
                ll tmp = S(i/p, k) + (p > pri[k]);
                ll tp = p % md;
                sp = (sp + tp * (tp - 1 + md) % md * tmp) % md;
            }
        }
        return sp;
    }
    

    最终目标是求出\(S(n, 0) + f(1)\),由积性函数性质知\(f(1)=1\),故答案为\((S(n, 0) + 1ll) \mod p\)

    洛谷ac代码如下

    #include<bits/stdc++.h>
    using namespace std;
    #define ll long long
    const ll maxn = 1e6 + 7, md = 1e9 + 7;
    bool isp[maxn];
    ll n, pri[maxn], idx1[maxn], idx2[maxn], sum1[maxn], sum2[maxn], w[maxn], sqrtn, g1[maxn], g2[maxn], inv2, inv3;
    int tot, totw;
    void sieve () {
        memset(isp, true, sizeof(isp));
        isp[1] = 0;
        for (int i = 2; i <= 1000000; i++) {
            if (isp[i]) {
                pri[++tot] = i;
                sum1[tot] = (sum1[tot-1] + pri[tot]) % md;
                sum2[tot] = (sum2[tot-1] + pri[tot] * pri[tot] % md) % md;
            }
            for (int j = 1; j <= tot && i * pri[j] <= 1000000; j++) {
                isp[i*pri[j]] = 0;
                if (i % pri[j] == 0) 
                    break;
            }
        }
    }
    ll ksm(ll a, int b) {
        ll res = 1;
        while (b) {
            if (b & 1) res = res * a % md;
            a = a * a % md;
            b >>= 1;
        }
        return res % md;
    }
    ll S(ll i, ll j) {
        //printf("i == %lld j == %lld\n", i, j);
        //S(n,j)=\sum_{i=2}^{n}f(i)[i的最小质因子比第j个质数大]
        if (pri[j] >= i) 
            return 0;
        //质数部分的前缀和,pri[j]在sqrtn以内,一定是由idx1表示。
        int id = i > sqrtn ? idx2[n/i] : idx1[i];
        ll sp = ((g2[id] - g1[id] + sum1[j] - sum2[j]) % md + md) % md;
        for (ll k = j + 1; k <= tot && pri[k] * pri[k] <= i; k++) { //超过sqrti的对答案没贡献
            for (ll p = pri[k]; p <= i; p *= pri[k]) {
                ll tmp = S(i/p, k) + (p > pri[k]);
                ll tp = p % md;
                sp = (sp + tp * (tp - 1 + md) % md * tmp) % md;
            }
        }
        return sp;
    }
    int main() {
        cin >> n;
        sqrtn = sqrt(n);
        inv2 = ksm(2, md - 2);
        inv3 = ksm(3, md - 2);
        sieve();
        for (ll l = 1, r; l <= n; l = r + 1) {  
            ll k = n / l;
            //printf("k == %lld\n", k);
            r = n / k;
            w[++totw] = k;
            //求g(i,0)
            //f'(i)=i^2-i
            //g(i,0) = \sum_{j=2}^{n}f'(j) = 一次幂的前缀和减去二次幂的前缀和减去f'(1)
            //这题f'(1)是0 不用减了。。。
            ll tk = k % md;
            g2[totw] = tk * (tk+1) % md * inv2 % md * (2*tk+1) % md * inv3 % md;
            g2[totw]--;
            g1[totw] = (tk + 1ll) * tk % md * inv2 % md; 
            g1[totw]--;    
            if (k <= sqrtn) {
                idx1[k] = totw;
            } else idx2[r] = totw;
        }
        //printf("tot == %d totw == %d\n", tot, totw);
        //g数组滚动掉了第二维,因此要先枚举第二维
        for (int i = 1; i <= tot; i++) {
            for (int j = 1; j <= totw; j++) {
                if (pri[i] * pri[i] > w[j]) break;
                ll k = w[j] / pri[i];
                ll pre = k <= sqrtn ? idx1[k] : idx2[n / k];
                g2[j] = (g2[j] - pri[i] * pri[i] % md * (g2[pre] - sum2[i-1] + md) % md + md) % md;
                g1[j] = (g1[j] - pri[i] * (g1[pre]-sum1[i-1] + md) % md + md) % md;
            }
        }
        cout << (S(n, 0) + 1ll) % md;
        return 0;
    }
    
  • 相关阅读:
    Codeforces Round #592 (Div. 2)C. The Football Season(暴力,循环节)
    Educational Codeforces Round 72 (Rated for Div. 2)D. Coloring Edges(想法)
    扩展KMP
    poj 1699 Best Sequence(dfs)
    KMP(思路分析)
    poj 1950 Dessert(dfs)
    poj 3278 Catch That Cow(BFS)
    素数环(回溯)
    sort与qsort
    poj 1952 buy low buy lower(DP)
  • 原文地址:https://www.cnblogs.com/YjmStr/p/15003872.html
Copyright © 2020-2023  润新知