あなたの範囲は 1,000 万までしかありませんが、これはこの種のものとしては小さいものです。2 つの提案があります。
1) 適切な間隔で pi(n) のテーブルを作成し、セグメント化されたエラトステネスのふるいを使用して、目的の値を囲む 2 つのテーブル エントリ間の素数を数えます。間隔のサイズによって、必要なテーブルのサイズと結果の計算速度が決まります。
2) Legendre の phi(x,a) 関数と Lehmer の素数計算式を使用して、結果を直接計算します。ファイ関数にはいくらかのストレージが必要ですが、正確な量はわかりません。
2 つのうち、問題のサイズを考えると、おそらく最初の選択肢を選択します。セグメント化されたエラトステネスの篩とLehmer の素数カウント関数の両方の実装は、私のブログで入手できます。
編集1：
振り返ってみると、3 番目の選択肢があります。
3) 対数積分を使用して pi(n) を推定します。これは単調に増加し、必要な間隔で常に pi(n) よりも大きくなります。しかし、差は小さく、200 を超えることはありません。したがって、1,000 万未満のすべての値の差を事前に計算し、200 の変化点の表を作成してから、要求に応じて対数積分を計算し、補正係数を調べることができます。テーブル。または、リーマンの R 関数で同様のことを行うこともできます。
3 番目の選択肢は最小量のスペースを必要としますが、最初の選択肢に必要なスペースはそれほど大きくないと思います。ふるい分けはおそらく対数積分を計算するよりも高速です。だから私は最初の推奨事項に固執します。私のブログには、対数積分とリーマン R 関数の両方の実装があります。
編集2：
コメントが示すように、それはうまくいきませんでした。私の 3 番目の提案は無視してください。
うまくいかない解決策を提案した私の過ちの罰として、私は pi(n) の値のテーブルとセグメント化されたエラトステネスのふるいを使用して n < 10000000 の pi(n) の値を計算するプログラムを書きました。元の投稿者が要求した C ではなく、Python を使用します。これは、Python が教育目的でより単純で読みやすいためです。
1,000 万の平方根未満のふるい素数を計算することから始めます。これらの素数は、pi(n) の値の表を作成するときと、最終的な答えを計算するふるいを実行するときの両方で使用されます。1000 万の平方根は 3162.3 です。ふるい素数として 2 を使用したくありません -- 奇数のみをふるい分け、2 を特別な場合として扱います -- しかし、次の素数を平方根よりも大きくしたいので、次のリストが素数のふるい分けが使い果たされることはありません (エラーの原因になります)。そこで、エラトステネスのふるいのこの非常に単純なバージョンを使用して、ふるい素数を計算します。
def primes(n):
b, p, ps = [True] * (n+1), 2, []
for p in xrange(2, n+1):
if b[p]:
ps.append(p)
for i in xrange(p, n+1, p):
b[i] = False
return ps
エラトステネスのふるいは 2 つの部分で機能します。最初に、2 から始めて、目標数よりも小さい数のリストを作成します。次に、最初の交差していない数から始めて、リストを繰り返し実行し、リストから数のすべての倍数を取り消します。最初は、2 が交差していない最初の数字なので、4、6、8、10 などを取り消します。次に、3 は次の交差していない数なので、6、9、12、15 などを取り消します。4 は 2 の倍数として取り消し線が引かれ、次の交差されていない数は 5 であるため、10、15、20、25 などを取り消します。交差していないすべての数が考慮されるまで続行します。交差していない数字は素数です。p のループは各数値を順番に考慮し、交差していない場合、i のループは倍数を取り消します。
このprimes関数は、2、3、5、7、11、13、...、3121、3137、3163 の 447 個の素数のリストを返します。リストから 2 を取り出し、446 個のふるい素数をグローバル ps 変数に格納します。
ps = primes(3163)[1:]
必要な主な関数は、範囲内の素数を数えます。count 関数の呼び出しごとに再割り当てされるのではなく、再利用できるように、グローバル配列に格納する sieve を使用します。
sieve = [True] * 500
このcount関数は、セグメント化されたエラトステネスのふるいを使用して、lo から hi までの範囲の素数を数えます (lo と hi は両方とも範囲に含まれます)。この関数には 4 つのforループがあります。最初のループはふるいをクリアし、最後のループは素数をカウントし、残りの 2 つは上記の単純なふるいに似た方法でふるい分けを実行します。
def count(lo, hi):
for i in xrange(500):
sieve[i] = True
for p in ps:
if p*p > hi: break
q = (lo + p + 1) / -2 % p
if lo+q+q+1 < p*p: q += p
for j in xrange(q, 500, p):
sieve[j] = False
k = 0
for i in xrange((hi - lo) // 2):
if sieve[i]: k += 1
return k
関数の中心はfor p in ps、ふるい分けを実行するループであり、各ふるい素数 p を順番に取得します。すべての素数がその時点で識別されるため、ふるい素数の 2 乗が範囲の制限よりも大きい場合、ループは終了します (平方根よりも大きい次の素数が必要な理由は、ふるい素数が存在するようにするためです)。ループを停止します)。神秘的な変数 q は、lo から hi の範囲の p の最小倍数のふるいへのオフセットです (範囲内の p の最小倍数ではなく、p の最小倍数のオフセットのインデックスであることに注意してください紛らわしいかもしれません）。このifステートメントは、完全平方の数を参照するときに q をインクリメントします。次に、j のループはふるいから p の倍数を打ちます。
countこの関数を 2 つの方法で使用します。最初の使用では、1000 の倍数で pi(n) の値のテーブルを作成します。2 番目の使用では、テーブル内で補間します。テーブルをグローバル変数 piTable に保存します。
piTable = [0] * 10000
メモリ使用量を 50 キロバイト以内に抑えるために、元の要求に基づいてパラメータ 1000 と 10000 を選択します。(はい、元の投稿者がその要件を緩和したことは知っています。しかし、とにかくそれを尊重することができます。) 1 万の 32 ビット整数は 40,000 バイトのストレージを必要とし、lo から hi までのわずか 1000 の範囲をふるいにかけるには 500 バイトしか必要としません。ストレージと非常に高速になります。他のパラメーターを試して、それらがプログラムのスペースと時間の使用にどのように影響するかを確認することをお勧めします。の構築は、関数を 1 万回piTable呼び出すことによって行われます。count
for i in xrange(1, 10000):
piTable[i] = piTable[i-1] + \
count(1000 * (i-1), 1000 * i)
この時点までのすべての計算は、実行時ではなくコンパイル時に行うことができます。私がideone.comでこれらの計算を行ったとき、約 5 秒かかりましたが、その時間はカウントされません。原則として、コードを実行時からコンパイル時に移動する機会を探して、プログラムを非常に高速に実行する必要があります。
あとは、n 以下の素数を実際に計算する関数を書くだけです。
def pi(n):
if type(n) != int and type(n) != long:
raise TypeError('must be integer')
if n < 2: return 0
if n == 2: return 1
i = n // 1000
return piTable[i] + count(1000 * i, n+1)
最初のifステートメントは型チェックを行います。2 番目のifステートメントは、ばかげた入力に対して正しい応答を返します。3 番目のifステートメントは 2 を特別に処理します。私たちのふるい分けは 1 を素数にし、2 を複合体にしますが、どちらも間違っているので、ここで修正します。次に、要求された n より小さい piTable の最大インデックスとして i が計算され、return ステートメントによって piTable からの値がテーブル値と要求された値の間の素数のカウントに追加されます。hi limit は n+1 です。そうしないと、n が素数の場合はカウントされないからです。例として、次のように言います。
print pi(6543223)
端末に 447519 という数字が表示されます。
pi関数は非常に高速です。ideone.comでは、pi(n) への 1000 のランダム呼び出しが約 0.5 秒で計算されたので、それぞれ約 0.5 ミリ秒です。これには、素数を生成して結果を合計する時間が含まれているため、実際に pi 関数を計算する時間は 0.5 ミリ秒未満です。これは、テーブルを構築するための投資に対するかなりの見返りです。
素数を使ったプログラミングに興味がある場合は、ブログでかなりの作業を行いました。是非ご来店ください。
ここでは、ヒューリスティック ハイブリッド モデルが機能することをお勧めします。n 番目ごとの素数を格納し、素数性テストを介して線形検索を実行します。高速化するために、高速で単純な素数性テスト ( を使用したフェルマー テストなどa==2) を使用して、偽陽性を事前に計算できます。入力の最大サイズとストレージの制約に基づいた微調整は簡単に行うことができます。
入力が素数であることがアプリオリにわかっている場合は、結果を丸めただけでは正しい値を取得できない場合に、素数の修正の小さなテーブルを使用して、近似 pi(n) ≈ n / log n を使用できる場合があります。 n. 遅い力ずくのアプローチを除けば、これがサイズの制約内での最善の策だと思います。
動作するコードを次に示します。試行割り算に基づく素数性テストを、入力範囲で機能する決定論的ミラーラビンテストに置き換える必要があります。適切な狭い範囲で素数を見つけるためにふるいにかけることは試行割り算よりもうまくいくでしょうが、それは間違った方向への一歩です。
#include <stdio.h>
#include <bitset>
using namespace std;
short smallprimes[549]; // about 1100 bytes
char in[19531]; // almost 20k
// Replace me with Miller-Rabin using 2, 7, and 61.
int isprime(int j) {
if (j<3) return j==2;
for (int i = 0; i < 549; i++) {
int p = smallprimes[i];
if (p*p > j) break;
if (!(j%p)) return 0;
}
return 1;
}
void init() {
bitset<4000> siv;
for (int i = 2; i < 64; i++) if (!siv[i])
for (int j = i+i; j < 4000; j+=i) siv[j] = 1;
int k = 0;
for (int i = 3; i < 4000; i+=2) if (!siv[i]) {
smallprimes[k++] = i;
}
for (int a0 = 0; a0 < 10000000; a0 += 512) {
in[a0/512] = !a0;
for (int j = a0+1; j < a0+512; j+=2)
in[a0/512] += isprime(j);
}
}
int whichprime(int k) {
if (k==2) return 1;
int a = k/512;
int ans = 1 + !a;
for (int i = 0; i < a; i++) ans += in[i];
for (int i = a*512+1; i<k; i+=2) ans += isprime(i);
return ans;
}
int main() {
int k;
init();
while (1 == scanf("%i", &k)) printf("%i\n", whichprime(k));
}
以下は、あなたが探しているもののように聞こえます。http://www.geekviewpoint.com/java/numbers/index_of_prime。そこにコードと単体テストがあります。あなたのリストは比較的小さいので（つまり10^7）、それを処理する必要があります。
基本的に、 と の間の2すべての素数を見つけてから、インデックスを見つけるnよりも少ないすべての素数を数えます。nまた、 がn素数でない場合、関数は を返します-1。
私はまさにこれを一度やった。n を指定すると、n 番目の素数 (n = 203542528 まで) をすばやく見つけることができるコードを作成しました。つまり、約 2e8 です。または、任意の数 n について逆方向に進むことができ、n より小さい素数がいくつあるかを知ることができます。
データベースが採用されています。すべての素数を特定のポイント (上限の sqrt) まで保存します。あなたの場合、これはすべての素数を sqrt(1e7) まで保存することを意味します。それらの数は 446 あり、その時点までの最大差は 34 しかないため、そのリストを圧縮形式で保存できます。その時点を超えると、すべての k 番目の素数を保存します (k のある値について)。その後、簡単なふるいで十分です。短い間隔ですべての素数を生成します。
したがって、MATLAB で 1e7 番目の素数を見つけるには、次のようにします。
nthprime(1e7)
ans =
179424673
または、1e7 未満の素数の数を見つけることができます。
nthprime(1e7,1)
ans =
664579
ポイントは、そのようなデータベースは構築と検索が簡単だということです。データベースが 50k 以下であれば問題ありません。
あなたが提案するものが最善です。10^7 未満の素数のリストを事前に計算 (またはダウンロード) してから、それらを二分探索します。
10^7 未満の素数は 664579 個しかないため、リストは最大 2.7 MB のスペースを消費します。各インスタンスを解決するための二分探索は非常に高速で、操作は 20 回までです。