PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 64 / 奇数周期の平方根

平方根は連分数の形で表したときに周期的であり,以下の形で書ける。

√N = a0 + 1 / (a1 + 1 / (a2 + 1 / (a3 + ...)))

たとえば√23を考えよう。

√23 = 4 + √23 - 4 = 4 + 1 / (1 / (√23 - 4)) = 4 + 1 / (1 + (√23 - 3) / 7)

この操作を続けていくと次のようになる。

√23 = 4 + 1 / (1 + 1 / (3 + 1 / (1 + 1 / (8 + ...))))

操作をまとめると以下のようになる。

  • a0 = 4, 1/(√23-4) = (√23+4)/7 = 1 + (√23-3)/7
  • a1 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
  • a2 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
  • a3 = 1, 7/(√23-4) = 7(√23+4)/7 = 8 + (√23-4)
  • a4 = 8, 1/(√23-4) = (√23+4)/7 = 1 + (√23-3)/7
  • a5 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
  • a6 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
  • a7 = 1, 7/(√23-4) = 7(√23+4)/7 = 8 + (√23-4)

よって,この操作は繰り返しになることが分かる。√23 = [4;(1,3,1,8)] と表す。(1,3,1,8)のブロックは無限に繰り返される項を表している。

最初の10個の無理数である平方根を連分数で表すと次のようになる。

  • √2=[1;(2)], period=1
  • √3=[1;(1,2)], period=2
  • √5=[2;(4)], period=1
  • √6=[2;(2,4)], period=2
  • √7=[2;(1,1,1,4)], period=4
  • √8=[2;(1,4)], period=2
  • √10=[3;(6)], period=1
  • √11=[3;(3,6)], period=2
  • √12= [3;(2,6)], period=2
  • √13=[3;(1,1,1,1,6)], period=5

N ≤ 13で奇数の周期をもつ平方根は4つある。N ≤ 10000 について奇数の周期をもつ平方根が何個あるか答えよ。

Problem 64 - Project Euler


無理数を連分数展開してみたら [4;(1,3,1,8)] の形で答えを返してくれたので助かりました。周期が奇数かどうかの判定はこうやります。

  1. ContinuedFraction で連分数に直す
  2. Last で循環する部分をとりだす
  3. Length で長さを求める
  4. OddQ で奇数かどうか調べる


ルートの中身が平方数のときは Last で循環する部分をとりだすところで失敗するので,! IntegerQ[Sqrt[n]] で除外します。


あとは Select して Length を求めれば終わりです。