Horner法的な変形による数列総和の短縮

ゴルフでは常識に囚われてはいけないのですね!


と,言いつつもタイムリーな問題が出てたので短縮パターンの紹介をしてみる.Horner法は次のような変形で,
a_0 + a_1x^1 + Ba_2x^2 + \cdots + Ba_nx^n = a_0 + x(a_1 + x(a_2 + \dots))
Haskellだとfoldで以下のように書く.

-- f(x) = a0 + a1 * x + a2 * x^2 + ... + an * x^n
f x=foldr(\a b->a+x*b)0[a0,a1,a2,(略),an]

係数が逆順になってればfoldlでどうぞ.

これの一体何がタイムリーかというとanarchy golf 391. Ellipse circumferenceがコレに類似した変形を使う問題だった.タイトル通り楕円の周長を求めさせる問題で,出力は整数に丸めろとのこと.アプローチとしては次の3ケースが考えられるかな.

  1. 真面目に数値積分する
  2. 解析的に積分した結果であるΣ計算を適当に打ち切って使う
  3. 単純な近似式を使う

まず,単純な近似式については,いくら整数に丸めるといっても題意を満たす精度にならないのでボツ.数値積分も台形公式なら最低2回積分対象の式を展開しなければならないこともあり,(2つ目と比べ結果的に)短くならないか精度が足りないか速度が足りないのいずれかに悩まされることになる.なので,ここは2つめの無限列のΣ計算を使う.列が十分急速に小さくなっていってくれれば,途中の適当な位置で列を打ち切って足せばよい.

無限列の加算としては次のいずれかを使う.

  1. 2\pi a\Sigma_{n=0}^{\infty}\frac{1}{1-2n}(1-\frac{b^2}{a^2})^n\Pi_{m=1}^{n}(\frac{2m-1}{2m})^2
  2. \pi(a + b)\Sigma_{n=0}^{\infty}(\frac{\Gamma(\frac{3}{2})}{\Gamma(n + 1)\Gamma(\frac{3}{2}-n)})^2(\frac{a-b}{a + b})^{2n}

収束は2番目のほうがずっと速いが,言語によってどっちを選ぶと有利かは変わってくる.最終的な目安は冪乗演算で,冪乗が安い言語python(**演算子)やhaskell(^演算子)では後者,冪乗が高い言語C(power関数)やJavaScript(Math.pow関数)では前者有利になる.

さて,前者には\Pi計算(Haskell的にはproduct)が,後者には\Gamma関数(Haskell的には…ガンマ関数あったっけ?)があるので,こっちの方向性で本当に短かくなるのかと最初は疑問に思われるところだろう.今回のネタであるHoner法っぽい形へ式変形することによってこれら長そうな処理を削除しながら短縮を図ることができる.ではそれぞれ変形してみよう.

  1. 2\pi a + (1-\frac{2}{1} + \frac{0.75}{1^2})(1-\frac{b^2}{a^2})(2\pi a + (1-\frac{2}{2} + \frac{0.75}{2^2})(1-\frac{b^2}{a^2})(2\pi a + \dots))
  2. \pi(a + b)(\frac{(1-1.5)(a-b)}{1(a + b)})^2(\pi(a + b)+(\frac{(2-1.5)(a-b)}{2(a + b)})^2(\pi(a + b) + \dots))

次の項を今の項で割ると各段で掛け算する式が出てくる.この割る際にproductやガンマ関数等のめんどくさいやつらが相殺してくれて,項間の差に相当する多少の後腐れだけ残して消滅してくれる.これらの形はHorner法と少し違い,係数に相当する部分は一定で,項に応じて変動するのは係数にあたる部分でなく変数にあたる部分になっている.

これらの式をベースにコードに落とすと,やはりfoldrで,

-- 1番
round$foldr(\n x->2*pi*a+(n-2+3/4/n)/n*(1-b^2/a^2)*x)0[1..4e3]
-- 2番 (こちらのほうが収束が速いので丸めに対して打ち切れる点が近くに取れる)
round$pi*foldr(\n x->a+b+((a-b)*(n-3)/n/(a+b))^2*x)0[2,4..98]

となる.

よって,Haskellのコードとしては以下のものを投稿することになった.

m@main=getLine>>=print.round.f.map read.words>>m
f[a,b]=pi*foldr(\n x->a+b+((a-b)*(n-3)/n/(a+b))^2*x)0[2,4..98]

残念ながら他のHaskell Golferが参戦してないため比較対象が無いが.まぁ,短いのではないかと思う.


おまけとして,この問題に関しては他のHaskell Golferがかまってくれなかったので,さびしくなって他の言語にカラみにいってみた.結果的に無双になってしまったが,これらの言語でのゴルフは不慣れなので,もっと短くなる表現があるのかもしれない.

まずPythonについては前述した通り冪乗が安いので,2番の式をベースにしたコードになっている.

while 1:a,b=map(int,raw_input().split());x=n=13;exec"x=3.1416*(a+b)+((1-1.5/n)*(a-b)/(a+b))**2*x;n-=1;"*n;print"%.f"%x

foldr相当になるように,nが大きいところからループを回している.細かい特筆点としては,初期値xと円周率の2点がある.初期値xをループ変数と同時に13で初期化しているが,この初期値の影響は反復計算と丸めの中に消せてしまう.円周率はmathをimportするのは重いし,どうせ丸められるので十分な近似値3.1416を使った.別途,Pythonに関しては1番の式が選べない理由があり,execする文字列に4000も掛けると大変なことになるためだ.

次にC.Cでは冪乗のコストが高い.1変数の2乗ならまだしも,長い式に対する2乗ではpowerが必要になってくる.なので,1変数の2乗だけで済ませることができる1番の式をベースにしたコードを作ることになる.

float x,a,b;main(n){for(;~scanf("%f%f",&a,&b);printf("%.f\n",x))for(n=4e3;--n;)x=6.2832*a+(n-2+.75/n)/n*(1-b*b/a/a)*x;}

これも初期値xと円周率についてはPython版と事情が同じになっている.ただ初期値xについてはよりいっそう雑な扱いになっていて,外のループ毎に初期化すらしていない.これでも消えてくれて大丈夫なので.

最後にJavaScript.実は初めてのJavaScript Golfだった.やはりCと同様に冪乗コストが高いため,1番の式をベースにする.

for(;[a,b]=readline().split(x=' ');print(x*a+.5|0))for(b/=a,n=4e3;--n;)x=6.2832+(n-2+.75/n)/n*(1-b*b)*x

こまかい点はCと同じ.だが,JavaScript版ではタイムアウトとの戦いになる.途中のb/=aが特に顕著で,これを廃止して内部ループ内を(1-b*b/a/a)に戻したほうが短いのは明らかだが,文字列から数字への変換が4000回のループそれぞれで入るのかタイムアウトになってしまう.このコードでanarchy golfのタイムリミットの3秒ギリギリだ.というより,実際にこのコード投稿時に何度かタイムアウトしており,そのうち運良く通ったという指運ゲーっぷりだった(何度か試せば通るだろうという目算はあったが).


今回のパターンは,一般項が複雑であっても隣と割るとシンプルになってくれる(ことがある)というところがミソで,割とイロイロな計算にあてはまりやすい.特にHaskell使いはfoldでカッコよくコードに落としこむことができる.他の言語であっても,加算結果と次項の情報を1変数で抱えて反復でき,また,十分に長く足せば反復時の初期値の影響も消し飛んでしまうので短くなりやすいようだ.anarchy golfでも他にこのパターンで縮む問題があった気がするので,気になる人は見直してみるのもいいかもしれない.