Prime Sum
(Math/primesum.hpp)

• View this file on GitHub
• Last update: 2023-01-16 20:41:46+09:00
• Include: #include "Math/primesum.hpp"

使い方

PrimeSum(ll n): テンプレートには

• T: 返り値の型
• T (*F)(ll): 完全乗法的関数 $f$ のprefix sumを返す関数 を指定。

T operator[](ll x): $\sum_{p \leq x:\mbox{prime}} f(p)$ を出力。 $x=\lfloor n/d \rfloor$ で表される必要がある。

Depends on

• Prime Sieve (Math/sieve.hpp)

Verified with

• Verify/LC_counting_primes.test.cpp
• Verify/YUKI_1781.test.cpp

Code

#pragma once
#include "Math/sieve.hpp"


template<typename T,T (*F)(ll)>struct PrimeSum{
    ll N,SQ;
    vector<T> lo,hi;
    PrimeSum(ll n=0):N(n),SQ(sqrtl(N)),lo(SQ+1),hi(SQ+1){
        rep(i,1,SQ+1){
            lo[i]=F(i)-1;
            hi[i]=F(N/i)-1;
        }
        auto ps=sieve(SQ);
        for(auto& p:ps){
            ll q=ll(p)*p;
            if(q>N)break;
            T sub=lo[p-1],fp=lo[p]-lo[p-1];
            ll L=min(SQ,N/q),M=SQ/p;
            rep(i,1,M+1)hi[i]-=fp*(hi[i*p]-sub);
            rep(i,M+1,L+1)hi[i]-=fp*(lo[double(N)/(i*p)]-sub);
            for(int i=SQ;i>=q;i--)lo[i]-=fp*(lo[double(i)/p]-sub);
        }
    }
    T operator[](ll x) {
        return (x<=SQ?lo[x]:hi[N/x]);
    }
};

/**
 * @brief Prime Sum
 * @docs docs/primesum.md
 */

#line 2 "Math/sieve.hpp"

template<int L=50101010>vector<int> sieve(int N){
    bitset<L> isp;
    int n,sq=ceil(sqrt(N));
    for(int z=1;z<=5;z+=4){
        for(int y=z;y<=sq;y+=6){
            for(int x=1;x<=sq and (n=4*x*x+y*y)<=N;++x){
                isp[n].flip();
            }
            for(int x=y+1;x<=sq and (n=3*x*x-y*y)<=N;x+=2){
                isp[n].flip();
            }
        }
    }
    for(int z=2;z<=4;z+=2){
        for(int y=z;y<=sq;y+=6){
            for (int x=1;x<=sq and (n=3*x*x+y*y)<=N;x+=2){
                isp[n].flip();
            }
            for(int x=y+1;x<=sq and (n=3*x*x-y*y)<=N;x+=2){
                isp[n].flip();
            }
        }
    }
    for(int y=3;y<=sq;y+=6){
        for(int z=1;z<=2;++z){
            for(int x=z;x<=sq and (n=4*x*x+y*y)<=N;x+=3){
                isp[n].flip();
            }
        }
    }
    for(int n=5;n<=sq;++n)if(isp[n]){
        for(int k=n*n;k<=N;k+=n*n){
            isp[k]=false;
        }
    }
    isp[2]=isp[3]=true;

    vector<int> ret;
    for(int i=2;i<=N;i++)if(isp[i]){
        ret.push_back(i);
    }
    return ret;
}

/**
 * @brief Prime Sieve
 */
#line 3 "Math/primesum.hpp"

template<typename T,T (*F)(ll)>struct PrimeSum{
    ll N,SQ;
    vector<T> lo,hi;
    PrimeSum(ll n=0):N(n),SQ(sqrtl(N)),lo(SQ+1),hi(SQ+1){
        rep(i,1,SQ+1){
            lo[i]=F(i)-1;
            hi[i]=F(N/i)-1;
        }
        auto ps=sieve(SQ);
        for(auto& p:ps){
            ll q=ll(p)*p;
            if(q>N)break;
            T sub=lo[p-1],fp=lo[p]-lo[p-1];
            ll L=min(SQ,N/q),M=SQ/p;
            rep(i,1,M+1)hi[i]-=fp*(hi[i*p]-sub);
            rep(i,M+1,L+1)hi[i]-=fp*(lo[double(N)/(i*p)]-sub);
            for(int i=SQ;i>=q;i--)lo[i]-=fp*(lo[double(i)/p]-sub);
        }
    }
    T operator[](ll x) {
        return (x<=SQ?lo[x]:hi[N/x]);
    }
};

/**
 * @brief Prime Sum
 * @docs docs/primesum.md
 */

Back to top page