【min_25筛】筛欧拉函数、莫比乌斯函数、积性函数
传送门:欧拉函数、莫比乌斯函数
花了一下午学了一下min25筛。
一知半解。不过还是学会了套板子。
mu函数要看成质数个数取反,在预处理g的时候不能带负数(我也不知道为什么)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e5 + 233;
ll primes[maxn], cnt;
bool st[maxn];
ll sp1[maxn], sp2[maxn], g1[maxn], g2[maxn], sm[maxn], g[maxn], w[maxn], wcnt;
ll id1[maxn], id2[maxn];
ll inv2, inv6, qn;
ll ksm(ll n, ll k)
{
ll res = 1;
while(k)
{
if(k & 1) res = res * n;
n = n * n ;
k >>= 1;
}
return res;
}
void get_primes(int n)
{
for(int i = 2; i <= n; i++) {
if(!st[i]) {
primes[++cnt] = i;
}
for(int j = 1; j <= cnt && i * primes[j] <= n; j++) {
st[i * primes[j]] = 1;
if(i % primes[j] == 0)
break;
}
}
for(int i = 1; i <= cnt; i++) {
sp2[i] = primes[i] + sp2[i - 1], sp1[i] = sp1[i - 1] + 1;
sm[i] = i;
}
}
ll n;
ll dfs(int x, int y)
{
if(primes[y] >= x) return 0;
int k = x;
if(k <= qn) k = id1[x]; else k = id2[n / x];
ll res = (g2[k] - g1[k]) - (sp2[y] - sp1[y]);
for(int i = y + 1; i <= cnt && primes[i] * primes[i] <= x; i++)
{
ll p = primes[i];
for(int e = 1; p <= x; e++, p *= primes[i])
{
res = (res + (p / primes[i] * (primes[i] - 1)) * (dfs(x / p, i) + (e != 1)));
}
}
return res;
}
ll dfs2(int x, int y)
{
if(primes[y] >= x) return 0;
int k = x;
if(k <= qn) k = id1[x]; else k = id2[n / x];
ll res = (g[k]) - (-sm[y]);
for(int i = y + 1; i <= cnt && primes[i] * primes[i] <= x; i++)
{
ll p = primes[i];
res = (res - (dfs2(x / p, i)));
}
return res;
}
int main()
{
get_primes(1e5 + 10);
int T; cin >> T;
while(T--)
{
cin >> n;
qn = sqrt(n);
// memset(g1, 0, sizeof g1);
// memset(g2, 0, sizeof g2);
// memset(g, 0, sizeof g);
// memset(w, 0, sizeof w);
wcnt = 0;
for(ll l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
w[++wcnt] = n / r;
ll t = (n / r);
g1[wcnt] = t - 1; //g(1)排除
g2[wcnt] = (t + 1) * t / 2 - 1;
g[wcnt] = t - 1;
if(n / r <= qn) id1[n / r] = wcnt; else id2[r] = wcnt;
}
for(int i = 1; i <= cnt; i++)
{
ll pp = primes[i] * primes[i];
for(int j = 1; j <= wcnt && pp <= w[j]; j++)
{
ll k = w[j] / primes[i];
if(k <= qn) k = id1[k]; else k = id2[n / k];
g2[j] = g2[j] - primes[i] * (g2[k] - sp2[i - 1]);
g1[j] = g1[j] - (g1[k] - sp1[i - 1]);
g[j] = g[j] - (g[k] - sm[i - 1]);
}
}
for(int i = 1; i <= wcnt; i++) g[i] = -g[i];
// for(int i = 1; i <= wcnt; i++) cerr << g[i] << endl;
printf("%lld %lld\n", dfs(n, 0) + 1, dfs2(n, 0) + 1);
}
}
传送门:min25模板题
只会套板子。。
// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 2e5 + 233;
ll mod = 1e9 + 7;
ll primes[maxn], phi[maxn], cnt;
ll w[maxn], wcnt; // 下标
bool st[maxn];//int st[N];
ll s_p1[maxn], s_p2[maxn], g1[maxn], g2[maxn], id1[maxn], id2[maxn];
// f(p) = p*p - p
ll inv2, inv6, qn;
ll ksm(ll n, ll k)
{
ll res = 1;
while(k)
{
if(k & 1) res = res * n % mod;
n = n * n % mod;
k >>= 1;
}
return res;
}
void get_primes(int n)
{
// mu[1] = phi[1] = 1;
phi[1] = 1;
for(int i = 2; i <= n; i++) {
if(!st[i]) {
primes[++cnt] = i;
phi[i] = i - 1;
// mu[i] = -1;
}
for(int j = 1; j <= cnt && i * primes[j] <= n; j++) {
st[i * primes[j]] = 1;
if(i % primes[j] == 0) {
phi[i * primes[j]] = phi[i] * primes[j];
// mu[i * primes[j]] = 0;
break;
}else {
phi[i * primes[j]] = phi[i] * (primes[j] - 1);
// mu[i * primes[j]] = -mu[i];
}
}
}
for(int i = 1; i <= cnt; i++)
s_p1[i] = (primes[i] + s_p1[i - 1]) % mod, s_p2[i] = (s_p2[i - 1] + primes[i] * primes[i]) % mod;
}
ll n;
ll dfs(ll x, int y)
{
if(primes[y] >= x) return 0;
ll k = x;
if(k <= qn) k = id1[x]; else k = id2[n / x];
ll res = ((g2[k] - g1[k]) - (s_p2[y] - s_p1[y])) % mod;
if(res < 0) res += mod;
for(int i = y + 1; i <= cnt && primes[i] * primes[i] <= x; i++)
{
ll p = primes[i];
for(int e = 1; p <= x; e++, p *= primes[i])
{
ll now = p % mod;
res = (res + now * (now - 1LL) % mod * (dfs(x / p, i) % mod + (e != 1)) % mod ) % mod;
}
}
return res;
}
int main()
{
inv2 = ksm(2, mod - 2); inv6 = ksm(6, mod - 2);
cin >> n;
qn = sqrt(n);
get_primes(1e5 + 10);
for(ll l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
w[++wcnt] = n / r;
ll t = (n / r) % mod;
g1[wcnt] = (t + 1) * t % mod * inv2 % mod - 1; //g(1)排除
g2[wcnt] = (2 * t + 1) * (t + 1) % mod * t % mod * inv6 % mod - 1;
if(n / r <= qn) id1[n / r] = wcnt; else id2[r] = wcnt;
}
for(int i = 1; i <= cnt; i++)
{
for(int j = 1; j <= wcnt && primes[i] * primes[i] <= w[j]; j++)
{
ll k = w[j] / primes[i];
if(k <= qn) k = id1[k]; else k = id2[n / k];
g1[j] = (g1[j] - primes[i] * (g1[k] - s_p1[i - 1]) % mod) % mod;
g2[j] = (g2[j] - primes[i] * primes[i] % mod * (g2[k] - s_p2[i - 1]) % mod) % mod;
if(g1[j] < 0) g1[j] += mod;
if(g2[j] < 0) g2[j] += mod;
}
}
printf("%lld\n", dfs(n, 0) + 1LL);
}
发表评论