結果

問題 No.1661 Sum is Prime (Hard Version)
コンテスト
ユーザー monnu
提出日時 2021-08-31 19:07:33
言語 C++14
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 1,775 ms / 3,000 ms
コード長 4,182 bytes
コンパイル時間 2,904 ms
コンパイル使用メモリ 190,140 KB
実行使用メモリ 5,248 KB
最終ジャッジ日時 2024-12-30 16:59:27
合計ジャッジ時間 18,604 ms
ジャッジサーバーID
（参考情報） 		judge2 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 22
権限があれば一括ダウンロードができます

ソースコード

diff # 

#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
using namespace std;
using ll=long long;

ll counting_primes(ll x){
  if(x<=70){
    vector<int> p_table(x+1,1);
    p_table[0]=0;
    p_table[1]=0;
    for(int i=2;i*i<=x;i++){
      if(p_table[i]==0){
        continue;
      }
      for(int j=2;i*j<=x;j++){
        p_table[i*j]=0;
      }
    }
    ll ret=0;
    for(int i=1;i<=x;i++){
      if(p_table[i]){
        ret++;
      }
    }
    return ret;
  }

  //P2(x,a)を計算
  int N=1;
  while(((ll)N+1)*((ll)N+1)*((ll)N+1)<=x){
    N++;
  }
  int NN=(int)pow(x,2.0/3.0);//誤差が怖い
  vector<int> p_table(N+1,1);
  vector<int> P;
  p_table[0]=0;
  p_table[1]=0;
  for(int i=2;i*i<=N;i++){
    if(p_table[i]==1){
      for(int j=2;i*j<=N;j++){
        p_table[i*j]=0;
      }
    }
  }
  for(int i=1;i<=N;i++){
    if(p_table[i]==1){
      P.push_back(i);
    }
  }
  int n=P.size();
  int pi1=n;
  vector<int> cnt(N+1);
  int pi2=-1;
  int x_half=1;
  while(((ll)x_half+1)*((ll)x_half+1)<=x){
    x_half++;
  }
  ll S=0;
  vector<int> I(N+1);
  for(int j=2;(j-1)*N+1<=NN;j++){
    for(int i=1;i<=N;i++){
      p_table[i]=1;
    }
    for(int i=0;i<n;i++){
      int index=((j-1)*N/P[i]+1)*P[i];
      while(index<=j*N&&index<=NN){
        p_table[index-(j-1)*N]=0;
        index+=P[i];
      }
    }
    cnt[0]=pi1;
    for(int i=1;i<=N&&(j-1)*N+i<=NN;i++){
      cnt[i]=cnt[i-1]+p_table[i];
      if((j-1)*N+i==x_half){
        pi2=cnt[i];
      }
    }
    ll left_tmp=max<ll>(x/((ll)min(j*N+1,NN+1)),(ll)N);
    ll right_tmp=min<ll>(x/(ll)((j-1)*N+1),(ll)x_half);
    if(left_tmp<right_tmp){
      int left=(int)left_tmp;
      int right=(int)right_tmp;
      for(int i=left+1;i<=right;i++){
        I[i-left]=1;
      }
      for(int i=0;i<n&&P[i]*P[i]<=x_half;i++){
        int index=(left/P[i]+1)*P[i];
        while(index<=right){
          I[index-left]=0;
          index+=P[i];
        }
      }
      for(int i=left+1;i<=right;i++){
        if(I[i-left]==1){
          S+=(ll)cnt[(int)(x/(ll)i)-(j-1)*N];
        }
      }
    }
    pi1=cnt[min(NN,j*N)-(j-1)*N];
  }
  S+=(ll)n*((ll)n-1)/2;
  S-=(ll)pi2*((ll)pi2-1)/2;

  //phi(x,a)の計算
  //S1
  ll S1=0;
  vector<int> mu(N+1,1);
  mu[0]=0;
  for(int i=0;i<n;i++){
    for(int j=P[i];j<=N;j+=P[i]){
      mu[j]*=-1;
    }
  }
  for(int i=0;i<n&&P[i]*P[i]<=N;i++){
    for(int j=P[i]*P[i];j<=N;j+=P[i]*P[i]){
      mu[j]=0;
    }
  }
  for(int i=1;i<=N;i++){
    S1+=(ll)mu[i]*(x/(ll)i);
  }

  //S2
  ll S2=0;
  vector<int> f(N+1,0);
  f[1]=0;
  for(int i=0;i<n;i++){
    for(int j=P[i];j<=N;j+=P[i]){
      if(f[j]==0){
        f[j]=P[i];
      }
    }
  }
  int tmp_N=N;
  int m=0;
  for(int i=0;tmp_N>0;i++,tmp_N>>=1){
    m++;
  }

  vector<ll> phi(n+1,0);
  vector<vector<int>> a(m+1,vector<int>(N+2,0));
  for(int k=1;(k-1)*N+1<=NN;k++){
    for(int i=0;i<=m;i++){
      for(int j=1;j<=(N/(1LL<<i))+1;j++){
        a[i][j]=min(1<<i,min(k*N,NN)-(k-1)*N-(j-1)*(1<<i));
      }
    }
    for(int b=0;b<n;b++){
      int left=(int)max<ll>(0,x/(((ll)min(k*N,NN)+1)*(ll)P[b]));
      int right=(int)min<ll>((ll)N,x/((ll)((k-1)*N+1)*(ll)P[b]));
      for(int i=left+1;i<=right;i++){
        if(f[i]>P[b]&&mu[i]!=0){
          int y=(int)(x/((ll)i*(ll)P[b]));
          int l=y-(k-1)*N;
          ll sum=phi[b];
          for(int i=0;l>0;i++,l>>=1){
            if((l&1)==1){
              sum+=(ll)a[i][l];
            }
          }
          //cout<<mu[i]<<' '<<sum<<'\n';
          S2-=(ll)mu[i]*sum;
        }
      }
      int l=N;
      for(int i=0;l>0;i++,l>>=1){
        if((l&1)==1){
          phi[b]+=(ll)a[i][l];
        }
      }
      int index=((k-1)*N/P[b]+1)*P[b];
      while(index<=k*N&&index<=NN){
        if(a[0][index-(k-1)*N]==1){
          int l=index-(k-1)*N;
          for(int i=0;i<=m;i++){
            a[i][(l+(1<<i)-1)/(1<<i)]--;
          }
        }
        index+=P[b];
      }
    }
  }
  ll ret=(ll)n-S+S1+S2-1;
  //cout<<n<<" "<<S<<" "<<S1<<" "<<S2<<'\n';
  return ret;
}

int main(){
  //ll x;
  //cin>>x;
  //cout<<counting_primes(x)<<'\n';
  ll L,R;
  cin>>L>>R;
  cout<<counting_primes(R)-counting_primes(L-1)+counting_primes(2*R)-counting_primes(2*L)<<'\n';

}
0