可能是一篇（抄来的）min25学习笔记
一个要求很多的积性函数
我们考虑有一个积性函数，这个函数满足可以快速计算质数处的值
且质数可以写成一个多项式的形式……而且这个多项式如果强行套在合数上，满足积性，我也不知道有没有除了(x^{k})别的多项式惹
假如(F(x) = x^{k})吧
我们想要计算这个东西
(g(n,j))表示前(n)个数里，质数的和，加上合数中最小质因子大于(P_{j})的和
那么，怎么求呢
我们考虑已经求好(g(n,j - 1))这个数组
那么如果(P_{j}^{2} > n)
那么很妙，就是(g(n,j) = g(n,j - 1))，因为没有以(P_{j})为最小质因子的合数
如果(P_{j}^{2} leq N)
那么考虑从(g(n,j - 1))到(g(n,j))减少了什么
例如，我显然减少了合数中最小质因子为(P_{j})的值，因为定义里是大于，(g(n,j))中的合数质因子至少是(P_{j + 1})
那么最小质因子为(P_{j})的值怎么算呢
我们先要提出一个(P_{j})出来，然后在剩下的(lfloor frac{n}{P_{j}} floor)中找出所有最小质因子大于等于(P_{j})的个数
好像(g(n,j - 1))就很不错！然而其中有一些小于(P_{j})的质数，考虑把它们扣掉，这个可以预处理，因为(P_{j})的范围在(sqrt{N})以内
记(sum[j] = sum_{i = 1}^{j} F(P_{i}))
那么我们有
同时我们发现，我们只需要用到(lfloor frac{n}{i} floor)的取值，一共(2sqrt{N})个位置，而且质数的话也只需要前(sqrt{N})个
初始化的话，对于(g(n,0))，我们把这个(F(x))多项式套在每个数上面，求一个(sum_{i = 2}^{n} x^{i})这个在k固定的时候有k + 1次多项式的公式可以快速计算
一般不会太大吧，如果拉格朗日插值的话可能和(sqrt{n})相乘会有点凉
g的用处
求完(g)就有很神奇的事情发生了！如果你认为(F(x) = 1)的话，我们甚至可以求出(n)以内质数的个数！在我不知道这些质数是什么的情况下！这简直是我这个菜鸡以前想都不敢想的事
不过，我们的征途是星辰大海1到n的前缀和
那么我们类似的定义一个(S(n,j))表示前(n)个数里，所有数的最小的质因子大于等于(P_{j})的方案数（注意这里就不包括小于(P_{j})的质数了）
我们先搞定(S(n,j))里的所有质数的(F(x))的和
方法是(g(n,|P|) - sum[j - 1])
那么我们枚举每个数的最小质因子
直接列个式子就是
意义就是我们枚举一个数的最小质因子是(P_{k}^{e})，然后乘上(lfloorfrac{n}{P_{k}^{e}} floor)中最小质因子大于(P_{j + 1})的
但是所有的(F(P_{k}^{e}),egeq 2)都会漏掉，所以要重算一下
这样的话复杂度就是……我抄的那个人没写……好像是(O(frac{n^{3/4}}{log n}))
好吧，这不重要……毕竟我……也不会证复杂度
这么看起来仿佛挺简单的好难
几道例题
LOJ#6053. 简单的函数
恍然觉得WC好像见过这道题……咋就变成板子了……板子都不会我也要废了qwq
考虑质数的函数值怎么表示，由于只有2是偶数，除了2以外剩下的质数值都是(f(p) = p - 1)，然后(f(2) = 3)
这个时候我们要大胆想象，我们就认为所有质数都是(f(p) = p - 1)，遇到2和1这个难受位置特判一下就好了
但是这样不满足积性
不过一个多项式嘛，我把每一块都拆开不就满足积性了吗
设(F(x) = x)，统计一个(g(n,|P|))
再来一个(F(x) = 1)，统计一个(h(n,|P|))
可以同时进行，方法都是一样的
于是我们的(S(n,j))初始的时候统计质数要加上(g(n,|P|) - h(n,|P|) - sum_{i = 1}^{j - 1}(P_{i} - 1))
如果(j = 1)的话，我们把2算成了1，那就再加上2
于是照着刚才的讨论，问题就完美解决了
最后加上1处的函数值是1
#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('
')
#define eps 1e-10
#define ba 47
#define MAXN 200005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template<class T>
void read(T &res) {
res = 0;T f = 1;char c = getchar();
while(c < '0' || c > '9') {
if(c == '-') f = -1;
c = getchar();
}
while(c >= '0' && c <= '9') {
res = res * 10 +c - '0';
c = getchar();
}
res *= f;
}
template<class T>
void out(T x) {
if(x < 0) {x = -x;putchar('-');}
if(x >= 10) {
out(x / 10);
}
putchar('0' + x % 10);
}
const int MOD = 1000000007;
int inc(int a,int b) {
return a + b >= MOD ? a + b - MOD : a + b;
}
int mul(int a,int b) {
return 1LL * a * b % MOD;
}
int fpow(int x,int c) {
int res = 1,t = x;
while(c) {
if(c & 1) res = mul(res,t);
t = mul(t,t);
c >>= 1;
}
return res;
}
void update(int &x,int y) {
x = inc(x,y);
}
int64 N,pos[MAXN],sqr;
int prime[MAXN],tot,sum[MAXN];
bool nonprime[MAXN];
int cnt,g[MAXN],h[MAXN];
int id[2][MAXN];
int getsum(int64 x) {
int u = x % MOD;
int res = 1LL * u * (u + 1) % MOD;
res = mul(res,(MOD + 1) / 2);
return res;
}
void process(int64 S) {
for(int i = 2 ; i <= S ; ++i) {
if(!nonprime[i]) {
prime[++tot] = i;
sum[tot] = sum[tot - 1];
update(sum[tot],i);
}
for(int j = 1 ; j <= tot ; ++j) {
if(prime[j] > S / i) break;
nonprime[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
}
int S(int64 x,int y) {
if(x <= 1 || prime[y] > x) return 0;
int res = 0;
if(y == 1) res += 2;
int k = x <= sqr ? id[0][x] : id[1][N / x];
update(res,inc(g[k],MOD - inc(inc(sum[y - 1],MOD - y + 1),h[k])));
for(int j = y ; j <= tot ; ++j) {
int64 t1 = prime[j],t2 = 1LL * prime[j] * prime[j];
if(t2 > x) break;
for(int e = 1 ; t2 <= x ; t1 = t2,t2 = t2 * prime[j],++e) {
int t = inc(mul(prime[j] ^ e,S(x / t1,j + 1)),prime[j] ^ (e + 1));
update(res,t);
}
}
return res;
}
int main(){
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
read(N);
sqr = sqrt(N);
process(sqr);
for(int64 i = 1 ; i <= N ; ++i) {
int64 r = N / (N / i);
pos[++cnt] = N / i;
h[cnt] = (N / i - 1) % MOD;
g[cnt] = getsum(N / i);update(g[cnt],MOD - 1);
if(N / i <= sqr) id[0][N / i] = cnt;
else id[1][N / (N / i)] = cnt;
i = r;
}
for(int j = 1 ; j <= tot ; ++j) {
for(int i = 1 ; i <= cnt && 1LL * prime[j] * prime[j] <= pos[i] ; ++i) {
int k = (pos[i] / prime[j]) <= sqr ? id[0][pos[i] / prime[j]] : id[1][N / (pos[i] / prime[j])];
int t = inc(g[k],MOD - sum[j - 1]);
t = mul(t,prime[j]);
update(g[i],MOD - t);
update(h[i],MOD - inc(h[k],MOD - (j - 1)));
}
}
int ans = inc(S(N,1),1);
out(ans);enter;
return 0;
}
LOJ#572. 「LibreOJ Round #11」Misaka Network 与求和
当时我打的时候不会min25
一年过去了
我竟然现在才学min25？！
这个，好像是类似min25的一种求和方法
首先我们熟练的把gcd转换成互质数对
这个时候我们发现要求的就是(f * mu(T))
然而……
回忆起学过的杜教筛，我们可以卷上一个(I(x) = 1)
(f * mu * I = f * varepsilon = f)
然后套到杜教筛里
可以得到
然后我们要做的就是求(f(n))的前缀和了
用一种类似min25的方法
设(s(n,j))为([1,n])里面，质数大于等于(j - 1)的和，(f(n) = S(n,1))，默认(prime[0] = 1)
再枚举第二大的质数的时候才统计贡献
考虑min25的时候，我们分了质数，和合数两部分统计
如果是质数的部分，统计了([1,lfloor frac{n}{prime[j - 1]} floor])中，大于(prime[j - 1])的一个质数，我们用(prime[j - 1])乘上这个质数，这个数第二大的质因子即为(prime[j - 1])
如果是合数的部分，从(j)到(tot)枚举最小质数的时候，显然前面加的是已经统计完次大质因子贡献的合数，所以系数是1，唯一漏掉的是(prime[j]^{e + 1})，加上即可
列一个式子就是
unordered_map比手写哈希慢了两倍不止……感觉不行……
#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int, int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('
')
#define eps 1e-10
#define ba 47
#define MAXN 200005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template <class T>
void read(T &res) {
res = 0;
T f = 1;
char c = getchar();
while (c < '0' || c > '9') {
if (c == '-')
f = -1;
c = getchar();
}
while (c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = getchar();
}
res *= f;
}
template <class T>
void out(T x) {
if (x < 0) {
x = -x;
putchar('-');
}
if (x >= 10) {
out(x / 10);
}
putchar('0' + x % 10);
}
u32 N, sqr, K;
u32 prime[MAXN], tot, h[MAXN], pos[MAXN];
u32 id[2][MAXN], cnt, pw[MAXN], rec;
bool nonprime[MAXN];
int c[10000005];
u32 getpos(u32 x) { return x <= sqr ? id[0][x] : id[1][N / x]; }
struct HASH {
int sumE,head[974711 + 5];u32 mo = 974711;
struct node {
pair<u32,u32> x;u32 val;int next;
}E[10000005];
u32 hash_key(pair<u32,u32> t) {
return t.fi << 16 | t.se;
}
void add(pair<u32,u32> x,u32 v) {
u32 u = hash_key(x) % mo;
E[++sumE].x = x;E[sumE].val = v;
E[sumE].next = head[u];head[u] = sumE;
}
bool count(pair<u32,u32> x) {
u32 u = hash_key(x) % mo;
for(int i = head[u] ; i ; i = E[i].next) {
if(E[i].x == x) return true;
}
return false;
}
u32 query(pair<u32,u32> x) {
u32 u = hash_key(x) % mo;
for(int i = head[u] ; i ; i = E[i].next) {
if(E[i].x == x) return E[i].val;
}
return 0;
}
}hsh;
u32 rf[MAXN];
bool vis[MAXN];
u32 fpow(u32 x, u32 c) {
u32 res = 1, t = x;
while (c) {
if (c & 1)
res = res * t;
t = t * t;
c >>= 1;
}
return res;
}
void Process() {
prime[0] = 1;
pw[0] = 1;
for (u32 i = 2; i <= sqr; ++i) {
if (!nonprime[i]) {
prime[++tot] = i;
pw[tot] = fpow(i, K);
}
for (int j = 1; j <= tot; ++j) {
if (prime[j] > sqr / i)
break;
nonprime[i * prime[j]] = 1;
if (i % prime[j] == 0)
break;
}
}
for (u32 i = 1; i <= N; ++i) {
u32 r = N / (N / i);
pos[++cnt] = N / i;
h[cnt] = N / i - 1;
if (N / i <= sqr)
id[0][N / i] = cnt;
else
id[1][N / (N / i)] = cnt;
i = r;
}
for (int j = 1; j <= tot; ++j) {
for (int i = 1; i <= cnt && 1LL * prime[j] * prime[j] <= pos[i]; ++i) {
int k = getpos(pos[i] / prime[j]);
h[i] = h[i] - (h[k] - (j - 1));
}
}
}
u32 SF(u32 x, u32 y) {
if (x <= 1 || prime[y] > x)
return 0;
if (hsh.count(mp(x, y)))
return hsh.query(mp(x, y));
++rec;
u32 res = pw[y - 1] * (h[getpos(x)] - (y - 1));
for (u32 j = y; j <= tot && 1LL * prime[j] * prime[j] <= x; ++j) {
int64 t1 = prime[j], t2 = 1LL * prime[j] * prime[j];
for (int e = 1; t2 <= x; t1 = t2, t2 *= prime[j], ++e) {
res += SF(x / t1, j + 1) + pw[j];
}
}
// out(x);space;out(y);enter;
// printf("%.3lf
",(db)clock() / CLOCKS_PER_SEC);
hsh.add(mp(x,y),res);
return res;
}
u32 FF(u32 x) {
if (vis[getpos(x)])
return rf[getpos(x)];
u32 res = SF(x, 1);
for (u32 i = 2; i <= x; ++i) {
u32 r = x / (x / i);
res -= (r - i + 1) * FF(x / i);
i = r;
}
vis[getpos(x)] = 1;
rf[getpos(x)] = res;
return res;
}
int main() {
#ifdef ivorysi
freopen("f1.in", "r", stdin);
#endif
read(N);
read(K);
sqr = sqrt(N);
Process();
u32 ans = 0;
u32 pre = 0;
for (u32 T = 1; T <= N; ++T) {
u32 r = N / (N / T);
u32 nw = FF(r);
ans += (N / T) * (N / T) * (nw - pre);
pre = nw;
T = r;
}
out(ans);
enter;
// out(rec);enter;
// printf("%.3lf
",(db)clock() / CLOCKS_PER_SEC);
return 0;
}
【51nod】1847 奇怪的数学题
这个怎么跟上面那个题好像啊！
所以直接给出杜教筛的部分
em……完全一样好吧
发现(sgcd(n) = frac{n}{min_{p}(n)})，就是n除以最小质数
这个和g函数就已经很类似了……
考虑到(g(n,j))求的是([1,n])中的质数和最小质因子大于(p_{j})的(x^{k})和
然后只要枚举一个质数(p_{i} * g(lfloor frac{n}{p_{i}} floor,i - 1))就搞定了，最后再统计一遍纯质数的答案
不需要后面的那部分了
至于自然数幂和的那部分可以用乘方转斯特林数推一个公式出来
（还傻乎乎求了一遍后面的我真是zz）
#pragma GCC optimize("Ofast,unroll-loops,no-stack-protector,fast-math")
#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('
')
#define eps 1e-10
#define ba 47
#define MAXN 100005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template<class T>
void read(T &res) {
res = 0;T f = 1;char c = getchar();
while(c < '0' || c > '9') {
if(c == '-') f = -1;
c = getchar();
}
while(c >= '0' && c <= '9') {
res = res * 10 +c - '0';
c = getchar();
}
res *= f;
}
template<class T>
void out(T x) {
if(x < 0) {x = -x;putchar('-');}
if(x >= 10) {
out(x / 10);
}
putchar('0' + x % 10);
}
u32 N,sqr,K;
u32 h[MAXN],g[MAXN],f[MAXN],pos[MAXN],id[2][MAXN],sum[MAXN],cnt;
u32 prime[MAXN],tot,S[55][55],rec[55];
bool nonprime[MAXN];
u32 getpos(u32 x) {
return x <= sqr ? id[0][x] : id[1][N / x];
}
u32 getsum(u32 x) {
u32 res = 0;
for(u32 j = 0 ; j <= K ; ++j) {
if(x < j) break;
u32 t = S[K][j];
bool f = 0;
for(u32 h = 0 ; h < j + 1 ; ++h) {
u32 p = (x + 1 - h);
if(!f && p % (j + 1) == 0) {f = 1;p /= j + 1;}
t = t * p;
}
res += t;
}
return res;
}
u32 fpow(u32 x,u32 c) {
u32 res = 1,t = x;
while(c) {
if(c & 1) res = res * t;
t = t * t;
c >>= 1;
}
return res;
}
void Process() {
for(u32 i = 2 ; i <= sqr ; ++i) {
if(!nonprime[i]) {
prime[++tot] = i;
sum[tot] = sum[tot - 1] + fpow(i,K);
}
for(int j = 1 ; j <= tot ; ++j) {
if(prime[j] > sqr / i) break;
nonprime[prime[j] * i] = 1;
if(i % prime[j] == 0) break;
}
}
S[0][0] = 1;
for(int i = 1 ; i <= K ; ++i) {
for(int j = 1 ; j <= i ; ++j) {
S[i][j] = S[i - 1][j - 1] + j * S[i - 1][j];
}
}
for(u32 i = 1 ; i <= N ; ++i) {
u32 r = N / (N / i);
pos[++cnt] = N / i;
if(N / i <= sqr) id[0][N / i] = cnt;
else id[1][N / (N / i)] = cnt;
h[cnt] = N / i - 1;
g[cnt] = getsum(N / i) - 1;
i = r;
}
for(u32 j = 1 ; j <= tot ; ++j) {
for(u32 i = 1 ; i <= cnt && 1LL * prime[j] * prime[j] <= pos[i] ; ++i) {
u32 k = getpos(pos[i] / prime[j]);
g[i] -= (sum[j] - sum[j - 1]) * (g[k] - sum[j - 1]);
h[i] -= h[k] - (j - 1);
f[i] += g[k] - sum[j - 1];
}
}
for(int i = 1 ; i <= cnt ; ++i) f[i] += h[i];
}
bool vis[MAXN];
u32 rs[MAXN];
u32 UT(u32 n) {
if(vis[getpos(n)]) return rs[getpos(n)];
u32 res = f[getpos(n)];
for(u32 i = 2 ; i <= n ; ++i) {
u32 r = n / (n / i);
res -= (r - i + 1) * UT(n / i);
i = r;
}
vis[getpos(n)] = 1;rs[getpos(n)] = res;
return res;
}
void Solve() {
read(N);read(K);sqr = sqrt(N);
Process();
u32 ans = 0,pre = 0,nw;
for(u32 i = 1 ; i <= N ; ++i) {
u32 r = N / (N / i);
nw = UT(r);
ans += (N / i) * (N / i) * (nw - pre);
pre = nw;
i = r;
}
out(ans);enter;
}
int main(){
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
Solve();
return 0;
}