自己同形数（Automorphic number）生成のアルゴリズム

【問題 1】10進数における100桁の自己同形数すべての数字和の合計を求めよ

いきなり問題ですが意味が分かりずらいので順番に説明します。

自己同形数（Automorphic number）とは

あまりなじみのない言葉ですが自己同形数(Wikipedia)とは、

平方したとき、下桁の数が自分自身と同じになる数の事である。 例えば $5^2 = 25, 6^2 = 36, 76^2 = 5776,$ そして $890625^2 = 793212890625$, それゆえ 5, 6, 76 , 890625はすべて自己同形数である。

自己同形数（Automorphic number）の生成

英語版のAutomorphic numberにより詳しい説明と以下のSample Programがあります。これをそのまま走らせると10桁までの自己同形数を求めることが出来ます。

基本的には以下の考え方で求めているようです。$k+1$桁の数は$k$桁の自己同形数に桁を加える形になっているので漸化的に求めることが出来ます。

\begin{align}
&k桁の自己同形数をxとすると、\\
&(x^2 - x) \equiv 0 \pmod{10^k} \\

&例えば x=625とするとk=3　\\

&(x^2 - x) = 390625 - 625 = 390000 \equiv 0 \pmod{1000}
\end{align}

def hensels_lemma(polynomial_function, base: int, power: int):
    """Hensel's lemma."""
    if power == 0:
        return [0]
    if power > 0:
        roots = hensels_lemma(polynomial_function, base, power - 1)
    new_roots = []
    for root in roots:
        for i in range(0, base):
            new_i = i * base ** (power - 1) + root
            new_root = polynomial_function(new_i) % pow(base, power)
            if new_root == 0:
                new_roots.append(new_i)
    return new_roots

base = 10
digits = 10

def automorphic_polynomial(x):
    return x ** 2 - x

for i in range(1, digits + 1):
    print(hensels_lemma(automorphic_polynomial, base, i))
#[0, 1, 5, 6]
#[0, 1, 25, 76]
#[0, 1, 625, 376]
#[0, 1, 625, 9376]
#[0, 1, 90625, 9376]
#[0, 1, 890625, 109376]
#[0, 1, 2890625, 7109376]
#[0, 1, 12890625, 87109376]
#[0, 1, 212890625, 787109376]
#[0, 1, 8212890625, 1787109376]

【問題 1】の解答

問題 1 の答えはこのプログラムを使えば以下のように簡単に求まります。（0,1は100桁でないので5,6から派生した自己同形数の合計を求めます）

def digitsum(n):
  return sum(map(int,str(n)))

autopoly = hensels_lemma(automorphic_polynomial, 10, 100)
print(sum(map(digitsum,autopoly[2:])))
#902

【問題 2】12進数における100桁の自己同形数すべての数字和の合計を12進数で求めよ

このプログラムの優れたところは10進数以外でもbaseを変更することで対応できるということです。

自己同形数の値は10進数で出てくるので、数字和の合計をn進数で求めるdsumNとn進数で表示するb10_nを追加で作りました。

dig = ["0","1","2","3","4","5","6","7","8","9","a","b","c","d","e","f"]
def b10_n(x, n):
  ret = ""
  while x > 0:
    ret, x  = dig[x%n]+ret, x//n
  return ret

def dsumN(x,n):
  ret = 0
  while x > 0:
    ret, x = ret+x%n, x//n
  return ret

N = 100
base = 12
autopoly = hensels_lemma(automorphic_polynomial, base, N)
print(b10_n(sum(map(lambda x:dsumN(x,base),autopoly[2:])),base))
# 77a

【問題 3】12進数における10000桁の自己同形数すべての数字和の合計を12進数で求めよ

桁数がこれだけ大きくなると上記のアルゴリズムでは無理になります。そこでWikipediaをよく見ると以下のアルゴリズムが見つかりました。これにより桁数が倍々となるので$O(log(n))$で生成することが出来ます。下の例では10進数で示していますがn進数でも適用できます。

\begin{align}
&n > 1 の k 桁の自己同形数が与えられたとき, 　\\
&高々 2k 桁の自己同形数 n' を発見することができる \\

&{\displaystyle n'=(3\cdot n^{2}-2\cdot n^{3})\ {\bmod {10^{2k}}}}
\end{align}

これをそのままコードにした関数amorphです。最初の1桁（12進数の場合は4と9）を渡すと必要な桁数になるまで倍々にして行きます。このコードでは$4$で初めて50桁の12進数の自己同形数を生成しています。

def amorph(n, N, base):
  k = 1
  while k <= N:
    n1 = (3*n**2-2*n**3) % base**(2*k)
    k, n = 2*k, n1
  return n % base**N

base, fp, N = 12, 4, 50
print(b10_n(amorph(fp,N,base),base))
# 3599830402b6468a133a81a3a16986b267b3452b21b61b3854

これを利用して、最初の1桁の（0,1以外の）数を見つける関数fplと、それをamorphに入力してその数字和を求めるamorphDsumを作れば解答が出ます。

def fpl(base):
  return [i for i in range(2,base) if i == (i**2)%base]

def amorphDsum(N,base):
  return b10_n(sum([dsumN(amorph(fp,N,base), base) for fp in fpl(base)] ),base)

base, N = 12, 50
print(amorphDsum(N,base))
# 3a0

このアルゴリズムはEuler Project Problem 284: Steady Squaresを解くのに役に立ちます。

（開発環境：Google Colab）