线性筛求约数个数以及约数和

RT 考试的时候口胡出来的,正确性不会证,不过貌似100000内的数都是对的(现在已经会了,在此鸣谢gyz大佬)
代码:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=5e7+10,maxm=3000000;
int a[maxn],prime[maxm],f[maxn];
bool is_not_prime[maxn];
int n,mod,tot;
ll ans;
void Solve(){
	scanf("%d%d",&n,&mod);
	is_not_prime[0]=is_not_prime[1]=f[1]=1;
	for(register int i=2;i<=n;++i){
		if(!is_not_prime[i]){
			prime[++tot]=i;
			f[i]=2;
		}
		for(register int j=1;j<=tot&&prime[j]*i<=n;++j){
			is_not_prime[i*prime[j]]=1;
			if(i%prime[j]!=0) f[i*prime[j]]=f[i]*f[prime[j]];
			else{
                              f[i*prime[j]]=f[i]*2-f[i/prime[j]];
                              break;
			}
		}
	}
	
}
int main(){
	freopen("in","r",stdin);
//	freopen("out2","w",stdout);
	Solve();
	return 0;
}

关于正确性证明(以下证明来自gyz大佬):
首先根据线性筛原理,当 (prime[j]) 可以整除 (i) 时,(prime[j])(i) 的最小质因子。
为了方便,用(p)代替(prime[j])
(i=p^k imes p_1^{k_1} imes p_2^{k_2} imes)……( imes p_n^{k_n})
根据唯一分解定理,有(d[i]=(k+1) imes (k_1+1) imes (k_2+1) imes)……( imes (k_n+1))
不妨设(t=(k_1+1) imes (k_2+1) imes)……( imes (k_n+1))
则有:
(d[i]=(k+1) imes t)
(d[i/p]=k imes t)
(d[i imes p]=(k+2) imes t)
所以
(d[i imes p]=2 imes (k+1) imes t-k imes t=d[i] imes 2-d[i/p])
和传统写法相比,少开一个数组记录每个数的最小质因子出现的次数。同时也省去一些对这个数组的维护。
oiwiki上的筛法:

void pre() {
  d[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!v[i]) v[i] = 1, p[++tot] = i, d[i] = 2, num[i] = 1;
    for (int j = 1; j <= tot && i <= n / p[j]; ++j) {
      v[p[j] * i] = 1;
      if (i % p[j] == 0) {
        num[i * p[j]] = num[i] + 1;
        d[i * p[j]] = d[i] / num[i * p[j]] * (num[i * p[j]] + 1);
        break;
      } else {
        num[i * p[j]] = 1;
        d[i * p[j]] = d[i] * 2;
      }
    }
  }
}

所以这份代码里面的num数组是可以省去的

线性筛求约数和:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=5e7+10,maxm=3000000;
int a[maxn],prime[maxm];
ll num[maxn],f[maxn];
bool is_not_prime[maxn];
int n,mod,tot;
ll ans;
void Solve(){
	scanf("%d%d",&n,&mod);
	is_not_prime[0]=is_not_prime[1]=f[1]=1;
	for(register int i=2;i<=n;++i){
		if(!is_not_prime[i]){
			prime[++tot]=i;
			f[i]=num[i]=i+1;
		}
		for(register int j=1;j<=tot&&prime[j]*i<=n;++j){
			is_not_prime[i*prime[j]]=1;
			if(i%prime[j]!=0){
				f[i*prime[j]]=f[i]*f[prime[j]];
				num[i*prime[j]]=prime[j]+1;
			}
			else{
				f[i*prime[j]]=f[i]/num[i]*(num[i]*prime[j]+1);
               	                num[i*prime[j]]*=prime[j];
               	                num[i*prime[j]]++;
               	                break;
			}
		}
	}
	
}
int main(){
	freopen("in","r",stdin);
//	freopen("out2","w",stdout);
	Solve();
	return 0;
}

(写的貌似没错吧。。。。没调试)
证明见这位神犇的博客
暂时没想出来去掉这份代码里面num数组的方法

原文地址:https://www.cnblogs.com/wwcdcpyscc/p/13662918.html