Puzzle DE Programming / 循環小数 (repeating decimal)

パズルの説明

分数を小数に直すことを考えます。1 / 2, 1 / 5, 1 / 8, 1 / 10 などは　0.5, 0.2, 0.125, 0.1 のように途中で割り切れて、有限な桁で表すことができます。これを「有限小数」といいます。ところが、1 / 3, 1 / 6, 1 / 7 は、次のように有限な桁では表すことができません。

\(\begin{array}{l} \dfrac{1}{3} = 0.333333\ldots \\ \dfrac{1}{6} = 0.166666\ldots \\ \dfrac{1}{7} = 0.142857142857142857\ldots \end{array}\)

1 / 3 は 3 を無限に繰り返し、1 / 6 は 0.1 のあと 6 を無限に繰り返し、1 / 7 は数列 142857 を無限に繰り返します。このような少数を「循環小数 (repeating decimal)」といい、繰り返される数字の列を「循環節」といいます。分数を小数に直すと、有限小数か循環小数のどちらかになります。

循環小数は次のように循環節の始めと終わりを点で示します。

\(\begin{array}{l} \dfrac{1}{3} = 0.\dot{3} \\ \dfrac{1}{6} = 0.1\dot{6} \\ \dfrac{1}{7} = 0.\dot{1}4285\dot{7} \end{array}\)

このほかに、循環節に下線を引いたり、カッコで囲む表記法もあります。

それでは問題です。

分数を循環小数に変換するプログラムを作ってください。
循環小数を分数に変換するプログラムを作ってください。
分子が 1 の分数を「単位分数」といいます。分母 d が 10000 以下の単位分数 1 / d のなかで、循環節が一番長くなる d を求めてください。

本ページでは Python3 (ver 3.8.10) を使ってプログラムを作ることにします。

解答１
解答２
解答３

●解答１

分数を循環小数に直すのは簡単です。筆算と同じように計算していくだけです。次の図を見てください。

     0 1 4 2 8 5 7
   ----------------
 7 ) 1
     0
    ----- 
     1 0  <-- 余り 1
       7
    ------- 
       3 0
       2 8
      -------
         2 0
         1 4
        -------
           6 0
           5 6
          -------
             4 0
             3 5
            -------
               5 0
               4 9
              -----
                 1  <-- 余り 1

7 で割った余り 1 が 2 回現れていますね。これから先は同じ数列を繰り返すことがわかります。7 の剰余は 0 から 6 まであり、剰余が 0 の場合は割り切れるので循環小数にはなりません。現れる余りの数が限られているので、割り切れなければ循環することになるわけです。また、このことから循環節の長さは分母の値よりも小さいことがわかります。

それでは Python (PyPy ver 4.0.1) でプログラムを作ってみましょう。

リスト : 循環小数を求める

# 探索
def position(n, xs):
    for i, x in enumerate(xs):
        if x == n: return i
    return -1

# 循環小数 m/n = ([...],[...])
def repdec(m, n):
    xs = []
    ys = []
    while True:
        p = m // n
        q = m % n
        if q == 0:
            ys.append(p)
            return ys, [0]
        else:
            x = position(q, xs)
            if x >= 0:
                ys.append(p)
                return ys[:x+1], ys[x+1:]
            xs.append(q)
            ys.append(p)
            m = q * 10

関数 repdec(m, n) は m / n を循環小数に変換します。返り値は 2 つのリストを格納したタプルで、先頭のリストが冒頭の小数部分を、次のリストが循環節の部分を表します。途中で割り切れる場合は循環節を [0] とします。これ以降、冒頭の小数部分を有限小数部分と記述することにします。

変数 xs が余りを格納するリスト、変数 ys が商を格納するリストです。最初に m と n の商と余りを計算して、変数 p と q にセットします。q が 0 ならば割り切れたので有限小数です。p を append() で ys に追加して、それと [0] をタプルに格納して返します。

割り切れない場合、余り q が xs にあるか関数 position() でチェックして、その位置を変数 x にセットします。同じ値を見つけた場合、ys の先頭から x 番目までの要素が有限小数部分で、それ以降が循環部になります。p を append() で ys に追加してから、ys を 2 つに分けて返します。見付からない場合は p, q を append() で xs, ys に追加して、m の値を q * 10 に更新してから割り算を続行します。

それでは実行してみましょう。

