実装準備

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

$\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
$

実装

こちらに置きました。

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)

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

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

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

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

参考資料

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

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

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

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

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

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

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