多是一篇（抄来的）min25学习笔记
多是一篇（抄来的）min25学习笔记
一个要求不少的积性函数
咱们考虑有一个积性函数，这个函数知足能够快速计算质数处的值html
且质数能够写成一个多项式的形式……并且这个多项式若是强行套在合数上，知足积性，我也不知道有没有除了\(x^{k}\)别的多项式惹node
假如\(F(x) = x^{k}\)吧c++
咱们想要计算这个东西数组
\(g(n,j)\)表示前\(n\)个数里，质数的和，加上合数中最小质因子大于\(P_{j}\)的和函数
那么，怎么求呢oop
咱们考虑已经求好\(g(n,j - 1)\)这个数组学习
那么若是\(P_{j}^{2} > n\)spa
那么很妙，就是\(g(n,j) = g(n,j - 1)\)，由于没有以\(P_{j}\)为最小质因子的合数code
若是\(P_{j}^{2} \leq N\)htm
那么考虑从\(g(n,j - 1)\)到\(g(n,j)\)减小了什么
例如，我显然减小了合数中最小质因子为\(P_{j}\)的值，由于定义里是大于，\(g(n,j)\)中的合数质因子至少是\(P_{j + 1}\)
那么最小质因子为\(P_{j}\)的值怎么算呢
咱们先要提出一个\(P_{j}\)出来，而后在剩下的\(\lfloor \frac{n}{P_{j}} \rfloor\)中找出全部最小质因子大于等于\(P_{j}\)的个数
好像\(g(n,j - 1)\)就很不错！然而其中有一些小于\(P_{j}\)的质数，考虑把它们扣掉，这个能够预处理，由于\(P_{j}\)的范围在\(\sqrt{N}\)之内
记\(sum[j] = \sum_{i = 1}^{j} F(P_{i})\)
那么咱们有
\[ g(n,j) =\left\{\begin{matrix} g(n,j - 1) & P_{j}^{2} > n\\ g(n,j - 1) - F(P_{j})(g(\lfloor \frac{n}{P_{j}} \rfloor,j - 1) - sum(j - 1)) & P_{j}^{2} \leq n \end{matrix}\right. \]
同时咱们发现，咱们只须要用到\(\lfloor \frac{n}{i} \rfloor\)的取值，一共\(2\sqrt{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]\)
那么咱们枚举每一个数的最小质因子
直接列个式子就是
\[ S(n,j) = g(n,|P|) - sum[j - 1] + \sum_{k = j}^{P_{k}^{2} \leq n}\sum_{e = 1}^{P_{k}^{e + 1} \leq n} [F(P_{k}^{e}) S(\lfloor \frac{n}{P_{k}^{e}}\rfloor,j + 1) + F(P_{k}^{e + 1})] \]
意义就是咱们枚举一个数的最小质因子是\(P_{k}^{e}\)，而后乘上\(\lfloor\frac{n}{P_{k}^{e}}\rfloor\)中最小质因子大于\(P_{j + 1}\)的
可是全部的\(F(P_{k}^{e}),e\geq 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('\n')
#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转换成互质数对
\[ ans = \sum_{d = 1}^{n} f(d)^{k} \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{t|i,j}\mu(t) \\ \sum_{d = 1}^{n}f(d)^{k} \sum_{t = 1}^{\lfloor \frac{n}{d} \rfloor}\mu(t) \lfloor \frac{n}{dt} \rfloor^{2} \\ \sum_{T = 1}^{n} \lfloor \frac{n}{T} \rfloor^{2} \sum_{d | T}f(d)\mu(\frac{T}{d}) \]
这个时候咱们发现要求的就是\(f * \mu(T)\)
然而……
回忆起学过的杜教筛，咱们能够卷上一个\(I(x) = 1\)
\(f * \mu * I = f * \varepsilon = f\)
而后套到杜教筛里
能够获得
\[ S(n) = \sum_{i = 1}^{n}f(n) - \sum_{t = 2}^{n}S(\lfloor \frac{n}{t}\rfloor) \]
而后咱们要作的就是求\(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]} \rfloor]\)中，大于\(prime[j - 1]\)的一个质数，咱们用\(prime[j - 1]\)乘上这个质数，这个数第二大的质因子即为\(prime[j - 1]\)
若是是合数的部分，从\(j\)到\(tot\)枚举最小质数的时候，显然前面加的是已经统计完次大质因子贡献的合数，因此系数是1，惟一漏掉的是\(prime[j]^{e + 1}\)，加上便可
列一个式子就是
\[ s(n,j) = prime[j - 1]^{k} * (g(n,j) - j + 1)\\+\sum_{j = 1}^{prime[j]^{2} \leq n} \sum_{e = 1}^{prime[j]^{e + 1} \leq n}s(\lfloor \frac{n}{prime[j]^{e}}\rfloor,j + 1) + prime[j]^{k} \]
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('\n')
#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\n",(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\n",(db)clock() / CLOCKS_PER_SEC);
return 0;
}
【51nod】1847 奇怪的数学题
这个怎么跟上面那个题好像啊！
因此直接给出杜教筛的部分
\[ S(n) = \sum_{i = 1}^{n} sgcd(n)^{k} - \sum_{i = 1}^{n}S(\lfloor \frac{n}{i} \rfloor) \]
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}}\rfloor,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('\n')
#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;
}
- 1. 「学习笔记」Min25筛
- 2. 【学习笔记】Min25筛
- 3. ajax学习摘抄笔记
- 4. min25筛学习小计
- 5. Java学习手抄笔记 —— Day04
- 6. Java学习手抄笔记 —— Day01~03
- 7. Lucene学习笔记(摘抄整理版)
- 8. RHCE的第一篇学习笔记
- 9. osworkflow学习笔记-这一步步是咋连起来的呢
- 10. 一篇很是详细的 Docker 学习笔记
- 更多相关文章...
- • 您已经学习了 XML Schema，下一步学习什么呢？ - XML Schema 教程
- • 我们已经学习了 SQL，下一步学习什么呢？ - SQL 教程
- • Tomcat学习笔记（史上最全tomcat学习笔记）
- • Kotlin学习（一）基本语法
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 「学习笔记」Min25筛
- 2. 【学习笔记】Min25筛
- 3. ajax学习摘抄笔记
- 4. min25筛学习小计
- 5. Java学习手抄笔记 —— Day04
- 6. Java学习手抄笔记 —— Day01~03
- 7. Lucene学习笔记(摘抄整理版)
- 8. RHCE的第一篇学习笔记
- 9. osworkflow学习笔记-这一步步是咋连起来的呢
- 10. 一篇很是详细的 Docker 学习笔记