>>> for x in range(2, 12): print(repdec(1, x))
...
([0, 5], [0])
([0], [3])
([0, 2, 5], [0])
([0, 2], [0])
([0, 1], [6])
([0], [1, 4, 2, 8, 5, 7])
([0, 1, 2, 5], [0])
([0], [1])
([0, 1], [0])
([0], [0, 9])
>>> repdec(355, 113)
([3], [1, 4, 1, 5, 9, 2, 9, 2, 0, 3, 5, 3, 9, 8, 2, 3, 0, 0, 8, 8, 4, 9, 5, 5, 7, 
5, 2, 2, 1, 2, 3, 8, 9, 3, 8, 0, 5, 3, 0, 9, 7, 3, 4, 5, 1, 3, 2, 7, 4, 3, 3, 6, 
2, 8, 3, 1, 8, 5, 8, 4, 0, 7, 0, 7, 9, 6, 4, 6, 0, 1, 7, 6, 9, 9, 1, 1, 5, 0, 4, 
4, 2, 4, 7, 7, 8, 7, 6, 1, 0, 6, 1, 9, 4, 6, 9, 0, 2, 6, 5, 4, 8, 6, 7, 2, 5, 6, 
6, 3, 7, 1, 6, 8])

正常に動作していますね。

●解答２

循環小数を分数に直すことも簡単にできます。循環小数 - Wikipedia によると、有限小数部分を a、循環節の小数表記を b、節の長さを n とすると、循環小数 x は次の式で表すことができる、とのことです。

\(\begin{eqnarray} x &=& a + b\left(1 + \dfrac{1}{10^n} + \dfrac{1}{10^{2n}} + \dfrac{1}{10^{3n}} + \cdots\right) \\ &=& a + b\dfrac{10^n}{10^n - 1} \end{eqnarray}\)

ここで、カッコの中は初項 1, 公比 1 / (10^n) の無限等比級数になるので、その和は次の公式より求めることができます。

\( \displaystyle \sum_{n=0}^{\infty} ar^n = \dfrac{a}{1-r}, \quad (|r| \lt 1) \)

簡単な例を示しましょう。

\(\begin{array}{l} 0.\dot{1}4285\dot{7} = 0.142857 \times \dfrac{10^6}{10^6 - 1} = \dfrac{142857}{999999} = \dfrac{1}{7} \\ 0.1\dot{6} = 0.1 + 0.06 \times \dfrac{10}{9} = \dfrac{1}{10} + \dfrac{6}{100} \times \dfrac{10}{9} = \dfrac{15}{90} = \dfrac{1}{6} \end{array}\)

プログラムを作る場合、次のように考えると簡単です。

有限小数部分を表すリストを xs とすると

有限小数部分 = p0 / q0
ただし p0 = xs を整数に変換
      q0 = 10 ^ (len(xs) - 1)

循環節を表すリストを ys とすると

循環節 = p1 / q1
ただし p1 = ys を整数に変換
      q1 = 10 ^ len(ys) - 1

            p0        p1       p0 * q1 + p1
循環小数 = ---- + --------- = --------------
            q0     q0 * q1        q0 * q1

冒頭の有限小数部分を分数に変換するのは簡単ですね。先頭が整数部分になるので、小数部分の桁はリスト xs の長さ - 1 になります。循環節だけを分数に変換する場合、たとえば 1 / 7 の循環節は [1, 4, 2, 8, 5, 7] になりますが、分子 p' は 0.142857 * 10^6 = 142857 となるので、循環節を表すリストを整数に変換するだけでよいことがわかります。有限小数部分がある場合は、その桁数だけ循環節の部分を右シフトすればよいので、p1 / q1 に 1 / q0 を掛けるだけです。

これを Python でプログラムすると、次のようになります。

リスト : 循環小数を分数に直す

from functools import reduce

# 最大公約数
def gcd(a, b):
    while b > 0:
        a, b = b, a % b
    return a

# 循環小数を分数に直す
def from_repdec(xs, ys):
    # 有限小数の部分を分数に直す
    p0 = reduce(lambda a, x: a * 10 + x, xs, 0)
    q0 = 10 ** (len(xs) - 1)
    # 循環節を分数に直す
    p1 = reduce(lambda a, x: a * 10 + x, ys, 0)
    q1 = (10 ** len(ys) - 1)
    # 有限小数 + 循環節
    a = q0 * q1
    b = q1 * p0 + p1
    c = gcd(a, b)
    return b // c, a // c

# 簡単なテスト
for x in range(2, 12):
    print from_repdec(*repdec(1, x))
print from_repdec(*repdec(355, 113))

アルゴリズムをそのままプログラムしただけなので、とくに難しいところは無いと思います。

それでは実行してみましょう。

>>> for x in range(2, 12): print(from_repdec(*repdec(1, x)))
...
(1, 2)
(1, 3)
(1, 4)
(1, 5)
(1, 6)
(1, 7)
(1, 8)
(1, 9)
(1, 10)
(1, 11)
>>> from_repdec(*repdec(355, 113))
(355, 113)

正常に動作していますね。

