Min_25 筛这个东西，彻底理解花了我很长的时间，因此写点东西来记录一些本身的理解。php
它能作什么
对于某个数论函数 f，若是知足如下几个条件，那么它就能够用 Min_25 筛来快速求出这个函数的前缀和。ios
1.它是一个积性函数c++
2.对于一个质数 p ，f(p) 的表达式必须是一个项数比较小的多项式。即 
数组
3.对于一个质数 p ，f(pk的表达式必须能够由 f(p) 快速获得。微信
看起来这些条件比较紧，其实仔细想一想许多函数仍是知足这些性质的。接下来讨论它的实现过程。函数
实现过程
Min_25 筛的核心思想就是对于 ≤n 的数中，除掉质数以后，每一个数的最小质因子大小不超过 √n 。即相似埃氏筛的方式，从每一个质因子筛掉一部分的贡献。整个过程能够分红三部分，即质数的部分的贡献之和与合数的部分的贡献和，以及 1 的特判。spa
记号约定
这里用 P 表示质数集合，Pi 表示第 i 小的质数。.net
质数部分
首先咱们不妨设 f(p)=pk 。对于更加通常的项数小的多项式形式，咱们能够分别记每个项的贡献，而后分别相加。注意到这里 f(p) 是一个彻底积性函数。对于 p 为合数的状况，为了方便，咱们先姑且一样定义 f(p)=pk ，到最后你会明白这样定义的缘由。3d
咱们再设g(n,i) 表示 ≤n 的全部数 x 中，x 的最小质因子 >Pi 或者 x 自己为质数的全部 f(x) 之和，这样 g(n,|P|) 就是全部质数的贡献之和了。而 g(n,0) 的答案也十分显然，即 code
。咱们只须要考虑，如何从 g(∗,i−1) 如何获得g(∗,i) 。
考虑咱们从i−1 推到 i 的过程当中，会减小哪一个部分的贡献，容易发现最小质因子刚好为 Pi 的贡献被减掉了。若是 Pi2>n ，那么这一部分是没有贡献的，即 g(n,i)=g(n,i−1)。
对于 Pi2≤n 的状况，咱们须要把全部 Pi 的倍数，且知足其最小质因子除去 Pi 以外 >Pi−1 的贡献所有减掉。若是直接减去
你会发现贡献是减多了的，缘由就是 g 函数的值还包含了全部质数的贡献之和，仔细分析一下，能够发现，对于
的贡献在 j 这个地方已经减过了，因此还要加回来，即
式子汇总起来也就是
回头来，咱们再来看最初咱们提出的一些问题。
咱们最开始为何要定义，对于一个合数 p ，f(p)=pk 呢？无非就是方便计算，方便利用彻底积性函数的性质。你发现到最后，合数的贡献仍是没有算进去，因此对于合数的定义是无所谓的。
咱们到底在哪里利用了 f(p) 是一个彻底积性函数的性质呢？这个性质是在上式中第二种状况中利用到的。注意到咱们在减贡献的时候，并无枚举 Pi 的次数，这也就意味着 中可能包含了某些数仍然是 Pi 的倍数。因此这一步须要有一个对于全部质数 p, 的条件。因为以前咱们假定了 f(p) 对于 p 为合数的状况也知足 f(p)=pk ，所以这个条件是获得知足了的。
递归实现很显然。考虑一下非递归的实现。首先对于这个步骤，咱们有用的 nn 只有全部不一样的 ⌊ni⌋,(1≤i≤n)⌊ni⌋,(1≤i≤n)，因此能够离散掉这些值，用一个大小为 2n−−√2n 的数组存下来。因为这个递推是一层一层推的，因此实现上能够不用开二维数组，能够直接离散以后用一个一维数组存下来，须要注意枚举顺序。
合数部分
考虑咱们在求质数部分的时候，除了质数的答案之和以外，还求出了些什么东西。事实上，对于全部的 ⌊n/i⌋,(1≤i≤n) ，咱们求出了 g([n/i⌋,|P|) 的值。
咱们设 S(n,i) 表示 ≤n 的全部数 x 中，x 的最小质因子≥Pi 的 f(x) 之和。注意这里的 f(x) 在 x 为合数的时候再也不沿用质数部分时的定义，也注意这里的定义与 g 的区别。
接下来咱们先算上全部知足条件的质数贡献之和，即。对于合数，咱们直接枚举其最小质因子以及质因子的个数，直接递归计算。注意在这里形如pk,(p∈P) 的贡献并无算进去，因此还要单独加一下。式子的形式以下：
直接对着式子看，应该仍是很好理解的。对这个式子，直接递归暴力计算便可，跑得会比第一部分快。
至此 Min_25 筛的全过程就讲完了。
附Min25筛求质数和代码
#include <math.h>#include <stdio.h>#include <assert.h>#define LINT long long
inline LINT V2IDX(LINT v, LINT N, LINT Ndr, LINT nv) { return v >= Ndr ? (N/v - 1) : (nv - v);}
LINT primesum(LINT N) { LINT *S; LINT *V;
LINT r = (LINT)sqrt(N); LINT Ndr = N/r;
assert(r*r <= N and (r+1)*(r+1) > N);
LINT nv = r + Ndr - 1;
V = new LINT[nv]; S = new LINT[nv];
for (LINT i=0; i<r; i++) { V[i] = N/(i+1); } for (LINT i=r; i<nv; i++) { V[i] = V[i-1] - 1; }
for (LINT i=0; i<nv; i++) { S[i] = V[i] * (V[i] + 1) / 2 - 1; }
for (LINT p=2; p<=r; p++) { if (S[nv-p] > S[nv-p+1]) { LINT sp = S[nv-p+1]; LINT p2 = p*p; for (LINT i=0; i<nv; i++) { if (V[i] >= p2) { S[i] -= p * (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp); } else { break; } } } }
return S[0];}
int main() { LINT N = 1000000000; printf("%lld\n", primesum(N));}
附加一道刚作过的题
http://acm.hdu.edu.cn/contests/contest_showproblem.php?pid=1002&cid=909
题目代码以下：
#include <iostream>#include <string.h>#include <stdlib.h>#include <stdio.h>#include <math.h>#include <assert.h>#include <bits/stdc++.h>using namespace std;const int INF = 0x3f3f3f3f;const int N = 1000010;typedef long long ll;const int mod =1e9+7;ll n,k;
long long ksm(long long x,long long y,long long mod){ long long ans=1; while(y){ if(y&1) ans=ans*x%mod; y>>=1; x=x*x%mod; } return ans;}
bool ip(ll x){ for(int i=2;i<sqrt(x);i++) { if(x%i==0) return false; } return true;}
ll zsh(ll n){ ll ans=0; for(int i=2;i<=n;i++) { if(ip(i)) ans+=i; } return ans;}
inline ll jf(ll x, ll y, ll p, ll q) { return x >= p ? (y/x - 1) : (q - x);}
ll solution(ll N) { ll *S; ll *V; ll r = (ll)sqrt(N); ll qq = N/r; assert(r*r <= N and (r+1)*(r+1) > N); ll maxn = r + qq - 1; V = new ll[maxn]; S = new ll[maxn]; for (ll i=0; i<r; i++) V[i] = N/(i+1); for (ll i=r; i<maxn; i++) V[i] = V[i-1] - 1; for (ll i=0; i<maxn; i++) S[i] = V[i] * (V[i] + 1) / 2 - 1; for (ll p=2; p<=r; p++) { if (S[maxn-p] > S[maxn-p+1]) { ll sp = S[maxn-p+1]; ll p2 = p*p; for (ll i=0; i<maxn; i++) { if (V[i] >= p2) { S[i] -= p * (S[jf(V[i]/p, N, qq, maxn)] - sp); } else { break; } } } } return S[0];}
int main(){ int t; cin>>t; while(t--) { scanf("%lld%lld",&n,&k); ll sum=0; cout<<(solution(n+1)%k-4+(n+3)%k*n%k*ksm(2,k-2,k)%k+k)%k<<endl; } return 0;}
本文分享自微信公众号 - WHICH工做室（which_cn）。
若有侵权，请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”，欢迎正在阅读的你也加入，一块儿分享。