前を計算せずにPiのN番目の符号を計算

最近、Pi番号を計算するためのエレガントな公式があります。これは、1995年にDavid Bailey、Peter Borwine、およびSimon Pluffによって最初に公開されました。

画像



それのように見える:それについての特別なこと-円周率を計算するための式がたくさんあります。学校モンテカルロ法から、とらえどころのないポアソン積分および中世後期のフランソワ・ベト式まで。 しかし、特別に注意を払う必要があるのはまさにこの式です。これにより、前のものを見つけることなくpiのn番目の符号を計算できます。 これがどのように機能するかについての情報と、1,000,000番目の文字を計算する既製のC言語コードについては、habracatを要求します。



PiのN番目の符号を計算するアルゴリズムはどのように機能しますか?

たとえば、Piの1000番目の 16進数記号が必要な場合は、式全体に16 ^ 1000を乗算し、括弧の前の係数を16 ^(1000-k)に反転します。 べき乗の場合、 バイナリべき乗アルゴリズムを使用します。または、以下の例に示すように、 べき乗剰余を使用します。 その後、シリーズのいくつかのメンバーの合計を計算します。 さらに、ロットを計算する必要はありません。kが増加すると、16 ^(Nk)はすぐに減少するため、後続の項は目的の桁の値に影響しません。 それはすべて魔法です-独創的でシンプル。



Bailey-Borwine-Puff式は、2000年に世紀のトップ10アルゴリズムリストに含まれていたPSLQアルゴリズムを使用してSimon Puffによって発見されました。 PSLQアルゴリズム自体は、ベイリーによって開発されました。 ここに数学者に関するメキシコのシリーズがあります。

ちなみに、アルゴリズムの実行時間はO(N)、メモリ使用量はO(log N)です。Nは目的の文字のシリアル番号です。



私は、アルゴリズムの作者であるデイビッド・ベイリーによって直接書かれたCのコードを引用するのが適切だと思います。



/* This program implements the BBP algorithm to generate a few hexadecimal digits beginning immediately after a given position id, or in other words beginning at position id + 1. On most systems using IEEE 64-bit floating- point arithmetic, this code works correctly so long as d is less than approximately 1.18 x 10^7. If 80-bit arithmetic can be employed, this limit is significantly higher. Whatever arithmetic is used, results for a given position id can be checked by repeating with id-1 or id+1, and verifying that the hex digits perfectly overlap with an offset of one, except possibly for a few trailing digits. The resulting fractions are typically accurate to at least 11 decimal digits, and to at least 9 hex digits. */ /* David H. Bailey 2006-09-08 */ #include <stdio.h> #include <math.h> int main() { double pid, s1, s2, s3, s4; double series (int m, int n); void ihex (double x, int m, char c[]); int id = 1000000; #define NHX 16 char chx[NHX]; /* id is the digit position. Digits generated follow immediately after id. */ s1 = series (1, id); s2 = series (4, id); s3 = series (5, id); s4 = series (6, id); pid = 4. * s1 - 2. * s2 - s3 - s4; pid = pid - (int) pid + 1.; ihex (pid, NHX, chx); printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n", id, pid, chx); } void ihex (double x, int nhx, char chx[]) /* This returns, in chx, the first nhx hex digits of the fraction of x. */ { int i; double y; char hx[] = "0123456789ABCDEF"; y = fabs (x); for (i = 0; i < nhx; i++){ y = 16. * (y - floor (y)); chx[i] = hx[(int) y]; } } double series (int m, int id) /* This routine evaluates the series sum_k 16^(id-k)/(8*k+m) using the modular exponentiation technique. */ { int k; double ak, eps, p, s, t; double expm (double x, double y); #define eps 1e-17 s = 0.; /* Sum the series up to id. */ for (k = 0; k < id; k++){ ak = 8 * k + m; p = id - k; t = expm (p, ak); s = s + t / ak; s = s - (int) s; } /* Compute a few terms where k >= id. */ for (k = id; k <= id + 100; k++){ ak = 8 * k + m; t = pow (16., (double) (id - k)) / ak; if (t < eps) break; s = s + t; s = s - (int) s; } return s; } double expm (double p, double ak) /* expm = 16^p mod ak. This routine uses the left-to-right binary exponentiation scheme. */ { int i, j; double p1, pt, r; #define ntp 25 static double tp[ntp]; static int tp1 = 0; /* If this is the first call to expm, fill the power of two table tp. */ if (tp1 == 0) { tp1 = 1; tp[0] = 1.; for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1]; } if (ak == 1.) return 0.; /* Find the greatest power of two less than or equal to p. */ for (i = 0; i < ntp; i++) if (tp[i] > p) break; pt = tp[i-1]; p1 = p; r = 1.; /* Perform binary exponentiation algorithm modulo ak. */ for (j = 1; j <= i; j++){ if (p1 >= pt){ r = 16. * r; r = r - (int) (r / ak) * ak; p1 = p1 - pt; } pt = 0.5 * pt; if (pt >= 1.){ r = r * r; r = r - (int) (r / ak) * ak; } } return r; }
      
      





これはどのような機会を与えますか? たとえば、Piの数を計算する分散コンピューティングシステムを作成し、計算精度のためにHabrtの新しいレコード全体を設定できます(ちなみに、これは小数点以下10兆桁です)。 経験的データによると、Pi番号の小数部分は通常の数値シーケンスです(まだ確実に証明されていませんが)。つまり、それからの数字シーケンスは、パスワードと乱数だけの生成、または暗号化アルゴリズム(ハッシュなど)で使用できることを意味します。 それを使用する多くの方法があります-あなたはあなたの想像力をオンにする必要があります。



このトピックの詳細については、アルゴリズムとその実装( pdf )について詳しく説明しているDavid Bailey自身の記事を参照してください。



そして、あなたはRunetでこのアルゴリズムに関する最初のロシア語の記事を読んだばかりのようです-他の人を見つけることができませんでした。



All Articles