●解答３

repdec() を使って単純にプログラムを作ると次のようになります。

リスト : 循環節が一番長い 1 / d を求める

import time

def solver(d):
    k = 0
    m = 0
    while k < d:
        _, xs = repdec(1, d)
        if len(xs) > k:
            k = len(xs)
            m = d
        d -= 1
    return m, k

def test(func, d):
    s = time.time()
    print(func(d))
    print(time.time() - s)

変数 k に循環節の長さの最大値、変数 m にそのときの分母の値 d を保存します。循環節の長さは分母 d よりも小さいので、d の大きな値から循環節の長さを調べていき、d が k 以下になったら、処理を終了して m, k の値を返します。

実行結果は次のようになりました。

>>> test(solver, 10000)
(9967, 9966)
5.3489086627960205

答えは 9967 で、循環節の長さは 9966 になりました。実行時間は 5.35 秒 (実行環境 : Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz) でした。けっこう時間がかかりますね。循環節の長さを求めるだけでよければ、repdec() よりも高速な方法があります。

●プログラムの高速化

今回は小数を 10 進数で表記しているので、d を素因数分解して 2 と 5 の因子しかない場合、1 / d は有限小数になります。d が 2 と 5 の因子を含んでいない場合は循環節だけになります。それ以外の場合は有限小数 + 循環小数の形になります。

簡単な例を示しましょう。

1 / 7  => ([0], [1, 4, 2, 8, 5, 7])
1 / 14 => ([0, 0], [7, 1, 4, 2, 8, 5])
1 / 21 => ([0], [0, 4, 7, 6, 1, 9])
1 / 28 => ([0, 0, 3], [5, 7, 1, 4, 2, 8])
1 / 35 => ([0, 0], [2, 8, 5, 7, 1, 4])
1 / 49 => ([0], [0, 2, 0, 4, 0, 8, 1, 6, 3, 2, 6, 5, 3, 0, 6, 1, 2, 2, 4, 4, 8, 9, 7, 9,
5, 9, 1, 8, 3, 6, 7, 3, 4, 6, 9, 3, 8, 7, 7, 5, 5, 1])
1 / 56 => ([0, 0, 1, 7], [8, 5, 7, 1, 4, 2])
1 / 70 => ([0, 0], [1, 4, 2, 8, 5, 7])

このように、循環節の長さは 2 と 5 以外の因子により決定されます。d に 2 と 5 の因子が含まれていると、循環節の長さは d よりも小さな値になります。そして、1 / d が循環節だけの場合、その長さを n とすると次の式が成り立ちます。

\( 10^n \equiv 1 \pmod d \)

\(\equiv\) は「合同式」を表す記号です。整数 a, b を整数 p で割った余りが等しいとき \(a \equiv b \pmod p\) と書いて、「a と b は p を法として合同である」といいます。

たとえば、10 を 7 で割ると 3 余りますが、24 を 7 で割っても 3 余りますよね。この場合、\(24 \equiv 10 \pmod 7\) と書くことができます。また、3 を 7 で割った余りも 3 なので、\(24 \equiv 3 \pmod 7\) や \(10 \equiv 3 \pmod 7\) と書くこともできます。つまり、式 \(10^n \equiv 1 \pmod d\) は、\(10^n\) を d で割った余りが 1 のとき、循環節の長さが n になることを表しているわけです。

この式は簡単に求めることができます。たとえば、1 / d = 0.(xs)(xs)(xs) ... としましょう。xs は長さ n の循環節を表します。両辺を \(10^n\) して、両辺から 1 / d を引くと次のようになります。

1 / d = 0.(xs)(xs)(xs) ... 
10ⁿ / d = (xs).(xs)(xs) ... 
10ⁿ / d - 1 / d = (xs).(xs)(xs) ... - 0.(xs)(xs) ...
(10ⁿ - 1) / d = xs

10ⁿ - 1 は d で割り切れるので、10ⁿ を d で割れば 1 余るわけです。簡単な例を示しましょう。

1 / 7 => ([0], [1, 4, 2, 8, 5, 7])
(10 ** 1) % 7 = 3
(10 ** 2) % 7 = 2
(10 ** 3) % 7 = 6
(10 ** 4) % 7 = 4
(10 ** 5) % 7 = 5
(10 ** 6) % 7 = 1

1 / 31 => ([0], [0, 3, 2, 2, 5, 8, 0, 6, 4, 5, 1, 6, 1, 2, 9])
(10 ** 1) % 31 = 10
(10 ** 2) % 31 = 7
(10 ** 3) % 31 = 8
(10 ** 4) % 31 = 18
(10 ** 5) % 31 = 25
(10 ** 6) % 31 = 2
(10 ** 7) % 31 = 20
(10 ** 8) % 31 = 14
(10 ** 9) % 31 = 16
(10 ** 10) % 31 = 5
(10 ** 11) % 31 = 19
(10 ** 12) % 31 = 4
(10 ** 13) % 31 = 9
(10 ** 14) % 31 = 28
(10 ** 15) % 31 = 1

