《Min_25筛》

其实在上一年就已经学过了Min_25，但当时理解得不是很好，现在重新来写一下。（其实是洲阁筛找不到阳间的板子，所以还是果断回来Min_25）

好了回归正题，其实也没有什么好讲的好像，就是推公式就行了。（雾）

可以使用Min25筛的前提是，该函数是个积性函数，且它在质数处的前缀和可以很快求出来，而且是关于质数p的一个低阶多项式。

积性函数：f(a) * f(b) = f(ab)，在a，b互质时。

$g(n,i) = \left\{\begin{matrix}
g(n,i - 1) & & (p[i-1])^{2} > n \\
g(n,i - 1) - F(p[i]) * (g(\left \lfloor \frac{n}{p[i]} \right \rfloor,i - 1) - \sum_{j = 1}^{i - 1}F(p[j]) )& & (p[i-1])^{2} <= n
\end{matrix}\right.$

$S(n,i) = g(n,|P|) - \sum_{j = 1}^{i - 1}f(p[j]) + \sum_{i = j}^{p[i] ^ 2 <= n}\sum_{e = 1}^{p[i] ^ {(e + 1)} <= n} (f(p[i]^e) * S(\left \lfloor \frac{n}{p[i]^e} \right \rfloor),i + 1) + f([p[i]]^{e+1}) )$

接下来是我个人的一些理解：

这里的F[i]函数并不是题目中给出的函数，而是我们要去找到的和f(i)在质数处取值相同的低阶的多项式。

而且这里系数要为1，也就是如果有系数那么就最后再乘上。

然后就是对于G的推导，这里的f就是题目中给出的f多项式，我们的处理是没有处理到f(1)的值的，所以最后要去加上。

关于g的预处理$g = \sum_{i = 2}^{n}F(i)$如果有系数那么先推出没有系数的最后在每次算S的时候，g都要乘上这个系数。

百做不厌了属于是.

假定当前函数为f(i) = 1，那么就是求i = prime,即g那部分。

所以这里只需要求g部分的值就可以了。

一点细节：此处f(1) = 0没有贡献，所以最后不需要+1.

// Author: levil
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef pair<int,int> pii;
const int N = 1e6 + 5;
const int M = 1e3 + 5;
const double eps = 1e-10;
const LL Mod = 1e9 + 7;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;

