min25筛

caoyang1123

已于 2024-10-13 06:18:27 修改

阅读量245

点赞数 1

分类专栏：笔记文章标签：拓扑学算法散列表

于 2021-08-13 23:20:37 首次发布

本文链接：https://blog.csdn.net/caoyang1123/article/details/119685375

版权

笔记专栏收录该内容

72 篇文章 0 订阅

订阅专栏

参考：【学习笔记】Min25筛_cz_xuyixuan的博客-CSDN博客_min25筛

求解积性函数前缀和，要求 $\sum [i \, is \, prime] * f(i)$ 能够由 $\sum [i \, is \, prime] * i^k$ 表示。

步骤：

首先对于每个x = [floor(N / i)]，求 $\sum_{i = 1}^x [i \, is \, prime] * i^k$ ，先预处理根号N以内的质数k次幂

前缀和，记为sum[i]，令g(N, i)表示埃氏筛第i轮没有被筛去的数的k次幂和，

可以发现所求就是i = cnt，cnt为根号N内质数个数。

可以发现

g(N,i) = g(N,i−1)，当prime[i]^2 > N

g(N,i) = g(N,i−1) − prime[i]^k∗(g(⌊N/prime[i]⌋,i−1)−sum[i−1]）

然后设S(N, i) 为1到N中最小质因数大于等于prime[i] 的f(i)之和，所求就是

S(N, 1) + f(1)

如何求S(N, i)

1.考虑质数贡献 $\sum_{j = 1}^N [j \, is \, a \, prime]f(j) - \sum_{j = 1}^{i - 1} f(prime_j)$

2.考虑合数贡献，通过枚举最小质因子以及出现幂次

$\sum_{j = i}^{prime_j^2<=N} \, \sum_{k = 1}^{prime_j^{k + 1}<=N} (S(\lfloor \frac{N}{prime_j^k}\rfloor,j+1) * f(prime_j^k) + f(prime_j^{k + 1} ))$

(第一个sigma下标是j = i，打错了)

模板题：

loj6053简单的函数

对于质数处的函数值，除了2是3，其它都是质数值-1，因此求质数上的函数值的和等价于质数和-质数个数，可以代入上述方法推。

#include<bits/stdc++.h>
#define pii pair<int,int>
#define fi first
#define sc second
#define pb push_back
#define ll long long
#define trav(v,x) for(auto v:x)
#define all(x) (x).begin(), (x).end()
#define VI vector<int>
#define VLL vector<ll>
using namespace std;
const int N = 1e6 + 100;
const ll mod = 1e9 + 7, inv2 = 5e8 + 4;

int pri[N], pnum = 0;
bool np[N];
ll sum[2][N], n, val[N], g[2][N];
int num = 0, ps1[N], ps2[N];

ll cal_s(ll x, int y)
{
    if(x < pri[y])
        return 0;
    ll res;
    int ps = (x <= n / x ? ps1[x] : ps2[n / x]);
    res = (g[1][ps] - g[0][ps] + mod) % mod;
    res = (res - sum[1][y - 1] + sum[0][y - 1] + mod) % mod;
    for(int j = y; j <= pnum && 1LL * pri[j] * pri[j] <= x; j++)
    {
        ll tmp = pri[j];
        for(ll k = 1; tmp * pri[j] <= x; tmp = tmp * pri[j], k++)
        {
            res = (res + cal_s(x / tmp, j + 1) * (pri[j] ^ k) + (pri[j] ^ (k + 1))) % mod;
        }
    }
    return res;
}

signed main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);
	for(int i = 2; i < N; i++)
    {
        if(!np[i])
        {
            pri[++pnum] = i;
            sum[0][pnum] = pnum;
            sum[1][pnum] = (sum[1][pnum - 1] + i) % mod;
        }
        for(int j = 1; j <= pnum && i * pri[j] < N; j++)
        {
            np[i * pri[j]] = 1;
            if(i % pri[j] == 0)
                break;
        }
    }
    cin >> n;
    for(ll l = 1, r = 0; l <= n; l = r + 1)
    {
        ll x = n / l;
        r = n / x;
        val[++num] = x;
        g[0][num] = (x - 1 + mod) % mod;
        g[1][num] = (x + 2) % mod * ((x - 1) % mod) % mod * inv2 % mod;
        if(x <= n / x)
            ps1[x] = num;
        else ps2[n / x] = num;
    }
    for(int i = 1; i <= pnum && 1LL * pri[i] * pri[i] <= n; i++)
    {
        for(int j = 1; j <= num && 1LL * pri[i] * pri[i] <= val[j]; j++)// j increase, val[j] is decrease
        {
            ll tmp = val[j] / pri[i];
            int bf = (tmp <= n / tmp ? ps1[tmp] : ps2[n / tmp]);
            g[0][j] = (g[0][j] - g[0][bf] + sum[0][i - 1] + mod) % mod;
            g[1][j] = (g[1][j] - pri[i] * (g[1][bf] - sum[1][i - 1] + mod) % mod + mod) % mod;
        }
    }
    int tmp = 1;
    if(n >= 2)
        tmp += 2;
    cout << (cal_s(n, 1) + tmp) % mod << '\n';
}