Amicable Pairs を見つけようとするプログラムを作成しました。これには、数の適切な約数の合計を見つける必要があります。
これが私の現在のsumOfDivisors()方法です：
int sumOfDivisors(int n)
{
int sum = 1;
int bound = (int) sqrt(n);
for(int i = 2; i <= 1 + bound; i++)
{
if (n % i == 0)
sum = sum + i + n / i;
}
return sum;
}
そのため、多くの因数分解を行う必要があり、それが私のアプリケーションの本当のボトルネックになり始めています。MAPLE に膨大な数を入力したところ、非常に高速に因数分解されました。
より高速な因数分解アルゴリズムの 1 つを教えてください。
メソッドは機能しますが、遅くなります。「あなたの数字はどれくらいですか？」使用する方法を決定します。
タイトル (および最後の行) の質問は、質問の実際の本文とはほとんど関係がないようです。友好的なペアを見つけようとしている場合、または多くの数の除数の合計を計算している場合、各数を個別に因数分解することは (可能な限り最速のアルゴリズムを使用しても) 絶対に非効率的な方法です。
除数の和関数, σ(n) = (sum of divisors of n), は乗法関数: 互いに素な m と n に対して、 がσ(mn) = σ(m)σ(n)あるので、
σ(p 1 k 1 …p r k r ) = [(p 1 k 1 +1 -1)/(p 1 -1)]…[(p r k r +1 -1)/(p r -1 )]。
したがって、任意の単純なふるい (たとえば、エラトステネスのふるいの拡張版) を使用して までの素数を見つけn、その過程で n までのすべての数を因数分解します。(たとえば、ふるいを行うときに、各 n の最小の素因数を保存します。その後n、繰り返して任意の数を因数分解できます。)これは、個別の因数分解アルゴリズムを数回使用するよりも (全体的に) 高速です。
ところで: 友好的なペアのいくつかの既知のリストが既に存在します (たとえば、こことMathWorldのリンクを参照してください)。記録を拡張しようとしていますか?
Shor のアルゴリズム: http://en.wikipedia.org/wiki/Shor%27s_algorithm
もちろん、量子コンピューターが必要ですが:D
Maple で使用されているのと同じアルゴリズムであるQuadratic Sieveから始めることをお勧めします。
これらのアルゴリズムの問題は、数値計算における多くの理論を実際に暗示していることです..
これは、Maple での整数因数分解の論文です。
「いくつかの非常に単純な命令から始めて — 「Maple で整数因数分解を高速化する」 — Maple と C の組み合わせで Quadratic Sieve 因数分解アルゴリズムを実装しました...」
1GB メモリ用の 2015 C++ バージョン 2 27ルックアップ テーブルの実装:
#include <iostream.h> // cerr, cout, and NULL
#include <string.h> // memcpy()
#define uint unsigned __int32
uint *factors;
const uint MAX_F=134217728; // 2^27
void buildFactors(){
factors=new (nothrow) uint [(MAX_F+1)*2]; // 4 * 2 * 2^27 = 2^30 = 1GB
if(factors==NULL)return; // not able to allocate enough free memory
int i;
for(i=0;i<(MAX_F+1)*2;i++)factors[i]=0;
//Sieve of Eratosthenese
factors[1*2]=1;
factors[1*2+1]=1;
for(i=2;i*i<=MAX_F;i++){
for(;factors[i*2] && i*i<=MAX_F;i++);
factors[i*2]=1;
factors[i*2+1]=i;
for(int j=2;i*j<=MAX_F;j++){
factors[i*j*2]=i;
factors[i*j*2+1]=j;
}
}
for(;i<=MAX_F;i++){
for(;i<=MAX_F && factors[i*2];i++);
if(i>MAX_F)return;
factors[i*2]=1;
factors[i*2+1]=i;
}
}
uint * factor(uint x, int &factorCount){
if(x > MAX_F){factorCount=-1;return NULL;}
uint tmp[70], at=x; int i=0;
while(factors[at*2]>1){
tmp[i++]=factors[at*2];
cout<<"at:"<<at<<" tmp:"<<tmp[i-1]<<endl;
at=factors[at*2+1];
}
if(i==0){
cout<<"at:"<<x<<" tmp:1"<<endl;
tmp[i++]=1;
tmp[i++]=x;
}else{
cout<<"at:"<<at<<" tmp:1"<<endl;
tmp[i++]=at;
}
factorCount=i;
uint *ret=new (nothrow) uint [factorCount];
if(ret!=NULL)
memcpy(ret, tmp, sizeof(uint)*factorCount);
return ret;
}
void main(){
cout<<"Loading factors lookup table"<<endl;
buildFactors(); if(factors==NULL){cerr<<"Need 1GB block of free memory"<<endl;return;}
int size;
uint x=30030;
cout<<"\nFactoring: "<<x<<endl;
uint *f=factor(x,size);
if(size<0){cerr<<x<<" is too big to factor. Choose a number between 1 and "<<MAX_F<<endl;return;}
else if(f==NULL){cerr<<"ran out of memory trying to factor "<<x<<endl;return;}
cout<<"\nThe factors of: "<<x<<" {"<<f[0];
for(int i=1;i<size;i++)
cout<<", "<<f[i];
cout<<"}"<<endl;
delete [] f;
x=30637;
cout<<"\nFactoring: "<<x<<endl;
f=factor(x,size);
cout<<"\nThe factors of: "<<x<<" {"<<f[0];
for(int i=1;i<size;i++)
cout<<", "<<f[i];
cout<<"}"<<endl;
delete [] f;
delete [] factors;
}
更新: または、単純さを犠牲にして、2 28を少し超えた範囲でもう少し範囲を広げます。
#include <iostream.h> // cerr, cout, and NULL
#include <string.h> // memcpy(), memset()
//#define dbg(A) A
#ifndef dbg
#define dbg(A)
#endif
#define uint unsigned __int32
#define uint8 unsigned __int8
#define uint16 unsigned __int16
uint * factors;
uint8 *factors08;
uint16 *factors16;
uint *factors32;
const uint LIMIT_16 = 514; // First 16-bit factor, 514 = 2*257
const uint LIMIT_32 = 131074;// First 32-bit factor, 131074 = 2*65537
const uint MAX_FACTOR = 268501119;
//const uint64 LIMIT_64 = 8,589,934,594; // First 64-bit factor, 2^33+1
const uint TABLE_SIZE = 268435456; // 2^28 => 4 * 2^28 = 2^30 = 1GB 32-bit table
const uint o08=1, o16=257 ,o32=65665; //o64=4294934465
// TableSize = 2^37 => 8 * 2^37 = 2^40 1TB 64-bit table
// => MaxFactor = 141,733,953,600
/* Layout of factors[] array
* Indicies(32-bit) i Value Size AFactorOf(i)
* ---------------- ------ ---------- ----------------
* factors[0..128] [1..513] 8-bit factors08[i-o08]
* factors[129..65408] [514..131073] 16-bit factors16[i-o16]
* factors[65409..268435455] [131074..268501119] 32-bit factors32[i-o32]
*
* Note: stopping at i*i causes AFactorOf(i) to not always be LargestFactor(i)
*/
void buildFactors(){
dbg(cout<<"Allocating RAM"<<endl;)
factors=new (nothrow) uint [TABLE_SIZE]; // 4 * 2^28 = 2^30 = 1GB
if(factors==NULL)return; // not able to allocate enough free memory
uint i,j;
factors08 = (uint8 *)factors;
factors16 = (uint16 *)factors;
factors32 = factors;
dbg(cout<<"Zeroing RAM"<<endl;)
memset(factors,0,sizeof(uint)*TABLE_SIZE);
//for(i=0;i<TABLE_SIZE;i++)factors[i]=0;
//Sieve of Eratosthenese
//8-bit values
dbg(cout<<"Setting: 8-Bit Values"<<endl;)
factors08[1-o08]=1;
for(i=2;i*i<LIMIT_16;i++){
for(;factors08[i-o08] && i*i<LIMIT_16;i++);
dbg(cout<<"Filtering: "<<i<<endl;)
factors08[i-o08]=1;
for(j=2;i*j<LIMIT_16;j++)factors08[i*j-o08]=i;
for(;i*j<LIMIT_32;j++)factors16[i*j-o16]=i;
for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i;
}
for(;i<LIMIT_16;i++){
for(;i<LIMIT_16 && factors08[i-o08];i++);
dbg(cout<<"Filtering: "<<i<<endl;)
if(i<LIMIT_16){
factors08[i-o08]=1;
j=LIMIT_16/i+(LIMIT_16%i>0);
for(;i*j<LIMIT_32;j++)factors16[i*j-o16]=i;
for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i;
}
}i--;
dbg(cout<<"Setting: 16-Bit Values"<<endl;)
//16-bit values
for(;i*i<LIMIT_32;i++){
for(;factors16[i-o16] && i*i<LIMIT_32;i++);
factors16[i-o16]=1;
for(j=2;i*j<LIMIT_32;j++)factors16[i*j-o16]=i;
for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i;
}
for(;i<LIMIT_32;i++){
for(;i<LIMIT_32 && factors16[i-o16];i++);
if(i<LIMIT_32){
factors16[i-o16]=1;
j=LIMIT_32/i+(LIMIT_32%i>0);
for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i;
}
}i--;
dbg(cout<<"Setting: 32-Bit Values"<<endl;)
//32-bit values
for(;i*i<=MAX_FACTOR;i++){
for(;factors32[i-o32] && i*i<=MAX_FACTOR;i++);
factors32[i-o32]=1;
for(j=2;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i;
}
for(;i<=MAX_FACTOR;i++){
for(;i<=MAX_FACTOR && factors32[i-o32];i++);
if(i>MAX_FACTOR)return;
factors32[i-o32]=1;
}
}
uint * factor(uint x, int &factorCount){
if(x > MAX_FACTOR){factorCount=-1;return NULL;}
uint tmp[70], at=x; int i=0;
while(at>=LIMIT_32 && factors32[at-o32]>1){
tmp[i++]=factors32[at-o32];
dbg(cout<<"at32:"<<at<<" tmp:"<<tmp[i-1]<<endl;)
at/=tmp[i-1];
}
if(at<LIMIT_32){
while(at>=LIMIT_16 && factors16[at-o16]>1){
tmp[i++]=factors16[at-o16];
dbg(cout<<"at16:"<<at<<" tmp:"<<tmp[i-1]<<endl;)
at/=tmp[i-1];
}
if(at<LIMIT_16){
while(factors08[at-o08]>1){
tmp[i++]=factors08[at-o08];
dbg(cout<<"at08:"<<at<<" tmp:"<<tmp[i-1]<<endl;)
at/=tmp[i-1];
}
}
}
if(i==0){
dbg(cout<<"at:"<<x<<" tmp:1"<<endl;)
tmp[i++]=1;
tmp[i++]=x;
}else{
dbg(cout<<"at:"<<at<<" tmp:1"<<endl;)
tmp[i++]=at;
}
factorCount=i;
uint *ret=new (nothrow) uint [factorCount];
if(ret!=NULL)
memcpy(ret, tmp, sizeof(uint)*factorCount);
return ret;
}
uint AFactorOf(uint x){
if(x > MAX_FACTOR)return -1;
if(x < LIMIT_16) return factors08[x-o08];
if(x < LIMIT_32) return factors16[x-o16];
return factors32[x-o32];
}
void main(){
cout<<"Loading factors lookup table"<<endl;
buildFactors(); if(factors==NULL){cerr<<"Need 1GB block of free memory"<<endl;return;}
int size;
uint x=13855127;//25255230;//30030;
cout<<"\nFactoring: "<<x<<endl;
uint *f=factor(x,size);
if(size<0){cerr<<x<<" is too big to factor. Choose a number between 1 and "<<MAX_FACTOR<<endl;return;}
else if(f==NULL){cerr<<"ran out of memory trying to factor "<<x<<endl;return;}
cout<<"\nThe factors of: "<<x<<" {"<<f[0];
for(int i=1;i<size;i++)
cout<<", "<<f[i];
cout<<"}"<<endl;
delete [] f;
x=30637;
cout<<"\nFactoring: "<<x<<endl;
f=factor(x,size);
cout<<"\nThe factors of: "<<x<<" {"<<f[0];
for(int i=1;i<size;i++)
cout<<", "<<f[i];
cout<<"}"<<endl;
delete [] f;
delete [] factors;
}
数字がどれだけ大きいかによって異なります。友好的なペアを探している場合、多くの因数分解を行っているため、重要なのはできるだけ早く因数分解することではなく、異なる呼び出し間でできるだけ多くの作業を共有することです。試行分割を高速化するには、メモ化、および/または気になる最大数の平方根までの素数の事前計算を検討できます。すべての数値について sqrt(n) までループするよりも、素因数分解を取得してからすべての因数の合計を計算する方が高速です。
2 ^ 64 より大きいなど、非常に大きな友好的なペアを探している場合、少数のマシンでは、因数分解がどれほど高速であっても、すべての数値を因数分解することはできません。候補者を見つけるために使用しているショートカットは、それらを因数分解するのに役立つ場合があります。
もちろん、Hal Mahutan 教授 (2021 年 2 月) による HAL アルゴリズムもあり、これは因数分解研究の最先端にあります。
ここで最新のアップデートをご覧ください
https://docs.google.com/document/d/1BlDMHCzqpWNFblq7e1F-rItCf7erFMezT7bOalvhcXA/edit?usp=sharing
公開鍵の 2 つの大きな素数を解くと、次のようになります...
任意の AxB = 公開鍵は、正の X 軸と Y 軸上に描くことができ、すべての非整数要素が公開鍵を解決する連続曲線を形成します。もちろん、これは役に立ちません。現時点では単なる観察です。
Hal の洞察は次のとおりです。A が整数である点、特に A が整数のときに B が示す点のみに関心があると主張する場合。
このアプローチによる以前の試みは、数学者が B の残りの明らかなランダム性、または少なくとも予測可能性の欠如と格闘したときに失敗しました。
Hal が言っていることは、a/b の比率が同じであれば、どの公開鍵でも予測可能性は普遍的であるということです。基本的に、一連の異なる公開鍵が分析のために提示される場合、a/b が一定である処理中に同じポイントを共有する場合、つまり同じ接線を共有する場合、それらはすべて同じように処理できます。
私が描いたこのスケッチを見て、マフタン教授がここで何をしているのかを説明してみてください.
だから、ここにハルの天才があります。Hal は強力なスーパーコンピューターを利用して、一連のハッシュテーブル (図では Q、R、S & T) を生成します。左側の 3 つの A x B = Key 曲線でわかることは、それらがすべて接線 T と S を共有していることです (ここでハイライト表示されているものだけ)。接線が同じである場合、その領域を主宰するハッシュテーブルを共有できます。
技術的なメモとして、明らかに AxB= Key の曲線では、AxB の値を変更すると、物事は継続的に変化するため、理論的には、ハッシュテーブルにマップされる共有接線は時代遅れになりますが、興味深いことはは、非常に大きなキーを使用することです (皮肉なことに、ハッシュテーブルが役立つ場所でより長い実行を共有するため、これによりクラックが容易になります)。因数分解と計算の進歩が加速するにつれて、鍵のサイズがはるかに大きくなると予想されるため、これは素晴らしいニュースです。実際に起こることは、ハッシュテーブルが適用される接線が発散し始めると、ハッシュテーブルの予測可能性が文字通り「焦点から外れる」ことです。幸いなことに、これは問題ではありません。新しい Tangent に適切にマップされた次のハッシュテーブルにジャンプするからです。
明確にするために、これまでに生成されたすべての公開鍵は常に同じハッシュテーブルのセットを使用するため、オンラインで保存できる一種の 1 回限りの投資であり、文字通り数百万テラバイトのルックアップ データです。接線比。
では、素数の検出を加速するためにハッシュテーブルは何をするのでしょうか。ハッシュテーブルは、公開鍵を B で割った余りで索引付けされます。基本的に、Hal は、すべての鍵について、A x B の任意の比率を調べることができると言っています。同じタンジェントを共有する異なるカーブ間の唯一の違いは、「コントロール カーブ」によって決定される異なるオフセットが必要なことです。「制御曲線」は、適切な係数を生成する任意の曲線です。「コントロール カーブ」の場合、キーが 15 で、マッピングされるタンジェントが B = 5 の場合、A は 3 で残りはゼロです。同じ接線を持つ別のマッピングで、キーが 16 になったとします。B の場合は 5.33、A の場合は 3.2 にある同じ接線を見つける必要があります。したがって、A の残りは .2 です。
では、ハッシュテーブルには何がありますか? オフセットによってインデックスが付けられ、値は常に AxB 曲線に沿った距離を返します。この距離では、B の別の整数が見つかりません。Hal が言っているのは、先にジャンプして、それらの数値が要因であることを確認しない方が安全だということです。そして、それは基本的にそれです。ハルは、チェックする必要のないカーブに穴を開け、ゲーム全体をスピードアップします。
マフタン先生ありがとう！
まだフォローしている方のために、作業メモの一部を以下に示します。
高速因数分解攻撃アルゴリズムの箇条書き
Bを予測するアルゴリズムの戦略に従います
もう 1 つの観察結果は、「キー」の値が十分に大きい場合、整数の増分で A の値を変更する際に、A/B の比率または接線角度がほぼ同じままであることを観察することです。
この観察の重要な側面は、キーのサイズが大きくなるにつれて、各反復でタンジェントがより一定のままになることです。基本的に、このプロパティを使用するアルゴリズムは、キーのサイズが大きくなるにつれて効率が向上することを意味します。これは、キーのサイズを大きくすると要因を推測することが指数関数的に難しくなる従来のアプローチとは逆です。これは非常に重要なポイントです... (これについて詳しく説明してください。ニックにお願いします)
アルゴリズム自体は次のとおりです。
ここでの重要な概念は、ルックアップが不正確になる (そして焦点が外れる) 前に、A/B (接線) の比率が十分に近い「キー」に対してルックアップ テーブルを共有できるということです。
(また、キー サイズが変化すると、新しいテーブル セットが必要になるか、既存の比率テーブルを再利用するために適切なマッピングを作成する必要があることに注意してください。)
もう 1 つの非常に重要な点 Nick は、すべてのキーが同じルックアップ テーブルのセットを共有できるということです。
基本的に、ルックアップ テーブルは Key / A の計算からの剰余をマップします。剰余に関心があるのは、剰余 = ゼロの場合、A が既に整数であるため、剰余を実行したからです。
実際の余りを計算しなくても先にスキャンできるように、十分な数のハッシュ テーブルを用意することをお勧めします。1 兆から始めるとしましょう。しかし、どれだけの計算能力を持っているかによって、それは明らかに変化する可能性があります。
適切に近い A/b 比のハッシュテーブルは次のとおりです。テーブルは適切な剰余で索引付け (キー付け) され、任意の剰余での値は、剰余がゼロに触れることなく A * B 曲線に沿って移動できる「安全な」距離です。
0 と 1 の間で (疑似ランダムに) 振動する剰余を視覚化できます。
アルゴリズムは曲線に沿って数値 A を選択し、安全な距離を調べて次のハッシュテーブルにジャンプするか、少なくともアルゴリズムは次のハッシュテーブルが利用可能になるまでこれらの要素のチェックを行います。十分な数のハッシュテーブルがあれば、ほとんどのチェックを回避できると思います。
ルックアップ テーブルに関する注意事項。
どのキーについても、A/B 比に応じて適切にカーブするテーブルを生成できます。テーブルの再利用は不可欠です。推奨されるアプローチ N の平方根 (適切なキー) から A/B の一連の制御曲線を生成し、途中で 2 乗します。それぞれが前のキーよりも 0.0001% 大きいとします テーブルのサイズも A/B 比の 1% にします 余素係数を計算するときは、キーに一致する最も近いルックアップ テーブルを選択します。ハッシュテーブルへのエントリ ポイントを選択します。これは、テーブル内のエントリ ポイントと実際のエントリ ポイントとのオフセットを覚えておくことを意味します。ルックアップ テーブルは、対応するコプライム エントリがゼロに非常に近く、手動でチェックする必要があるエントリ ポイントの一連のポイントを提供する必要があります。シリーズのポイントごとに、適切な数式を使用して実際のオフセットを計算します。（これは積分計算になります。数学者に見てもらう必要があります）なぜですか？A/B が Key の平方根であるときに制御テーブルが計算されたためです。曲線に沿って移動するとき、適切に間隔をあける必要があります。おまけとして、数学者はキー スペースを A/B 曲線上の点に折りたたむことができます。その場合、1 つのテーブルを生成するだけで済みます。
重要な概念
数値 A x B = Key は次をプロットします。
X 軸マップ A と Y 軸マップ B。
A 軸および B 軸に対する曲線の近さは、公開鍵のサイズ (A x B = キー) に依存します。基本的にキーが大きくなるにつれてカーブが右にシフトしていきます。
さて、私があなたに消化したい、または質問したい考えは
これは機能しますが、遅すぎます。HAL アルゴリズムで行っていることは、上記の基本アルゴリズムを改善して、剰余がゼロに近づきすぎないことがわかっている曲線のチャンクをジャンプできるようにすることです。
ジャンプできる数が多いほど、アルゴリズムはより効率的になります。
改善された HAL アルゴリズムに入る前に、重要な概念を確認しましょう。
Key の値が非常に大きい場合 (A x B = Key を思い出してください)、A/B の比率はほとんど一定で、RSA キーは 2 乗 4096 であり、これは大きな値です。
特定の (大まかな) サイズのキー用に最適化された、クラウドに事前にロードされた一連のテーブルを作成したと仮定します。
この作品は重要なコンセプトです。