Python

√2を高精度で計算する

昨日はcos20°を求めましたが、今日は√2を求めます。100万桁を目標にしましょう。 √nの循環連分数になります。√2は、 √2 - 1 = 1 / (2 + 1 / (2 + ...)) = [ 0; 2, 2, ... ] となります。連分数を途中で打ち切ったときの分数をan / bnと置くと、 an+2 = 2an+1…

cos20°を高精度で計算する

数学ガール/ガロア理論 (数学ガールシリーズ 5)作者: 結城浩出版社/メーカー: SBクリエイティブ発売日: 2012/05/30メディア: 単行本購入: 8人 クリック: 310回この商品を含むブログ (57件) を見る昨日から読んでいる。P.164に、 cos が扱える計算機を使えば…

関数のリスト

例えば、0に等しいか調べる関数、...、3に等しいか調べる関数をリストにしたいとします。 preds = [ lambda n: n == 0, lambda n: n == 1, lambda n: n == 2, lambda n: n == 3 ] print preds[0](3) # False print preds[1](3) # False print preds[2](3) # …

結合則が成り立たない場合のfold系の高速化

Modegramming Style: Scalazの型クラスこれを読んでいてふと思ったこと。 結合則が成り立つ場合、畳込み(fold)を行うときリストを分割して並列で計算させれば速くなります。例えばsumですね。1〜100を足し合わせるときに、1〜50と51〜100に分けてそれぞれ…

Project Euler 144 namedtuple

Problem 144ただの高校数学ですが、Pythonのnamedtupleを説明するために使います。namedtupleは中身はタプルでもクラスのように使えるものです。 from collections import namedtuple Point = namedtuple("Point", "x y") これだけでPointクラスのように使え…

Project Euler 141(3)

Problem 141sがtで割り切れるから、 p = s / t と置きます。そうすると元の式は、 m2 = pt3(u3 + p) となります。 ここでt = 1 p = 1とすると、 m2 = u3 + 1 これは因数分解ができるので簡単に解が求められますが、こういうケースは稀なので考えません。 次…

Project Euler 141(2)

Problem 141rを平方部分とそうでない部分に分けます。例えば、24なら22と6です。一般に r = s2t(tは平方成分を含まない) と書けます。そうすると、q, dは、 q = sut d = u2t 元の式は、 n = m2 = qd + r = st(tu3 + s) ここからtはsで割り切れることを示し…

Project Euler 141(1)

Problem 141 n = m2 = qd + r とします。rはdより小さいから、r2 = qdまたはq2 = rdとなります。前者なら、 m2 = r2 + r r2と(r + 1)2の間に平方数は無いので前者はありえません。 さて、 q2 = rd が成り立つので、例えば、r = 1なら、公比は2, 3, 4, …とな…

Project Euler 137

