*/, /* x < 0.5のとき収束速度が落ちるので、   factorial(data, 1000); Algorithm for Calculating Individual Hexadecimal Digits of Pi. であることを使う.$x<-2$のときも$e^{-2}$を考えて符号を変えて同様にできる.    (2t+1)2 - (2s)2 2/(2n+1) = 1   πの値を求めるには、いろいろな式があるが、arctan(逆正接関数)を無限級数で表した式がよく使われる。式が単純な形で扱いやすく、ある程度の桁数の値を得るには適している。しかし欠点は、非常に長い桁数を得るには時間がかかりすぎることである。 $f^{(1)}(x)=f'(x),f^{(2)}(x)=f''(x)$である. 級数展開出来ればあとは任意項まで計算することで関数値を近似することができる(かもしれない).   変数名の整理: p = ak  (a、k、p はすべて正の整数) とかく.数列{$a_n$}に対して,その形式的な和$a_0+a_1+a_2+\cdots+a_n+\cdots$を級数といい. 有限個の数,$a_0,a_1,a_2,\cdots,a_n$の和を. である.つまり$n$項までの和を計算してから$n\to\infty$とすればよい.例えば, 関数$f(x)$は区間$J$で定義されているとし、$a\in J$ とする.いま, $a\in I$ $\subset J$を満たすある開区間$I$において,$f(x)$が$x-a$のベキ級数で表されているとする.つまり. </問題 #10241>, この問題の難しいところはやはり大きさの一点につきるだろう。10進36桁の整数はいまのプログラミング環境ではかなり厳しい。, <番外編> 最終的に、 具体的にはmath.hで宣言されている関数とほぼ遜色ないレベルで近似できることが目標である. これは級数の収束と関わりがある. What is going on with this article? * e^x = e^2 * e^(x-2) を考える. ソースコード 2.1.      π/4 = 12 arctan (1/49) + 32 arctan (1/57) - 5 arctan (1/239) + 12 arctan (1/110443)   1. p に対し、 p ← 1 と初期値を代入。 テイラー級数はベキ級数であり,$e^x$と$\log_ex$の収束半径はそれぞれ$\infty$,$|x-1|<1$である.$e^x$はすべての実数でテイラー級数が収束するが,$\log_ex$は$0 < x < 2$でしか収束しないのでより考察する必要がある. $f(x)$が$x=a$でテイラー展開できるとは, ($(e^x)'=e^x$が分かっているなら,中段の本題まで読み飛ばしてよい)      d = 2 ⇒ a = 3, b = 5  上の例では、級数の末尾項が999100なら、個々の項の指数部の値をe+300に統一して、仮数部の値を合計する。 下方からの理論構築を頭の中ですることはやはり重要だと思う。.   ここで使う式は、最もよく知られている、J.    X = s2        (2.2) p = p*3 = 81*3 = 243    (= 35) 当たりから、その後の合計計算を打ち切る。基数やkの値が大きければ大きいほどすぐに打ち切ることになる。, また、合計の計算では、級数の末尾から項ひとつひとつ足していくが、個々の項の指数部は、末尾項の指数部の値に(つまり、最大指数部)に合わせて、仮数部の小数点をシフトさせる。        (2.1) p = p*p = 3*3 = 9      (=32) が得られる。, 整数のべき乗、320等、多少高速にする計算。  整数dが与えられたら、素数 a とそのつぎの素数 b との間隔 (=b-a) がd以上であるような、最小素数 a, bペアを計算せよ。ただし、a, bの大きさは10進35桁まで。 <用例>   2 の平方根を小数点以下1000桁で求める例   int data[1000];   squareRoot(data, 1000, 2); 2 の平方根の値 (小数点以下1万桁)   3 の平方根の値 (小数点以下1万桁)   5 の平方根の値 (小数点以下1万桁)   6 の平方根の値 (小数点以下1万桁)   7 の平方根の値 (小数点以下1万桁)   8 の平方根の値 (小数点以下1万桁)   10 の平方根の値 (小数点以下1万桁), なお、初期値 p0、q0 は ペル方程式 p02 - M q02= 4 を満たす整数解。ただし、M は計算しようとする平方根の2乗の値。, となる。したがって平方根の計算において、本アルゴリズムの方法ではニュートン法等ほかの方法よりも効率がよく、数回の計算で必要とする精度のものが得られる。, をまともにと計算していくと、かなりのコンピュータでも、すぐにオーバフローしてしまう。ということで、本関数のようにそれなりの工夫が必要。, <関数名> 指数関数$e^x$は全ての$x\in\mathbb{R}$で微分可能である. Machinの式 * logx = log1/(1/x) = log1-log1/xを使って0 < x < 2の範囲に落とし込む. ベキ級数の収束に関して次の3つの場合しか起こらない. </番外編>, ここで考えている問題では、n=0の場合が Square Triangular Number にあたる。, 問題解くための数式を整理してみる。T = t (t+1)/2、S = s2 とおくと、   int data[1000]; (1)のとき収束半径は0, (3)のとき収束半径は$\infty$とする.        (2.1) p = p*p = 243*243 = 59049   (= 310) What is going on with this article? 1. が成り立つ.したがって$f'(x)=e^x$を得た.    (t2+t) (2n+1) / 2 = s2 の値を多桁数で算出する, <引数>   三角数 (Triangular Number) とはつぎの数列で表される数のこと、つまり、 1, 3, 6, 10, 15, 21, 28 など、k( k+1)/2 を満たす数のことをいう。一方、平方数(四角数ともいう) (Square number) とは 1, 4, 9, 16, 25 のような数のことをいう。   クリンジェンシェルナの公式   2. 20の2進数は 10100 であるので、 1:いま,$x>2$とする.対数関数の基本性質$(b)$を使うことによって次のようにできる. 実際に実行してみると. 以下の順序で話していく. とかく.$\sum_{k=0}^\infty a_k$が意味をもつのは,収束するときだけであることを忘れてはならない.    (2t+1)2-1 = 8 s2/(2n+1)      これで終了。, <計算量> 世の中にはこんな不思議な式があります. </オリジナル問題>, ということで、プログラミング問題を解いたといっても、数学と違って、いろいろなレベルがある。つまり、計算量を探求するならば、多くの問題はいまだに解けられていない、と強引に主張しても間違いない。しかし、コンピュータの世界では、そのための言い訳というか、責任逃れの逃げ道が理論的にある程度用意されている。NP完全問題がその好例。誰もが解けないから、私ができなくても当然という、おかしな論理が、コンピュータの世界では堂々と成り立つ。, <問題 #10241> 数学科卒。数論好き。   高野喜久雄の公式 C/C++の数値の0,NULL,空文字('\0'),空文字列("")の違いがよくわからなくなったので整理する。 内部的な値 まず,これらの内部的な値を以下のプログラムで確認する。 null.c/// \file null.c#include #define PRINT(x) printf(#x":%x\n */, デバイスでのパフォーマンス分析を自動化する新しいツールArm Mobile Studio, 想定読者:高校数学がある程度理解できていれば何とかなると思う.(論理ギャップは埋めた).   n   (入力) 階乗 n の値. 入力した自然数の階乗を求める.    X = (2n+1) t (t+1)/2 $(e):$ $(a)$において$p=0$とすると$\log_ea^0=\log_e1=0$を得る. 実際,$f(x)=e^x$とおくと,$f(x)$は全ての実数で定義されていて,微分の定義より, $p(h)=e^h-1,q(h)=h$ とおくと $h=0$で$p(0),q(0)$は明らかに連続で, $p(0)=q(0)=0$ , $q'(0)=1$より ロピタルの法則 が使えて. の2条件が成り立つことである.    ((2t+1)2-1) (2n+1) / 8 = s2   998100 仮数部 0.8185668046 指数部 e+300    Y(1) = 1, Y(2) = 36      (2.2) そのビットが1なら、さらに a を掛ける p ← p*a, <例> $x=a$におけるテイラー級数の収束半径が0でない. 級数の収束速度.極端な値をとると多くの項を計算しなければ正確な近似値が得られない.    Y(N) = 34 * Y(N-1) - Y(N-2) + 2 Why not register and get more from Qiita? 2.   keta (入力) 小数点以下の桁数の指定, の式を利用して計算。さらに、右辺の第2項目以降の計算では、計算回数の少ないつぎの式に変形して使う。, (((...((1/n + 1) 1/(n-1) + 1) 1/(n-2) + 1)...) 1/3 + 1) 1/2 + 1, <引数>   factorial ---- 階乗 n! 単純計算では、上の例だと、3を20回乗算すればいいわけだが、以下のようにすると5回の乗算で済む。, セコイ考えだが、64ビット整数でもべき乗の指数は大きく取れない(2のべき乗でも、最大64以下だから)ので、実質はあまり意味のないことかもしれない。考え方やC言語プログラムを以下に示す(メモ程度)。, <考え>   arctanをによる古典的な公式にほかに以下のものがある。 かなり野暮ったい級数が出てきたが,まずやることはひたすら微分することである. (2)の$r$のことを,ベキ級数$\sum_{n=0}^\infty c_nx^n$の収束半径という.      π/4 = 4 arctan (1/5) - arctan (1/239) $(3):\sum_{n=0}^\infty c_nx^nはすべてのxで収束する.$ <形式>   int squareRoot(int *data, int keta, int n); <引数>   data (出力) 平方根の値が10進4桁ずつ入った整数の配列   keta (入力) 小数点以下の桁数の指定   n   (入力) 計算しようとする平方根の2乗の値 (1 〜 99), <関数値>   正常に計算できたときは、計算結果となる data の長さ。エラーのときは 0。.        (2.1) p = p*p = 1*1 = 1    t (t+1) (2n+1) / 2 = s2  合計値は級数の後ろの項によってほぼ決まるので、級数の末尾から計算することにしよう。, 例えば、999100、998100の値を $(1):\sum_{n=0}^\infty c_nx^nは,x=0でしか収束しない.$ $(2):あるr,(0 < r <\infty)が存在して,\sum_{n=0}^\infty c_nx^nは,|x| < rならば収束し,|x| > rならば発散する.$ $(3):\sum_{n=0}^\infty c_nx^nはすべてのxで収束する.$      d = 6 ⇒ a = 23, b = 29 $(c)$から,$\log_eab^{-1}=\log_ea+\log_eb^{-1}.$また,$(a)$より$\log_ea+\log_eb^{-1}=\log_ea-\log_eb$を得る.      左から2ビット目のビット0に対し、 上式の通り、ごく簡単な、ベキ級数の合計を算出する問題。変数 l, h, k の値はともに1~15000000、hとlとの差は1000以内。, kの値が10以内なら、探せば、それなりの計算式が出てくるが、1500万ではそんな式は存在しないだろう。, 力任せで、C言語の数学関数をそのまま使うと、無限大になってしまい、計算不能になる。, 考え方 テイラー級数が意味をもてば,よく分からない関数値が実数係数多項式という非常に簡単なもので近似出来るのが,テイラー展開の醍醐味である.   999100 仮数部 0.9047921471 指数部 e+300 という形の級数を$x$のベキ級数という.ベキ級数がある$x_0\neq0$で発散すれば,$|x|>|x_0|$を満たす全ての$x$で発散する. 入力値が自然数前提のため,負数入力は考慮していない.必要に応じて,条件分岐で対処すること., 配列を用いることで,long int型に格納できないような,より大きな数を表現する..     X = (2n+1) * T  (ただし、n = 0, 1, ..., 7)  例: X = 1, 9, 36, 196, 225, 324, 900, 1225, 4225, 11025, 41616, 53361, 88209, 108900, ... なので、   970100 仮数部 0.0475525075 指数部 e+300, ベキ乗 xk → 仮数部 pow(10, modf(k*log10(x), &e)) 指数部 左の e, 最後に、元問題の検証データを載せておく(入力データの各行は左から、l, h, kの値), 16進数100万桁表示の円周率計算プログラムについて考える。実行速度は2秒以内としたい。効率のよい計算式と、FFTによる乗算で高速化にチャンレンジしたい。, ここのブログでも紹介している多倍長整数を利用するが、16進数に合わせて、1語の大きさを2進32ビット、つまり、16進8桁に直す。10進数と違って、除算ではビットシフト、剰余計算ではビットマスクだけでできるので、さらなる高速化が期待できそう。, 整数一つの配列はこんなものになるだろう。オーバーフローが起きて、64ビット整数にしないといけない可能性もあるが。, ネット上を調べたら、全ての桁をいっぺんに計算するのではなく、特定の桁だけを算出する方法もあるらしい。今回の用途に合致するかもしれない。, 英文 Algorithm for Calculating Individual Hexadecimal Digits of Pi, 複素数や無理数にまで拡張したら面白くないので、今回は、2次多項式の係数すべては整数ということに制限する。定理により、整数係数の多項式が有理数体で因数分解できるなら、因数の係数もすべて整数にすることができることから、因数の係数すべても整数とする。, 前回の記事 「大きさとの戦い」 では結局解き方が判らなくて、そのまま、中途半端に終わってしまった。, その時はまだ簡単に解けそうな問題が結構残っていたのだが、ここに来ると一問一問解くのに益々時間がかかる気がしてきて、再チャレンジしてみた訳だ。, (2t+1)2 - (2s)2 2/(2n+1) = 1  (n = 0, 1, 2, 3, 4, 5, 6, 7), というものだった。X = 2t+1とし、n に 0~7の値をそれぞれ代入すると、つぎの8つの方程式が得られる。    X2 - 8s2 = 1    X2 - 8/3 s2 = 1    X2 - 8/5 s2 = 1    X2 - 8/7 s2 = 1    X2 - 8/9 s2 = 1    X2 - 8/11 s2 = 1    X2 - 8/13 s2 = 1    X2 - 8/15 s2 = 1, 分母を払うために、s にそれぞれ、以下のものを代入する。すると見事にペル方程式にたどり着く。, 右列の Ans 部分は、ペル方程式の一般解 (X, Y) のうちの Y の値を、問題の解答に対応させるためのもの。, 結論的にいうと、方程式の数は9個と多いが、それぞれに簡単に最小解が得られる。さらに、多倍長整数とペル方程式の解の漸化式を使えば、一般解もすぐに判る。, たとえば、2番目のペル方程式 X2 - 24Y2 = 1 の最小解は (5, 1)。一般解は, Xn+1 = 5 Xn + 24 Yn  Yn+1 = 5 Yn + Xn  X1 = 5  Y1 = 1, になり、最小解(初期解)から漸化的に計算していくと、  (X, Y) = (5, 1), (49, 10), (485, 99), (4801, 980), ...  Ans = 9Y2 = 9, 900, 88209, 8643600, ... が得られる。, 解答をすべてここに書くのはまずいが、条件を満たす数は計 112 個。ラストの10個は以下となっている。, 73702684612392858376746072080400 151046383493325234090009219613956 663812918895887474609694358375209 1359417451439927106810082976525604 5131130648390546663702275158894481 8758618149740884101499623707907225 35524541072381125235676320801486025 46180175835514919973320476430050329 60153048103383854522005973075150400 65046891745147247744323041849672900, 出題側として、難しい問題にするには、自分や仲間達の最新研究成果を使うが、問題を難しくするには、計算すべき数値の範囲(大きさ)や数値の量を大幅に増やす、といった安易なアプローチもある。, 例えば、ごく簡単な足し算。数学では、a と b が与えられたら、その和は a+b で計算できるから、それ以上のことを深く考える必要はない。しかし、コンピュータの世界では、10進20桁未満の整数同士なら簡単だが、それ以上の大きさ、例えば、10進200桁同士の整数 a と b を足すにはひと工夫が必要。, 近々、2進128桁までの整数はC言語でもふつうに扱えるようになると予想できるが、ある種の計算ではいままで以上に難しくなる。例えば、素数表の生成には、エラトステネスのふるいが効率がもっともよいのでよく利用されるが、128桁になると、あまり役に立たなくなるかもしれない。10秒では計算が終了しないから。, <オリジナル問題>(素数の分布)      左からの5ビット目のビット0に対し、   data (出力) 円周率の値が10進4桁ずつ入った整数の配列 対数関数$\log_ex$は$x>0$で微分可能である. $f''(x)=-\frac{1}{x^2},f^{(3)}(x)=\frac{2}{x^3},f^{(4)}(x)=-\frac{6}{x^4},f^{(5)}(x)=\frac{24}{x^5},\cdots,f^{(n)}(x)=(-1)^{n-1}\frac{(n-1)!      T = T - X * (Y-A)2、X = 2 * X, すると、πの値は(A + B)2/(4 * T) となる。ただし、A、B、T、Yの各値は、求めようとする精度以上で計算しておく必要がある。, 1 から 99 の間の平方根の値を多桁数で算出する関数。多くの桁数の平方根が必要な時に利用するといいだろう。, 平方根の算出は一般的に難しいけど、結果があっているかどうかの検証はものすごく簡単。2乗してみればすぐに判るからだ。円周率 π とは大違い。, <関数名>   squareRoot ---- 1 から 99 の間の平方根の値を多桁数で算出する。.