

Google Colaboratoryを使っています。

相異な自然数の逆数の和が１  
`1/x_1 + 1/x_2 +‥‥+1/x_n = 1 (x_1>x_2>‥‥>x_n)`  
自然数ｎのそれぞれの値に対する(x_1〜x_n)の組の個数や各組の`ｘ`について調べたいと思っています。

色々調べてみましたが逆数ゆえかエラーばかり吐かれてしまいます。  
コードの例など教えてほしいです。

また相異な`ｘ`をa,b‥,ｚとして導入していましたが上記のように（x_1>x_2>‥‥>x_n)とするなど、うまい書き方は可能でしょうか。

#### 試した例１

エラーは出ませんでしたが、すべての解が求められていません

[code]

    for i in range(2, 999): 
     for j in range(i+1, 999):
      for k in range(j+1, 999):
       for l in range(k+1, 999):
         for m in range(l+1, 999):
          if 1/i + 1/j + 1/k + 1/l + 1/m == 1:
           print(i, j ,k ,l, m)
    
[/code]

#### 例２

無効な構文と表示されます。

[code]

    import sympy as sym
    from sympy.abc import a, b, c, d, e, n 
    from sympy.solvers.diophantine import diophantine
      
    f = sym.Eq(a**-1 + b**-1 + c**-1 + d**-1 + e**-1, n)
    diophantine(f.subs(n, 1))
    
[/code]


CC BY-SA 4.0
2022-11-09T10:37:18.743
92108
2022-11-14T04:34:56.953
2022-11-09T12:17:46.830
3060
55291
question
python google-colaboratory
逆数和の不定方程式を解きたい


<https://docs.python.org/ja/3/library/fractions.html>  
python には有理数クラスってのがはじめからあるので  
これを使うと無理数のでてこない計算はエラー誤差を気にせず正確にできます

[code]

    from fractions import Fraction
    
    n = 5
    results = []
    
    # x_list のうち fix_index までを固定して残りの変数を変化させて検索
    # rest_value はまだ決まってない項で使える合計値(有理数インスタンスを渡す)
    def iter(x_list, fix_index, rest_value):
        if fix_index >= n: # 全部の変数が決まった
            if rest_value == 0: # １とぴったり一致した
                results.append(x_list[:]) # 副作用を受けないようにコピーして結果に追加
            return
    
        # 残りで使える値がすでにマイナスだったら検索終了
        if rest_value <= 0:
            return
    
        # 次の変数を決めて探索する
        # e.g. x1 x2 を固定 | x3 x4 x5 を変化させて探す場合を考える
        # x1 < x2 < x3 ... を仮定 (質問とは逆順だけど最大値不明なので 2からはじめてだんだん増やしていきたいから低い順に)
        # なので検索範囲最小値は x2 + 1 (x1 は１つ前がないので 2 から開始)
        min_x = x_list[fix_index - 1] + 1 if fix_index else 2
        # 最大値は 1/x3 + 1/x3 + 1/x3 > rest_value となる最大の x3
        # なぜなら x3 < x4 < x5 に対して 1/x3 + 1/x3 + 1/x3 > 1/x3 + 1/x4 + 1/x5 = rest_value
        # つまり x3 < (rest_value / 決まってない変数の数)の逆数 = 決まってない変数の数 / rest_value
        max_x = int((n - fix_index) / rest_value)
        for x in range(min_x, max_x + 1) :
            # x3 を決定して再度残りを検索する
            x_list[fix_index] = x
            iter(x_list, fix_index + 1, rest_value - Fraction(1,x)) # 1/x の有理数インスタンス
    
    
    def main():
        x_list = [0] * n;
        iter(x_list, 0, Fraction(1)) # 最初は固定なしで全部の合計が 1 のものを検索
        print(results)
    
    
    if __name__ == '__main__':
        main()
    
[/code]

実行結果

[code]

    > python test.py
    [[2, 3, 7, 43, 1806], [2, 3, 7, 44, 924], [2, 3, 7, 45, 630], [2, 3, 7, 46, 483], [2, 3, 7, 48, 336], [2, 3, 7, 49, 294], [2, 3, 7, 51, 238], [2, 3, 7, 54, 189], [2, 3, 7, 56, 168], [2, 3, 7, 60, 140], [2, 3, 7, 63, 126], [2, 3, 7, 70, 105], [2, 3, 7, 78, 91], [2, 3, 8, 25, 600], [2, 3, 8, 26, 312], [2, 3, 8, 27, 216], [2, 3, 8, 28, 168], [2, 3, 8, 30, 120], [2, 3, 8, 32, 96], [2, 3, 8, 33, 88], [2, 3, 8, 36, 72], [2, 3, 8, 40, 60], [2, 3, 8, 42, 56], [2, 3, 9, 19, 342], [2, 3, 9, 20, 180], [2, 3, 9, 21, 126], [2, 3, 9, 22, 99], [2, 3, 9, 24, 72], [2, 3, 9, 27, 54], [2, 3, 9, 30, 45], [2, 3, 10, 16, 240], [2, 3, 10, 18, 90], [2, 3, 10, 20, 60], [2, 3, 10, 24, 40], [2, 3, 11, 14, 231], [2, 3, 11, 15, 110], [2, 3, 11, 22, 33], [2, 3, 12, 13, 156], [2, 3, 12, 14, 84], [2, 3, 12, 15, 60], [2, 3, 12, 16, 48], [2, 3, 12, 18, 36], [2, 3, 12, 20, 30], [2, 3, 12, 21, 28], [2, 3, 14, 15, 35], [2, 4, 5, 21, 420], [2, 4, 5, 22, 220], [2, 4, 5, 24, 120], [2, 4, 5, 25, 100], [2, 4, 5, 28, 70], [2, 4, 5, 30, 60], [2, 4, 5, 36, 45], [2, 4, 6, 13, 156], [2, 4, 6, 14, 84], [2, 4, 6, 15, 60], [2, 4, 6, 16, 48], [2, 4, 6, 18, 36], [2, 4, 6, 20, 30], [2, 4, 6, 21, 28], [2, 4, 7, 10, 140], [2, 4, 7, 12, 42], [2, 4, 7, 14, 28], [2, 4, 8, 9, 72], [2, 4, 8, 10, 40], [2, 4, 8, 12, 24], [2, 4, 9, 12, 18], [2, 4, 10, 12, 15], [2, 5, 6, 8, 120], [2, 5, 6, 9, 45], [2, 5, 6, 10, 30], [2, 5, 6, 12, 20], [3, 4, 5, 6, 20]]
    
