求十亿内所有素数之和与个数

求十亿内所有质数的和,怎么做最快? - PlanarG的回答 - 知乎
Min25筛实现,大约每100ms算一个


const int N = 1e6 + 10;//对应素数精度可以到1e10

typedef __int128 LL;//防止1e10 * 1e10爆long long

namespace Min25
{

    int prime[N], id1[N], id2[N], flag[N], ncnt, m;

    LL g[N], sum[N], a[N], T, n;

    inline int ID(LL x)
    {
        return x <= T ? id1[x] : id2[n / x];
    }

    inline LL calc(LL x)
    {
        return x * (x + 1) / 2 - 1;
    }

    inline LL f(LL x)
    {
        return x;
    }

    inline void init()
    {
        ncnt = m = T = 0;

        T = sqrt(n + 0.5);
        for (int i = 2; i <= T; i++)
        {
            if (!flag[i])
                prime[++ncnt] = i, sum[ncnt] = sum[ncnt - 1] + i;
            for (int j = 1; j <= ncnt && i * prime[j] <= T; j++)
            {
                flag[i * prime[j]] = 1;
                if (i % prime[j] == 0)
                    break;
            }
        }
        for (LL l = 1; l <= n; l = n / (n / l) + 1)
        {
            a[++m] = n / l;
            if (a[m] <= T)
                id1[a[m]] = m;
            else
                id2[n / a[m]] = m;
            g[m] = calc(a[m]);
        }
        for (int i = 1; i <= ncnt; i++)
            for (int j = 1; j <= m && (LL)prime[i] * prime[i] <= a[j]; j++)
                g[j] = g[j] - (LL)prime[i] * (g[ID(a[j] / prime[i])] - sum[i - 1]);
    }

    inline LL solve(LL x)
    {
        if (x <= 1)
            return x;
        n = x;
        init();
        return g[ID(n)];
    }

} // namespace Min25

int main() {
    LL n; scanf("%lld", &n);
    printf("%lld\n", Min25::solve(n));
}

Meissel Lehmer Algorithm 求前n个数中素数个数 【模板】

#include <bits/stdtr1c++.h>
  
#define MAXN 100
#define MAXM 10001
#define MAXP 40000
#define MAX 400000
#define clr(ar) memset(ar, 0, sizeof(ar))
#define read() freopen("lol.txt", "r", stdin)
#define dbg(x) cout << #x << " = " << x << endl
#define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))
#define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))
#define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))
  
using namespace std;
  
namespace pcf{
    long long dp[MAXN][MAXM];
    unsigned int ar[(MAX >> 6) + 5] = {0};
    int len = 0, primes[MAXP], counter[MAX];
  
    void Sieve(){
        setbit(ar, 0), setbit(ar, 1);
        for (int i = 3; (i * i) < MAX; i++, i++){
            if (!chkbit(ar, i)){
                int k = i << 1;
                for (int j = (i * i); j < MAX; j += k) setbit(ar, j);
            }
        }
  
        for (int i = 1; i < MAX; i++){
            counter[i] = counter[i - 1];
            if (isprime(i)) primes[len++] = i, counter[i]++;
        }
    }
  
    void init(){
        Sieve();
        for (int n = 0; n < MAXN; n++){
            for (int m = 0; m < MAXM; m++){
                if (!n) dp[n][m] = m;
                else dp[n][m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];
            }
        }
    }
  
    long long phi(long long m, int n){
        if (n == 0) return m;
        if (primes[n - 1] >= m) return 1;
        if (m < MAXM && n < MAXN) return dp[n][m];
        return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);
    }
  
    long long Lehmer(long long m){
        if (m < MAX) return counter[m];
  
        long long w, res = 0;
        int i, a, s, c, x, y;
        s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);
        a = counter[y], res = phi(m, a) + a - 1;
        for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;
        return res;
    }
}
  
long long solve(long long n){
    int i, j, k, l;
    long long x, y, res = 0;
  
    for (i = 0; i < pcf::len; i++){
        x = pcf::primes[i], y = n / x;
        if ((x * x) > n) break;
        res += (pcf::Lehmer(y) - pcf::Lehmer(x));
    }
  
    for (i = 0; i < pcf::len; i++){
        x = pcf::primes[i];
        if ((x * x * x) > n) break;
        res++;
    }
  
    return res;
}
  
int main(){
    pcf::init();
    long long n, res;
  
    while (scanf("%lld", &n) != EOF){
        //res = solve(n);
        printf("%lld\n",pcf::Lehmer(n));
        //printf("%lld\n", res);
    }
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值