Problem 137 AF(x) = xF1 + x2F2 + x3F3 + ... xAF(x) = x2F1 + x3F2 + x4F3 + ... x2AF(x) = x3F1 + x4F2 + x5F3 + ... 引き算して、 (1 - x - x2)AF(x) = xF1 + x2(F2 - F1) + x3(F3 - F2 - F1) + ... Fk = Fk-1 + Fk-2 F1 = F2 = 1 を使って、 (1 - x - x…

Project Euler 136

Problem 136これは前々回の面倒な方法が使えます。nが2を含まないとき、4で割って余り3の素数pなら約数は1とpですが、 y = p 4d - y = 1 なら z = p - (p + 1) / 4 > 0 で成り立ちますが、 y = 1 4d - y = p なら z = 1 - (p + 1) / 4 ≤ 0 だから成り立ちま…

Project Euler 135(2)

Problem 135難しいことを考えずに大きさNのリストaを用意して、 y(4d - y) = n のyと4d - yを振ってその結果得られたnに対し、 a[n] += 1 とすればよいです。 0.8sでした。こんな実質10行に満たない単純な方法の方が速いんですね。

Project Euler 135(1)

Problem 135公差をdとすると、x = y + d, z = y - dだから、 x2 - y2 - z2 = -y2 + 4dy = y(4d - y) = n nを素因数分解して約数を生成します。その約数をd1としてd2 = n / d1とすると、 y = d1 4d - y = d2 4d = d1 + d2 だから、d1 + d2が4で割り切れてd1が…

Project Euler 131

Problem 131 n3 + n2p = m3 とおきます。ここで、nがpの倍数と仮定しましょう。n = kpとすると、 p3 + p = (m/k)3 p3と(p + 1)3の間に立方数はないので矛盾。よってnとpは互いに素です。 さて、nは立方数です。なぜなら、nに立方数でない因子があるとすると…

Project Euler 133

Problem 1335より大きい素数pに対し、A(p)が2と5のみから成り立っていなければなりません。そこで、p - 1から2と5を取り出して、例えばp = 73ならp - 1 = 23 * 32だから23 = 8ですが、 108 ≡ 1(mod 73) なので、A(p)は8の約数で、十分に大きな10nを割り切る…

Project Euler 132

Problem 1322, 3, 5以外の素数を列挙します。A(p)(すなわちB(p))が109を割り切る素数pを列挙すればよいです。109を割り切るかどうかの判定は比較的簡単です。まず、p - 1の2と5の成分を取り出します。例えばp = 1920001として、p - 1 = 210 * 3 * 54だから…

Project Euler 130

Problem 130前回の定義でB(n)は、 10e ≡ 1(mod n) を満たす最小のe(> 0)としました。これをnの既約剰余類群での10の位数と言います。これを効率的に求めればよいです。nを素因数分解して、 n = p1e1...pmem とすると、 B(n) = lcm(B(p1e1), ... , B(pmem)) …

Project Euler 129

Problem 129A(n)は、 (10e - 1) / 9 ≡ 0 (mod n) となる1以上で最小のeです。 nが2または5で割り切れるとき、この式を満たすeは無いので考慮に入れません。 A(n)は、nが3で割り切れないとき、 10e ≡ 1 (mod n) となる最小のeです。7の剰余なら、10から10倍し…

Project Euler 127

Problem 127 9秒でした。

Project Euler 124(2)

Problem 124もっと速い方法を考えましょう。10万で9ms、1000万で1.1sでした。

Project Euler 124(1)

Problem 124何の工夫もしていませんが、0.4sです。

Project Euler 122(2)

Problem 122200で1.9s、400で20sでした。

Project Euler 121(2)

Problem 1212000で5秒くらいでした。

Project Euler 121(1)

Problem 121この問題は解くだけなら簡単です。

ラマヌジャンの高度合成数

英語ではHighly composite numberです。これより小さい数は全て、約数の個数が自分のそれより小さい自然数です。なぜこれが重要なのかは知りません。これはProblem 110のコードを少し改変するだけで得られます。1分以内に約数の個数が1012以内の高度合成数な…

Project Euler 107(3) ブルーフカ法

Problem 107今度はブルーフカ法です。日本語での解説ページが見当たらなかったので少し丁寧に説明します。まず、クラスカル法と同様にノード一つの木を全てのノードについて作ります。次に全ての木について、他の木と繋がる辺の中で最も重みの小さい辺を選び…

Project Euler 107(2) クラスカル法

Problem 107次はクラスカル法です。まず、一つのノードのみの木を全てのノードについて作ります。すなわち例では木が7つ作成されることになります。ノードはルート(根)がすぐにわかるようにしておきます。ノードに親を持たせれば、親を辿っていくとルート…

Project Euler 107(1) プリム法

Problem 107この問題は「最小全域木」を求める問題です。最小というのは重みの和が最小ということで、全域というのは元のグラフの全てのノードを持つグラフということで、木というのはループが無い連結グラフということです。この用語さえ知っていればあとは…

PriorityQueueで比較法を指定する

Problem 70はフォーラム向けにまずPythonで解きました。そのとき、最初はPriorityQueueにつっこむ値をdoubleにしていました。しかし、これは正確ではありません。doubleで表現すると同じ値でも実は違う値ということもあり得ます。もっともそれが両方とも並べ…

ビットごとのエラトステネスのふるい

Pythonで何も考えずにエラトステネスのふるいのコードを書くとこんな感じでしょうか。 def sieve(max_n): a = [ True ] * L for p in takewhile(lambda n: n * n < L, (n for n in xrange(2, L) if a[n])): for k in xrange(p * 2, L, p): a[k] = False Bool…

setで違う値を同じとみなす

タイトルの意味がわかりにくいですが、 s = set([((2, 3), 1), ((2, 4), 2), ((2, 3), 3)]) print s set([((2, 4), 2), ((2, 3), 1), ((2, 3), 3)])となりますが、タプルの第2項は無視して第1項だけで同じ値か否かを判定したいときにこれでは困るという場合…