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

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

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

### 自己同形数（Automorphic number）とは

あまりなじみのない言葉ですが[自己同形数(Wikipedia)](https://ja.wikipedia.org/wiki/%E8%87%AA%E5%B7%B1%E5%90%8C%E5%BD%A2%E6%95%B0)とは、

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


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

英語版の[Automorphic number](https://en.wikipedia.org/wiki/Automorphic_number)により詳しい説明と以下のSample Programがあります。これをそのまま走らせると10桁までの自己同形数を求めることが出来ます。

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

```math 
\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}
```


```python 
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から派生した自己同形数の合計を求めます）

```python 
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**を追加で作りました。

```python 
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進数でも適用できます。

```math
\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進数の自己同形数を生成しています。

```python 
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**を作れば解答が出ます。

```python 
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）