bool vis[N];
int tot = 0,id1[N],id2[N],sq,sw = 0;
LL sum[N],prime[N],g[N],n,w[N];
void init() {
    sq = sqrt(n);
    for(int i = 2;i <= sq;++i) {
        if(!vis[i]) {
            prime[++tot] = i;
            sum[tot] = sum[tot - 1] + 1;
        }
        for(int j = 1;j <= tot && prime[j] * i <= sq;++j) {
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void get_g() {
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        w[++sw] = n / L;//块的值
        g[sw] = w[sw] - 1;
        if(w[sw] <= sq) id1[w[sw]] = sw;
        else id2[r] = sw;
    }
    for(int i = 1;i <= tot;++i) {
        for(int j = 1;j <= sw && w[j] >= prime[i] * prime[i];++j) {
            LL p = w[j] / prime[i];
            p = (p <= sq ? id1[p] : id2[n / p]);
            g[j] = g[j] - 1LL * (g[p] - sum[i - 1]);
        }
    }
}
void solve() {
    scanf("%lld",&n);
    init();
    get_g();
    int p = (n <= sq) ? id1[n] : id2[n / n];
    LL ans = g[p] - 0; 
    printf("%lld\n",ans);
}
int main() {
    solve();
    system("pause");
    return 0;
}

View Code

这题乍一看没什么思路，但是仔细看可以发现，对于c = 1的时候就是各个素数，那么我们知道除了2之外所有的素数都是奇数，即p ^ 1 = p - 1

所以对于除二之外的所有素数 f[i] = i ^ 1 - i ^ 0，这里也要分成两段来处理。

然后我们在底层y = 1有2的影响的时候加上2的贡献即可。

// Author: levil
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef pair<int,int> pii;
const int N = 5e5 + 5;
const int M = 1e3 + 5;
const double eps = 1e-10;
const LL Mod = 1e9 + 7;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;

bool vis[N];
int tot = 0,id1[N],id2[N],sq,sw = 0;
LL sum1[N],sum2[N],prime[N],g1[N],g2[N],n,inv6 = 166666668,inv2 = 500000004,w[N];
LL ADD(LL x,LL y) {return (x + y) % Mod;}
LL MUL(LL x,LL y) {return x * y % Mod;}
LL DEC(LL x,LL y) {return ((x - y) % Mod + Mod) % Mod;}
void init() {
    sq = sqrt(n);
    for(int i = 2;i <= sq;++i) {
        if(!vis[i]) {
            prime[++tot] = i;
            sum1[tot] = ADD(sum1[tot - 1],i);
            sum2[tot] = ADD(sum2[tot - 1],1);
        }
        for(int j = 1;j <= tot && prime[j] * i <= sq;++j) {
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void get_g() {
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        w[++sw] = n / L;//块的值
        LL x = w[sw] % Mod;
        g1[sw] = MUL(MUL(x,x + 1),inv2);//预处理边界 - g(x,0)
        g1[sw] = DEC(g1[sw],1);
        g2[sw] = DEC(x,1);
        if(n / r <= sq) id1[n / r] = sw;
        else id2[r] = sw;
    }
    for(int i = 1;i <= tot;++i) {
        LL val = 1LL * prime[i] * prime[i];
        for(int j = 1;j <= sw && w[j] >= val;++j) {
            LL p = w[j] / prime[i];
            p = (p <= sq ? id1[p] : id2[n / p]);
            g1[j] = DEC(g1[j],MUL(prime[i],DEC(g1[p],sum1[i - 1])));//dp转移
            g2[j] = DEC(g2[j],MUL(1,DEC(g2[p],sum2[i - 1])));
        }
    }
}
LL get_S(LL x,int y) {
    if(prime[y] > x || x <= 1) return 0;
    int p = (x <= sq) ? id1[x] : id2[n / x];
    LL res = DEC(DEC(g1[p],g2[p]),DEC(sum1[y - 1],sum2[y - 1]));
    if(y == 1) res = ADD(res,2);
    for(int i = y;i <= tot && prime[i] * prime[i] <= x;++i) {
        for(LL e = 1,sp = prime[i],pp = prime[i] * prime[i];pp <= x;sp *= prime[i],pp *= prime[i],e++) {
            LL ma1 = (prime[i] ^ e) % Mod;
            LL ma2 = (prime[i] ^ (e + 1)) % Mod;
            res = ADD(res,ADD(MUL(ma1,get_S(x / sp,i + 1)),ma2));
        }
    }
    return res;
}
void solve() {
    scanf("%lld",&n);
    init();
    get_g();
    LL ans = get_S(n,1) + 1;
    printf("%lld\n",ans % Mod);
}
int main() {
    solve();
    system("pause");
    return 0;
}

View Code

UOJ：188. 【UR #13】Sanrd

求区间次大质因子的和，即求

虽然对于质数的值 = 0，很容易得出，但是在这里这个函数并不是积性函数，所以不好操作。

但是我们可以从Min25的推导出发，由S函数的推导，我们可以去新建一种dp的思路。

因为我们枚举的是minp > pj，和e次，所以对于每个段，它要满足pj是次大质因子的话。

那每个pj^k的贡献就是$\sum_{p[j]}^{\left \lfloor \frac{n}{p[j]} \right \rfloor} i (i~is~a~ prime)$

且这里的f（pj ^ k） = pj。

所以可得转移式子$S(n,i) = S(n / pj,k + 1) + p[j] * \sum_{i = p[j]}^{\left \lfloor \frac{n}{p[j]} \right \rfloor}(i~is~a~ prime)$

// Author: levil
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef pair<int,int> pii;
const int N = 1e6 + 5;
const int M = 1e3 + 5;
const double eps = 1e-10;
const LL Mod = 1e9 + 7;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;

bool vis[N];
int tot = 0,id1[N],id2[N],sq,sw = 0;
LL sum[N],prime[N],g[N],L,r,w[N];
void init(LL n) {
    sq = sqrt(n);
    tot = sw = 0;
    for(int i = 2;i <= sq;++i) {
        if(!vis[i]) {
            prime[++tot] = i;
            sum[tot] = sum[tot - 1] + 1;
        }   
        for(int j = 1;j <= tot && prime[j] * i <= sq;++j) {
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void get_g(LL n) {
    memset(g,0,sizeof(g));
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        w[++sw] = n / L;//块的值
        LL x = w[sw];
        g[sw] = x - 1;
        if(n / r <= sq) id1[n / r] = sw;
        else id2[r] = sw;
    }
    for(int i = 1;i <= tot;++i) {
        LL val = 1LL * prime[i] * prime[i];
        for(int j = 1;j <= sw && w[j] >= val;++j) {
            LL p = w[j] / prime[i];
            p = (p <= sq ? id1[p] : id2[n / p]);
            g[j] = g[j] - (g[p] - sum[i - 1]);
        }
    }
}
LL get_S(LL x,int y,LL n) {
    if(prime[y] > x || x <= 2) return 0;
    LL res = 0;
    for(int i = y;i <= tot && prime[i] * prime[i] <= x;++i) {
        for(LL e = 1,sp = prime[i],pp = prime[i] * prime[i];pp <= x;sp *= prime[i],pp *= prime[i],e++) {
            LL p = x / sp;
            p = (p <= sq ? id1[p] : id2[n / p]);
            res += get_S(x / sp,i + 1,n) + prime[i] * (g[p] - sum[i - 1]);
        }
    }
    return res;
}
LL cal(LL n) {
    init(n);
    get_g(n);
    return get_S(n,1,n);
}
void solve() {
    scanf("%lld %lld",&L,&r);
    LL ans = cal(r) - cal(L - 1);
    printf("%lld\n",ans);
}
int main() {
    solve();
    system("pause");
    return 0;
}

View Code

我们可以想到对于每个i，他有d(i)个因子，那么显然j 和 j + k都要是它的两个因子，且不能一样，即C(d(i),2)

所以本题即是求$ans = \sum_{i = 1}^{n}C(d(i),2)$

化简一下$ans = \sum_{i = 1}^{n}C(d(i),2) = \sum_{i = 1}^{n} d(i) * (d(i) - 1) / 2 = \frac{1}{2}\sum_{i = 1}^{n} d(i)^{2} - d(i)$

可以发现，对于前后两个都是积性函数，那么就可以进行两次Min25筛即可。

但是我们分析后面的可以发现，我们转化成去求每个i有几个倍数，即为$\sum_{i = 1}^{n}\left \lfloor \frac{n}{i} \right \rfloor$

那么我们对后面部分只需要进行一次整除分块即可了。

这里可以发现对于f(p) = 4 = 4 * (p ^ 0) ，所以我们处理除了p ^ 0的值，算的时候还是g(n)要乘上4.

// Author: levil
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef pair<int,int> pii;
const int N = 5e5 + 5;
const int M = 1e3 + 5;
const double eps = 1e-10;
const LL Mod = 998244353;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;

bool vis[N];
int tot = 0,id1[N],id2[N],sq,sw = 0;
LL prime[N],g[N],n,inv2 = 500000004,w[N];
LL ADD(LL x,LL y) {return (x + y) % Mod;}
LL MUL(LL x,LL y) {return x * y % Mod;}
LL DEC(LL x,LL y) {return ((x - y) % Mod + Mod) % Mod;}
void init() {
    sq = sqrt(n);
    for(int i = 2;i <= sq;++i) {
        if(!vis[i]) {
            prime[++tot] = i;
        }
        for(int j = 1;j <= tot && prime[j] * i <= sq;++j) {
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void get_g() {
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        w[++sw] = n / L;//块的值
        LL x = w[sw] % Mod;
        g[sw] = DEC(x,1);
        if(n / r <= sq) id1[n / r] = sw;
        else id2[r] = sw;
    }
    for(int i = 1;i <= tot;++i) {
        LL val = 1LL * prime[i] * prime[i];
        for(int j = 1;j <= sw && w[j] >= val;++j) {
            LL p = w[j] / prime[i];
            p = (p <= sq ? id1[p] : id2[n / p]);
            g[j] = DEC(g[j],MUL(1,DEC(g[p],i - 1)));//dp转移
        }
    }
}
LL get_S(LL x,int y) {
    if(prime[y] > x || x <= 1) return 0;
    int p = (x <= sq) ? id1[x] : id2[n / x];
    LL res = DEC(g[p],DEC(y,1));//这里是到y - 1的前缀和.
    res = MUL(res,4);
    for(int i = y;i <= tot && prime[i] * prime[i] <= x;++i) {
        for(LL e = 1,sp = prime[i],pp = prime[i] * prime[i];pp <= x;sp *= prime[i],pp *= prime[i],e++) {
            LL f1 = MUL(e + 1,e + 1);
            LL f2 = MUL(e + 2,e + 2);
            res = ADD(res,ADD(MUL(f1,get_S(x / sp,i + 1)),f2));
        }
    }
    return res;
}
LL block(LL n) {
    LL sum = 0;
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        LL len = DEC(r,L) + 1;
        sum = ADD(sum,MUL(len,n / L));
    }
    return sum;
}
LL quick_mi(LL a,LL b) {
    LL re = 1;
    while(b) {
        if(b & 1) re = re * a % Mod;
        a = a * a % Mod;
        b >>= 1;
    }
    return re;
}
void solve() {
    scanf("%lld",&n);
    init();
    get_g();
    LL ans = DEC(get_S(n,1) + 1,block(n));
    ans = ans * quick_mi(2,Mod - 2) % Mod;
    printf("%lld\n",ans % Mod);
}
int main() {
    solve();
    system("pause");
    return 0;
}

View Code

我们对题目的式子进行莫比乌斯反演

$ans = \sum_{d = 1}^{n}\sum_{i = 1}^{n}\sum_{j = 1}^{n}f(d)^{k} * [gcd(i,j) = d] = \sum_{d = 1}^{n}\sum_{i = 1}^{\frac{n}{d}}\sum_{j = 1}^{\frac{n}{d}}f(d)^{k} *[gcd(i,j) = 1] = \sum_{d = 1}^{n} f(d)^{k} \sum_{t = 1}^{\frac{n}{d}} \mu (t) \left \lfloor \frac{n}{dt} \right \rfloor^{2}$

令K = td

$\sum_{d = 1}^{n}f(d) ^ k\sum_{k = 1}^{n}\mu (\frac{k}{d})\left \lfloor \frac{n}{k} \right \rfloor ^ 2 = \sum_{k = 1}^{n}\left \lfloor \frac{n}{k} \right \rfloor ^ 2 \sum_{d | k}^{}f(d) ^ k \mu (\frac{k}{d})$

这里为什么d的枚举会变成d | k，因为本来我k那一个地方枚举的是d的倍数，所以保证了<=[n /d],现在它扩展到了1 ~ n，那就是去枚举倍数，所以我们的d就要变成去枚举因子了。

可以发现对于后面的东西，就是一个狄利克雷卷积。我们设F[x] = f(x) ^ k。

那么就是$F * \mu $很显然我们可以卷上一个恒等函数I。那么得$F * \mu * I = F * \epsilon = F$

我们一开始用了后面两个去推然后发现杜教筛不了，感觉自己像个傻子。

我们令$S(n) = \sum_{i = 1}^{n}\sum_{d | i}^{}f(d) ^ k\mu (\frac{i}{d})$

显然这个东东就是$F * \mu $所以最后的式子就是$S(n) = \sum_{i = 1}^{n}F(i) - \sum_{i = 2}^{n}I * S(n / i) = \sum_{i = 1}^{n}F(i) - \sum_{i = 2}^{n}S(n / i)$

可以发现前面这个F的前缀和我们可以Min25筛出，剩下的部分杜教筛即可。

有一个细节就是这里的f(x)和UOJ那题有点不一样就是对于质数它的值应该是1，而之前是0，所以我们还是当成0来处理但是最后要去加上质数的个数。

// Author: levil
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned int ULL;
typedef pair<int,int> pii;
const int N = 1e5 + 5;
const int M = 1e3 + 5;
const double eps = 1e-10;
const LL Mod = 1e9 + 7;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;

bool vis[N];
int tot = 0,id1[N],id2[N],sq,sw = 0;
LL n,k;
ULL g[N],prime[N],w[N],f[N];
ULL quick_mi(ULL a,LL b) {
    ULL re = 1;
    while(b) {
        if(b & 1) re = re * a;
        b >>= 1;
        a = a * a;
    }
    return re;
}
void init() {
    sq = sqrt(n);
    for(int i = 2;i <= sq;++i) {
        if(!vis[i]) {
            prime[++tot] = i;
            f[tot] = quick_mi(i,k);
        }
        for(int j = 1;j <= tot && prime[j] * i <= sq;++j) {
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void get_g() {
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        w[++sw] = n / L;//块的值
        LL x = w[sw];
        g[sw] = x - 1;
        if(n / r <= sq) id1[n / r] = sw;
        else id2[r] = sw;
    }
    for(int i = 1;i <= tot;++i) {
        LL val = prime[i] * prime[i];
        for(int j = 1;j <= sw && w[j] >= val;++j) {
            LL p = w[j] / prime[i];
            p = (p <= sq ? id1[p] : id2[n / p]);
            g[j] = g[j] - (g[p] - (i - 1));
        }
    }
}
ULL get_S(LL x,int y) {
    if(prime[y] > x || x <= 1) return 0;
    ULL res = 0;
    for(int i = y;i <= tot && prime[i] * prime[i] <= x;++i) {
        for(LL e = 1,sp = prime[i],pp = prime[i] * prime[i];pp <= x;sp *= prime[i],pp *= prime[i],e++) {
            LL p = x / sp;
            p = (p <= sq ? id1[p] : id2[n / p]);
            res = res + get_S(x / sp,i + 1) + f[i] * (g[p] - (i - 1));
        }
    }
    return res;
}
unordered_map<LL,ULL> mp;
ULL cal(LL x) {
    if(mp.count(x)) return mp[x];
    int p = (x <= sq ? id1[x] : id2[n / x]);
    ULL ans = get_S(x,1) + g[p];
    for(LL L = 2,r = 0;L <= x;L = r + 1) {
        r = x / (x / L);
        ans -= cal(x / L) * (r - L + 1);
    }
    return mp[x] = ans;
}
void solve() {
    scanf("%lld %lld",&n,&k);
    init();
    get_g();
    ULL ans = 0;
    for(LL L = 1,r;L <= n;L = r + 1) {
        r = n / (n / L);
        ans += 1ULL * (n / L) * (n / L) * (cal(r) - cal(L - 1));
    }
    printf("%u\n",ans);
}
int main() {
    solve();
    system("pause");
    return 0;
}

View Code

posted @ 2021-08-24 16:37 levill 阅读(65) 评论(0) 编辑收藏举报

刷新页面返回顶部

levill

《Min_25筛》

LOJ：6235. 区间素数个数

LOJ：6053. 简单的函数

UOJ：188. 【UR #13】Sanrd

LOJ：6682. 梦中的数论

LOJ：572. 「LibreOJ Round #11」Misaka Network 与求和

公告