確かに 10ⁿ ≡ 1 (mod d) を満たす n が循環節の長さになっています。これをプログラムすると次のようになります。

リスト : 循環節が一番長い 1 / d を求める (2)

# 循環節の長さを求める
def repeat_len(d):
    n = 1
    m = 10 % d
    while m != 1:
        n += 1
        m = (m * 10) % d
    return n

def solver1(d):
    k = 0
    m = 0
    if d % 2 == 0: d -= 1
    while k < d:
        if d % 5 != 0:
            n = repeat_len(d)
            if n > k:
                k = n
                m = d
        d -= 2
    return m, k

プログラムのポイントは循環節の長さを求める関数 repeat_len() です。Python は多倍長整数をサポートしていますが、多倍長整数の計算には時間がかかるので、\(10^n \mod d\) で余りを求めると、実行時間はかえって遅くなってしまいます。この場合、次に示す合同式の性質を使います。

\(a \equiv b \pmod p\) であれば \(a \times c \equiv b \times c \pmod p\)

たとえば、\(10 \equiv 3 \pmod 7\) は両辺を 10 倍すると \(100 \equiv 30 \equiv 2 \pmod 7\) になります。さらに 10 倍すると \(1000 \equiv 20 \equiv 6 \pmod 7\) になります。\(10^k \mod d\) の値を \(m^k\) とすると、\(10^{k+1} \mod d\) の値は \((m^k \times 10) \mod d\) で求めることができます。つまり、\(10^n\) を実際に計算しなくても、\(10^n \mod d\) を求めることができるわけです。repeat_len() は剰余 m が 1 になるまで、この処理を繰り返すだけです。

実行結果を示します。

>>> test(solver1, 10000)
(9967, 9966)
0.0043566226959228516

とても速くなりましたね。ちなみに、d < 100000 の範囲で最長の循環節を求めてみたところ、結果は次のようになりました。

>>> test(solver1, 100000)
(99989, 99988)
0.017964601516723633

答えは 1 / 99989 で、循環節の長さは 99988 になりました。solver1() は 0.018 秒と瞬時に答えを求めることができました。

●参考文献, URL

遠山啓, 『数学入門 (上) (下)』, 岩波新書, 1959
循環小数 - Wikipedia

●プログラムリスト

#
# repdec.py : 循環小数
#
#             Copyright (C) 2016-2022 Makoto Hiroi
#
import time
from functools import reduce

# 探索
def position(n, xs):
    for i, x in enumerate(xs):
        if x == n: return i
    return -1

# 循環小数 m/n = ([...],[...])
def repdec(m, n):
    xs = []
    ys = []
    while True:
        p = m // n
        q = m % n
        if q == 0:
            ys.append(p)
            return ys, [0]
        else:
            x = position(q, xs)
            if x >= 0:
                ys.append(p)
                return ys[:x+1], ys[x+1:]
            xs.append(q)
            ys.append(p)
            m = q * 10

# 最大公約数
def gcd(a, b):
    while b > 0:
        a, b = b, a % b
    return a

# 循環小数を分数に直す
def from_repdec(xs, ys):
    # 有限小数の部分を分数に直す
    p0 = reduce(lambda a, x: a * 10 + x, xs, 0)
    q0 = 10 ** (len(xs) - 1)
    # 循環節を分数に直す
    p1 = reduce(lambda a, x: a * 10 + x, ys, 0)
    q1 = (10 ** len(ys) - 1)
    # 有限小数 + 循環節
    a = q0 * q1
    b = q1 * p0 + p1
    c = gcd(a, b)
    return b // c, a // c

def solver(d):
    k = 0
    m = 0
    while k < d:
        _, xs = repdec(1, d)
        if len(xs) > k:
            k = len(xs)
            m = d
        d -= 1
    return m, k

# 循環節の長さを求める
def repeat_len(d):
    n = 1
    m = 10 % d
    while m != 1:
        n += 1
        m = (m * 10) % d
    return n

def solver1(d):
    k = 0
    m = 0
    if d % 2 == 0: d -= 1
    while k < d:
        if d % 5 != 0:
            n = repeat_len(d)
            if n > k:
                k = n
                m = d
        d -= 2
    return m, k

def test(func, d):
    s = time.time()
    print(func(d))
    print(time.time() - s)

初版 2016 年 1 月 2 日
改訂 2022 年 10 月 22 日