[/code]


CC BY-SA 4.0
2022-11-09T15:24:33.027
92115
2022-11-11T18:59:32.240
2022-11-11T18:59:32.240
92108
answer


条件判定式の両辺に `i * j * k * l * m` を掛けた式で判定することにしました。整数なので精度の考慮は必要ありません。なお，

[code]

    i < j < k < l < m
    1/5 + 1/5 + 1/5 + 1/5 + 1/5 = 1
    1/2 + 1/8 + 1/8 + 1/8 + 1/8 = 1
    1/2 + 1/3 + 1/18 + 1/18 + 1/18 = 1
    
[/code]

より，`i < 5, j < 8, k < 18` が成り立ちます。  
これを使って計算時間を減らしましたが，うわさどおり（?） Python の for ループは遅いようで，試す場合は `stop = 100`
から始めて徐々に増やすことをお勧めします。

[code]

    stop = 100             # range() does not include 'stop'
    # stop = 1000
    nv = 5                 # number of variables(i, j, k, l, m)
    stop_i = nv            # 1/i * nv > 1 -> i < nv
    stop_j = (nv - 1) * 2  # 1/2 + 1/j * (nv - 1) > 1 -> j < (nv - 1) * 2
    stop_k = (nv - 2) * 6  # 1/2 + 1/3 + 1/k * (nv - 2) > 1 -> k < (nv - 2) * 6
    
    for i in range(2, stop_i):
        for j in range(i+1, stop_j):
            for k in range(j+1, stop_k):
                for l in range(k+1, stop):
                    for m in range(l+1, stop):
                        if (j * k * l * m + i * k * l * m
                            + i * j * l * m + i * j * k * m
                            + i * j * k * l == i * j * k * l * m):
                            print(f'{i:4d} {j:4d} {k:4d} {l:4d} {m:4d}')
    
[/code]

[code]

       2    3    7   78   91
       2    3    8   32   96
       2    3    8   33   88
       2    3    8   36   72
       2    3    8   40   60
       2    3    8   42   56
       2    3    9   22   99
       2    3    9   24   72
       2    3    9   27   54
       2    3    9   30   45
       2    3   10   18   90
       2    3   10   20   60
       2    3   10   24   40
       2    3   11   22   33
       2    3   12   14   84
       2    3   12   15   60
       2    3   12   16   48
       2    3   12   18   36
       2    3   12   20   30
       2    3   12   21   28
       2    3   14   15   35
       2    4    5   28   70
       2    4    5   30   60
       2    4    5   36   45
       2    4    6   14   84
       2    4    6   15   60
       2    4    6   16   48
       2    4    6   18   36
       2    4    6   20   30
       2    4    6   21   28
       2    4    7   12   42
       2    4    7   14   28
       2    4    8    9   72
       2    4    8   10   40
       2    4    8   12   24
       2    4    9   12   18
       2    4   10   12   15
       2    5    6    9   45
       2    5    6   10   30
       2    5    6   12   20
       3    4    5    6   20
    
[/code]

（追記，御参考）  
**numba** モジュールの導入（`pip install numba`）により計算時間の短縮を試みました。私の環境（macOS 13(M1),
Python 3.10.8, numba 0.56.4）では，`stop = 1000` で約40分の1に短縮できました。

[code]

    from numba import jit
    
    @jit(nopython=True)
    def rp_search(stp_i, stp_j, stp_k, stp_l, stp_m):
        rslt = []
        for i in range(2, stp_i):
            for j in range(i+1, stp_j):
                for k in range(j+1, stp_k):
                    for l in range(k+1, stp_l):
                        for m in range(l+1, stp_m):
                            if (j * k * l * m + i * k * l * m
                                + i * j * l * m + i * j * k * m
                                + i * j * k * l == i * j * k * l * m):
                                rslt.append([i, j, k, l, m])
        return rslt
    
    
    stop = 1000
    nv = 5                 # number of variables
    stop_i = nv
    stop_j = (nv - 1) * 2
    stop_k = (nv - 2) * 6
    
    result = rp_search(stop_i, stop_j, stop_k, stop, stop)
    
    for r in result:
        print(f'{r[0]:4d} {r[1]:4d} {r[2]:4d} {r[3]:4d} {r[4]:4d}')
    
[/code]


CC BY-SA 4.0
2022-11-11T09:29:43.693
92160
2022-11-14T04:34:56.953
2022-11-14T04:34:56.953
54588
54588
92108
answer