一个筛子
一个筛子
做法:
-
求:\(\displaystyle \sum_{i=2}^{n} f(i)\),\(i\in prime\)
-
初始时:S[x]=\(\displaystyle \sum_{i=2}^{x} f(i)\)
-
for prime:primelist
- for x in range[n,sqrt(p)]
- S[x]-=(S[x/p]-S[p-1])*\(f(x)\)
- for x in range[n,sqrt(p)]
代码(HDU 6889):
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cstring>
#include <vector>
#include <map>
#include <queue>
#include <iterator>
#include <utility>
#include <algorithm>
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;
typedef pair<int, int> pii;
#define mp make_pair
#define fi first
#define se second
#define All(x) (x).begin(), (x).end()
#define Y1 "YES"
#define N1 "NO"
#define ENDL '\n'
#define count2(x) __builtin_popcount(x)
#define countleadingzero(x) __builtin_clz(x)
inline ll read() { // not solve LLONG_MIN LMAX=9,223,372,036,854,775,807
ll s = 0, w = 1;
char ch = getchar();
while (ch < '0' || ch > '9') {
if (ch == '-')
w = -1;
ch = getchar();
}
while (ch >= '0' && ch <= '9') s = s * 10 + ch - '0', ch = getchar();
return s * w;
}
typedef __uint128_t returnT;
typedef unsigned long long argT;
typedef unsigned int sqargT;
constexpr argT maxm=1e5+1;
sqargT primelist[maxm],pcnt=0;
bool isnprime[maxm];
double inv[maxm];
inline void init(){
constexpr argT m=maxm-1;
// constexpr argT mz=m/log10(m);
primelist[++pcnt]=2;
for(argT i = 1;i <= m; ++i)
inv[i] = 1.0 / i;
// i*i may overflow
for(argT j = 4; j <=m; j+=2)
isnprime[j] = true;
for (argT i = 3; i <= m; i+=2){
if(!isnprime[i]){
primelist[++pcnt]=i;
// i*i may overflow
for(argT j = i * i; j <=m; j+=i)
isnprime[j] = true;
}
}
return ;
}
argT a[maxm];
argT b[maxm];
//int ppt=0;
inline argT n34divlognsieve(const argT n) {
if (n <= 1) return 0;
const sqargT m = sqrt(n);
constexpr double eps = 1e-7;
for (sqargT i = 2; i <= m; ++i)
a[i] = a[i - 1] + i;
for (sqargT i = 1; i <= m; ++i) {
returnT sum_n = n / i;
b[i] = (sum_n + sum_n * sum_n) / 2 - 1;
}
const sqargT pt=upper_bound(primelist+1,primelist+1+pcnt,m)-primelist;
for (sqargT r=1;r<pt;++r) {
const sqargT prime=primelist[r];
const argT t = a[prime - 1];
const double invprime=inv[prime];
const sqargT limmp = m * invprime + eps;
for (sqargT i = 1; i <= limmp; ++i)
b[i] -= (b[i * prime] - t) * prime;
const sqargT limnpp = n / prime / prime;
const argT limnp = n / prime;
const sqargT lim2=min(limnpp,m);
double *ptr=&inv[limmp + 1];
for (sqargT i = limmp + 1; i <= lim2 ; ++i,++ptr)
b[i] -= (a[(unsigned int)(limnp * (*ptr) + eps)] - t) * prime;
const sqargT limpp = prime * prime;
for (sqargT i = m; i >= limpp; --i)
a[i] -= (a[(unsigned int)(i * invprime + eps)] - t) * prime;
}
return b[1];
}
int main() {
init();
int t=read();
while(t--){
argT n=read();
int k=read();
returnT ans = (returnT)n*(3+n)/2+n34divlognsieve(n+1)-4 ;
printf("%d\n",(int)(ans%k));
}
return 0;
}
耗时
使用i7-10710U,输入的n为1e10时,ans函数执行平均需要约83ms,稍作优化即可达到48ms,HDU 6889实测514MS!
代码(HDU 5901)
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
typedef unsigned long long ull;
typedef unsigned int uint;
constexpr uint maxm=sqrt(1e11)+1;//maxn=1e12
uint primelist[maxm],pcnt=0;
bool isnprime[maxm];
double inv[maxm];
inline void init(){
constexpr uint m=maxm-1;
primelist[++pcnt]=2;
for(uint i=1;i<=m;++i)inv[i]=1.0/i;
for(uint j=4;j<=m;j+=2)isnprime[j]=true;
for(uint i=3;i<=m;i+=2){
if(!isnprime[i]){
primelist[++pcnt]=i;
// i*i may overflow
for(ull j=(ull)i*i;j<=m;j+=i)
isnprime[j]=true;
}
}
return ;
}
uint a[maxm];
ull b[maxm];
inline ull n34divlognsieve(const ull n) {
if (n<=1) return 0;
const uint m=sqrt(n);
constexpr double eps=1e-7;
for(uint i=2;i<=m;++i)a[i]=a[i-1]+1;//sum
for(uint i=1;i<=m;++i)b[i]=n/i-1;//sum
const uint pt=upper_bound(primelist+1,primelist+1+pcnt,m)-primelist;
for(uint r=1;r<pt;++r){
const ull prime=primelist[r];
const uint limmp=m*inv[prime]+eps;
for (uint i=1;i<=limmp;++i)
b[i]-=(b[i*prime]-a[prime-1]);
double *ptr=&inv[limmp+1];
for (uint i=limmp+1;i<=min((uint)(n/prime/prime),m);++i,++ptr)
b[i]-=(a[(uint)(n/prime*(*ptr)+eps)]-a[prime-1]);
for (uint i=m;i>=prime*prime;--i)
a[i]-=(a[(uint)(i*inv[prime]+eps)]-a[prime-1]);
}
return b[1];
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
init();
ull n;
while(cin>>n){
cout<<n34divlognsieve(n)<<'\n';
}
return 0;
}