一个筛子

一个筛子

做法：

• 求：(displaystyle sum_{i=2}^{n} f(i)),(iin 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))

代码(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 '
'
#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
",(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)<<'
';
    }
    return 0;
}