【Project Euler】Problem 26: 循環小数

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題 26：循環小数

原文 Problem 26: Reciprocal cycles

問題の要約：分母が1000以下の単位分数（分子が１のもの）で循環小数となる循環節が最も長いものの分母を求めよ

循環小数(Wikipedia)に詳細がありますが、分数は有限小数（Repeating decimal）か循環小数（Rerminating decimal）になり、ここでは$1/d$の形の単位分数で循環小数になるものの循環節の長さを求める問題です。

まず単位分数を小数に変換するプログラムdecUnitFracです。Wikipediaの筆算にあるように「余りに同じ数が現れた時点で、繰り返しに入ったことがわかる」ので余りをすべてリストに入れてチェックします。小数部分の文字列(dec)と循環節の始まりの位置(rs)を返します。有限小数の時はrs=-1とします。問題の例の続きのd=11~23を表示して正しそうなことを確認します。

# Input: 
#   d - denominator of a unit fraction
# Output: 
#   dec - decimal string
#   rs - starting positon of reptend, rs<0 if terminating decimal:
def decUnitFrac(d):  
  m, mlist, dec = 1, [], ""
  while True:
    if m in mlist: # check the same remainder in the list)
      return (dec, mlist.index(m))
    mlist.append(m)               # remember all reminders
    dec += str((m*10) // d)       # repeating part string
    m = (m*10) % d                # m: reminder
    if m==0: return (dec, -1)  # terminating decimal

for d in range(11, 23+1):
  (dec, rs) = decUnitFrac(d)
  if rs>=0:
    print(f"1/{d}: 0.{dec[:rs]}({dec[rs:]}) ")
  else:
    print(f"1/{d}: 0.{dec}")
#1/11: 0.(09) 
#1/12: 0.08(3) 
#1/13: 0.(076923) 
#1/14: 0.0(714285) 
#1/15: 0.0(6) 
#1/16: 0.0625
#1/17: 0.(0588235294117647) 
#1/18: 0.0(5) 
#1/19: 0.(052631578947368421) 
#1/20: 0.05
#1/21: 0.(047619) 
#1/22: 0.0(45) 
#1/23: 0.(0434782608695652173913)

分母が素数の単位分数の循環小数

問題を解くだけであればこれを1000まで回して循環節の長さの最大値を求めれば良いのですが。いろいろ気付いたのでちょっと考察してみます。

まず循環小数の結果を眺めてみると、分母が素数の時に循環節の長くなっているのが分かります。そこで素数の場合だけ表示してみます。$rc$が循環節の長さです。

1/13: rc=6, 0.(076923) 
1/17: rc=16, 0.(0588235294117647) 
1/19: rc=18, 0.(052631578947368421) 
1/23: rc=22, 0.(0434782608695652173913) 
1/29: rc=28, 0.(0344827586206896551724137931) 
1/31: rc=15, 0.(032258064516129) 
1/37: rc=3, 0.(027) 
1/41: rc=5, 0.(02439) 
1/43: rc=21, 0.(023255813953488372093) 
1/47: rc=46, 0.(0212765957446808510638297872340425531914893617) 

すると素数でも長いもの(17,19,23,29,47)と短いもの(13,31,37,41,43)があり、長いものは必ず$rc=d-1$となっていることが分かります。

###循環節の長さと位数・原始根

これはなにか規則性がと思って調べてみたら、英語のWikipedia Repeating_decimal-Fractions with prime denominatorsに「1/pの循環節の長さは10 mod pの位数(order)に等しい」との記述がありました。位数はかなり難しい概念ですが「原始根の定義と具体例（高校生向け）」の説明が比較的分かりやすいと思います

$r,r^2 ,r^3 ,⋯ $と $r$ のベキを増やしていき，はじめて「$p$ で割った余りが $1$ となる」ような指数のことを $r$ の位数と言います。

これをもとに以下のようにプログラムで確認すると、位数（初めて１になる指数）はorder(10)で表示されていますが、循環節の長さと一致していることが分かります。

for p in [13,17,19,23,29,31,37,41,41,43]:
  pm = [pow(10,i,p) for i in range(1,p)]
  print(f"p={p}, order(10)={pm.index(1)+1}, {pm}")

p=13, order(10)=6, [10, 9, 12, 3, 4, 1, 10, 9, 12, 3, 4, 1]
p=17, order(10)=16, [10, 15, 14, 4, 6, 9, 5, 16, 7, 2, 3, 13, 11, 8, 12, 1]
p=19, order(10)=18, [10, 5, 12, 6, 3, 11, 15, 17, 18, 9, 14, 7, 13, 16, 8, 4, 2, 1]
p=23, order(10)=22, [10, 8, 11, 18, 19, 6, 14, 2, 20, 16, 22, 13, 15, 12, 5, 4, 17, 9, 21, 3, 7, 1]
p=29, order(10)=28, [10, 13, 14, 24, 8, 22, 17, 25, 18, 6, 2, 20, 26, 28, 19, 16, 15, 5, 21, 7, 12, 4, 11, 23, 27, 9, 3, 1]
p=31, order(10)=15, [10, 7, 8, 18, 25, 2, 20, 14, 16, 5, 19, 4, 9, 28, 1, 10, 7, 8, 18, 25, 2, 20, 14, 16, 5, 19, 4, 9, 28, 1]
p=37, order(10)=3, [10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1, 10, 26, 1]
p=41, order(10)=5, [10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1]
p=41, order(10)=5, [10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1, 10, 18, 16, 37, 1]
p=43, order(10)=21, [10, 14, 11, 24, 25, 35, 6, 17, 41, 23, 15, 21, 38, 36, 16, 31, 9, 4, 40, 13, 1, 10, 14, 11, 24, 25, 35, 6, 17, 41, 23, 15, 21, 38, 36, 16, 31, 9, 4, 40, 13, 1]

ここで元の問題に戻って、なるべく循環節の長いものすなわち位数が最大（p-1）のものを探せばよいということになりますが、前出のリンクの以下の定義によると結論として10が素数pの原始根となるものを探せばよいことが分かります。

原始根は 「位数が $p-1$ であるもの」と簡潔に表現できます。

10が原始根になる素数p

原始根の求め方は位数を求めて$p-1$となるものでも良いのですが、sympyにis_primitive_rootというその物があるのでそれを使います。したがって問題の答えは、

• 1000から大きい順に素数をリストして初めて10が原始根になる素数が、循環節の長さが最大のものということになります。

素数の逆順の生成にやはりsympyのprevprimeを使うと以下のようになります。これであれば1000が$10^6$になっても1秒以内に答えが出せます。

from sympy.ntheory import is_primitive_root
from sympy import prevprime
p = 1000
while True:
  p = prevprime(p)
  if is_primitive_root(10, p):
    break
print(f"Answer : {p}")

1/nの循環小数の長さを求める

ここまでで基本的なことは分かったので$\Large \frac{1}{n}$を小数にした時の循環節の長さ$L(n)$を求めるコードを以下の手順で作っていきます。

1. $n$が素数のとき
2. $n$が素数のべき乗のとき
3. $n$が合成数のとき

1. nが素数のとき

これは既に出た「1/pの循環節の長さは10 mod pの位数(order)に等しい」」なので、sympyのn_order関数を使って求めます。

import sympy as sp
for p in [3,11,13,17]:
  L = sp.n_order(10, p)
  print(f"1/{p}の循環節の長さ={L}")
# 1/3の循環節の長さ=1
# 1/11の循環節の長さ=2
# 1/13の循環節の長さ=6
# 1/17の循環節の長さ=16

2. nが素数のべき乗のとき

この時も位数(order)を調べれば良いので求めてみると。

$p$	$L(1/p)$	$L(1/p^2)$	$L(1/p^3)$	$L(1/p^4)$	$L(1/p^5)$
3	1	1	3	9	27
7	6	42	294	2058	14406
11	2	22	242	2662	29282
13	6	78	1014	13182	171366
17	16	272	4624	78608	1336336

なんとなく規則性がありそうですが分かりづらいので積の形で書くと規則性が見えてきます。

$p$	$L(1/p)$	$L(1/p^2)$	$L(1/p^3)$	$L(1/p^4)$	$L(1/p^5)$
3	1	1	1x3	1x3x3	1x3x3x3
7	6	6x7	6x7x7	6x7x7x7	6x7x7x7x7
11	2	2x11	2x11x11	2x11x11x11	2x11x11x11x11
13	6	6x13	6x13x13	6x13x13x13	6x13x13x13x13
17	16	16x17	16x17x17	16x17x17x17	16x17x17x17x17

すなわち、

L(1/p^k) = L(1/p) \times p^{k-1}

と表せそうです。ただし$p=3$の時は$p$が一つ少ないのが気になります。
これは3だけが例外なのか調べると以下に記述がありました。$3, 487, 56598313$の３つだけが例外ということです。

Reciprocals of composite integers coprime to 10
There are three known primes (3, 487, and 56598313) for which this is not true, and for those the period of $1/p^2$ is the same as the period of $1/p$ because $p^2$ divides $10^{p−1}−1$.

3. nが合成数のとき

最後に合成数の場合です、これも同じリンクに記述があります。結論から言うと以下のようにそれぞれの素因数の循環節の長さを求めて、その最小公倍数（LCM）を求めれば良いということです。

L(1/pq) = LCM(L(1/p),L(1/q))

1/nの循環小数の長さを求めるプログラム

以上の3つをもとにしたpythonの関数 L_recip です。意外に簡潔に書けました。$3 \le n \le 99$の結果を表にしたものを添付しました。

import sympy as sp
import math
def L_recip(n):
  fct = sp.factorint(n)
  ret = 0
  for p in fct.keys():
    if p in [2,5]:continue
    d = [1,2][p in [3, 487]]     # 2: 3, 487, 1: others
    rc =  sp.n_order(10, p)* p**max((fct[p]-d),0)
    ret = math.lcm(max(ret,1), rc)
  return ret

	0	1	2	3	4	5	6	7	8	9
0				1	0	0	1	6	0	1
10	0	2	1	6	6	1	0	16	1	18
20	0	6	2	22	1	0	6	3	6	28
30	1	15	0	2	16	6	1	3	18	6
40	0	5	6	21	2	1	22	46	1	42
50	0	16	6	13	3	2	6	18	28	58
60	1	60	15	6	0	6	2	33	16	22
70	6	35	1	8	3	1	18	6	6	13
80	0	9	5	41	6	16	21	28	2	44
90	1	6	22	15	46	18	1	96	42	2

（開発環境：Google Colab）

この考え方は Project Euler Problem 417: Reciprocal Cycles IIを解くのに役に立ちます。