Python 有理数のp進展開を計算する

# 実装準備

今回実装したアルゴリズムです。

$\alpha = \displaystyle\frac{1}{15},\ p = 5$

(1) $p$ べきと $p$ の因子を持たない部分に分ける

$\displaystyle \alpha = \frac{1}{3} \cdot 5^{-1}$

(2) $p$ 因子を持たない部分を ${\rm mod}\ p$ で考える

$\displaystyle \frac{1}{3} \equiv 2 \ ({\rm mod} \ 5)$

(3) 式変形をする

$
\displaystyle \alpha = \frac{1}{3} \cdot 5^{-1}
 = 2 \cdot 5^{-1} + \left(\frac{1}{3} - 2 \right) \cdot 5^{-1} \\\\
\displaystyle \ \ = 2 \cdot 5^{-1} + \left(-\frac{5}{3} \right) \cdot 5^{-1}
 = 2 \cdot 5^{-1} + \left(-\frac{1}{3} \right)
$

(4) (1) ～ (3) を繰り返す

$
\displaystyle \alpha = \frac{1}{3} \cdot 5^{-1}
 = 2 \cdot 5^{-1} + \left(\frac{1}{3} - 2 \right) \cdot 5^{-1} \\\\
\displaystyle \ \ = 2 \cdot 5^{-1} + (-\frac{5}{3}) \cdot 5^{-1}
 = 2 \cdot 5^{-1} + \left(-\frac{1}{3} \right) \\\\
\displaystyle \ \ = 2 \cdot 5^{-1} + (-\frac{1}{3})
= 2 \cdot 5^{-1} + 3 \cdot 5^0 + \left(-\frac{1}{3} - 3 \right) \cdot 5^0 \\\\
\displaystyle \ \ = 2 \cdot 5^{-1} + 3 \cdot 5^0 + \left(-\frac{2}{3} \right) \cdot 5^1 
= \cdots
$


# 実装

こちらに置きました。

https://github.com/quryu-1/padic

```python
import math

def get_num_divisible(n, p):
    i = 1
    while n % p ** i == 0:
        i += 1
    return i - 1

def get_inv_elem(n, p):
    if n % p == 0:
        return False
    for i in range(1, p):
        if i * n % p == 1:
            return i

def get_p_adic(frac, p, prec=10):
    if frac[1] == 0:
        return False
    index_dict = {}
    p_pow = get_num_divisible(frac[0], p) - get_num_divisible(frac[1], p)
    i = p_pow
    while i < prec:
        if frac[0] == 0:
            index_dict[i] = 0
            i += 1
            continue
        p_pow = get_num_divisible(frac[0], p) - get_num_divisible(frac[1], p)
        if p_pow == i:
            if p_pow > 0:
                beta = [int(frac[0] / p ** p_pow), frac[1]]
            elif p_pow < 0:
                beta = [frac[0], int(frac[1] / p ** - p_pow)]
            else:
                beta = frac
            x = get_inv_elem(beta[1], p)
            a = beta[0] * x % p
            index_dict[i] = a
            frac = [int((beta[0] - a * beta[1]) * p ** p_pow), beta[1]]
        else:
            index_dict[i] = 0
        i += 1
    return index_dict

def convert_to_formula(index_dict, p):
    if not index_dict:
        return False
    f = ''
    for k, v in index_dict.items():
        if v > 1:
            f += str(v) + '*' + str(p) + '^' + str(k) + ' + '
        elif v == 1:
            f += str(p) + '^' + str(k) + ' + '
        elif v  == 0:
            pass
    f += f'O({p}^{list(index_dict.keys())[-1] + 1})'
    return f

if __name__ == '__main__':
    p = 3
    prec = 10
    frac = [2, 5]
    index_dict = get_p_adic(frac, p, prec)
    f = convert_to_formula(index_dict, p)
    print('p = {}, alpha = {}/{} のp進展開:\n  {}'.format(p, frac[0], frac[1], f))
```

# 実行例

素数 $p$ と 有理数 $\alpha$ について、$\alpha$ の $\mathbb{Q}_p$ における $p$ 進展開を求めてみます。

## 実行例1

$p = 3,\ \alpha = \displaystyle \frac{2}{5}$

```
p = 3, alpha = 2/5 のp進展開:
  3^0 + 3^1 + 2*3^2 + 3^3 + 3^5 + 2*3^6 + 3^7 + 3^9 + O(3^10)
```

**PARI/GP**

```
gp > 2/5 + O(3^10)
%24 = 1 + 3 + 2*3^2 + 3^3 + 3^5 + 2*3^6 + 3^7 + 3^9 + O(3^10)
```

## 実行例2

$p = 5,\ \alpha = \displaystyle \frac{1}{15}$

```
p = 5, alpha = 1/15 のp進展開:
  2*5^-1 + 3*5^0 + 5^1 + 3*5^2 + 5^3 + 3*5^4 + 5^5 + 3*5^6 + 5^7 + 3*5^8 + 5^9 + 3*5^10 + 5^11 + O(5^12)
```

**PARI/GP**

```
gp > 1/15 + O(5^12)
%25 = 2*5^-1 + 3 + 5 + 3*5^2 + 5^3 + 3*5^4 + 5^5 + 3*5^6 + 5^7 + 3*5^8 + 5^9 + 3*5^10 + 5^11 + O(5^12)
```

## 実行例3

$p = 7,\ \alpha = \displaystyle \frac{49}{5}$

```
p = 7, alpha = 49/5 のp進展開:
  3*7^2 + 7^3 + 4*7^4 + 5*7^5 + 2*7^6 + 7^7 + 4*7^8 + 5*7^9 + 2*7^10 + O(7^11)
処理時間: 0.0
```

**PARI/GP**

```
gp > 49/5 + O(7^11)
%26 = 3*7^2 + 7^3 + 4*7^4 + 5*7^5 + 2*7^6 + 7^7 + 4*7^8 + 5*7^9 + 2*7^10 + O(7^11)
```

# PARI/GP

今回比較に使用した PARI/GP は数論の計算ソフトです。

`x + O(p^k)` で p 進数の近似値を求めることができます。

```
gp > 1/15 + O(5^10)
%19 = 2*5^-1 + 3 + 5 + 3*5^2 + 5^3 + 3*5^4 + 5^5 + 3*5^6 + 5^7 + 3*5^8 + 5^9 + O(5^10)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/ac866f78-39bc-89f1-a3e3-6638427b8816.png)

上の画像のように SageMath の CoCalc を使って Python で PARI/GP を使うことができます。

https://www.sagemath.org/

以下の赤枠で囲ったアイコンをクリックします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/b67aa721-dc03-d578-9c75-25f8abf9a596.png)

「Run CoCalc Now」をクリックします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/bb05ce22-beaf-989c-0e32-56f5a4359a93.png)

「SageMath」をクリックします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/d8e93818-bbbc-3b16-0e65-ead3a9099052.png)



# 参考資料

PARI/GPの公式ホームページです。

http://pari.math.u-bordeaux.fr/index.ja.html

こちらのスライドがわかりやすく解説されております。

https://www.sci.u-toyama.ac.jp/~iwao/pari-gp/hokuriku-nt-06-july.pdf

コマンドはこちらのリファレンスにまとめられています。

https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.11.0/refcard.pdf

こちらの記事は非常に細かく、かつ丁寧に整理されています。

https://sehermitage.web.fc2.com/etc/parigp_func.html

Qiita にも丁寧に解説されているのを見つけました。

https://qiita.com/mather314/items/0399668ef3fbbfb4951e

https://qiita.com/iwaokimura/items/a7e71a154f8473a3edeb

こちらにはSageMathについて詳しく書かれています。

https://www.kurims.kyoto-u.ac.jp/~kenkyubu/bessatsu/open/B32/pdf/B32_009.pdf

SageMath のクイックリファレンスです。コマンドが整理されていて見やすいです。

http://math.shinshu-u.ac.jp/~nu/nora/mldoc/ml/1406/0022/content/quickref-algebra-ja.pdf