公式

O - 数列と素数 / Number Sequence and Prime 解説 by physics0523


制約より、 (等差数列の最大の項) \(\le 10^{12}\) が分かります。よって、素数判定をするためには \(10^6\) までの素数で試し割りすればよいことになります。しかし、この解法では実行時間制限に間に合いません。

ここで、エラトステネスの篩に似たことを考えましょう。 \(2\) 以上 \(10^6\) 以下の素数で等差数列内の要素を篩っていくことを考えます。(これらの素数は予め列挙しておきます)
素数 \(p\) の倍数を探す際、以下の要領で篩っていけばよいです。

  • まず、等差数列のうち最小の \(p\) の倍数を探す。これは以下の要領で探すことができます。
    • \(t \equiv B ( \rm{mod}\) \(A)\) かつ \(t \equiv 0 (\rm{mod}\) \(p)\) なる整数 \(t\) を探す。
      • \(Ax+py=\gcd(A,p)\) なる \(x,y\) を拡張ユークリッド互除法で求める。
      • 両辺に適切な係数 \(k\) をかけ、右辺を \((B\%A)\) にする。これができない場合等差数列に \(p\) の倍数が含まれないことが示せる。
      • このとき、とりうる \(t\) はある整数 \(m\) を用いて \(t = ( B \% A - kAx ) + \rm{lcm}\)\((A,p) \times m\) と表現できる。
    • 上記の条件を満たす \(t\) のうち、等差数列の中に含まれる最小のものを求める。それが求めたい最小の \(p\) である。
  • そこから \(\rm{lcm}\)\((A,p)/A\) 項ごとに篩っていく。

この篩を実装することによりこの問題に正解できます。 \(\rm{lcm}\)\((A,p)/A \neq p\) となる \(p\) は定数個しかないため、篩う部分の計算量は通常のエラトステネスの篩と等しいこともわかり、この問題に十分高速に正解できます。

実装例 (C++):

#include<bits/stdc++.h>

using namespace std;

long long extgcd(long long a,long long b,long long &x,long long &y){
  if(b==0){
    x=1; y=0;
    return a;
  }
  long long res=extgcd(b,a%b,y,x);
  y-=(a/b)*x;
  return res;
}

vector<long long> p_list(){
  vector<bool> fl(1e6+5,true);
  vector<long long> res;
  for(int i=2;i<1e6+5;i++){
    if(fl[i]){
      res.push_back(i);
      for(int j=2*i;j<1e6+5;j+=i){
        fl[j]=false;
      }
    }
  }
  return res;
}

long long llceil(long long a,long long b){if(a%b==0){return a/b;}return (a/b)+1;}

int main(){
  long long n,a,b;
  cin >> n >> a >> b;
  vector<long long> pl=p_list();
  vector<bool> fl(n,true);
  for(auto &p : pl){
    // t \equiv B (mod A)
    // t \equiv 0 (mod p)
    long long x,y;
    long long g=extgcd(a,p,x,y);
    long long l=(a*p)/g; // lcm
    // ax + py = g
    // kyp = kg - kxa
    // kg == B%A
    // then, t = (kg - kxa) + l * (integer)
    if((b%a)%g){continue;}
    long long k=(b%a)/g;
    // minimum positive t
    long long t=(k*((a*p+g-x*a)%l))%l;
    long long fe;
    if(b-t>0){fe=t+l*llceil(b-t,l);}
    else{fe=t;}
    long long fi=(fe-b)/a;
    for(long long i=fi;i<n;i+=(l/a)){
      if((a*i+b)==p){continue;}
      fl[i]=false;
    }
  }

  int res=0;
  for(int i=0;i<n;i++){
    if(fl[i]){res++;}
  }
  cout << res << '\n';
  return 0;
}

投稿日時:
最終更新: