開平法を用いてn乗根を求める

はじめに

こちらの記事(https://qiita.com/Hisa_/items/57f976f218a24bbf373c )で平方根を開平法で求めるコードを書きました。これを用いてn乗根を求められればなぁと考えていたら、参考になる記事を発見したのでその方法をコードに落としてみることにしました。

手法

手法はこちらの記事(https://wakara.co.jp/mathlog/20201119) を参考にしました。詳細な方法は上の記事を見ていただくとして、こちらの記事で使われているn乗根の求め方について簡単にまとめると、求めたいn乗根を$\sqrt[n]{a} = a^{\frac{1}{n}}$とすると$a^{\frac{1}{n}} = a^{\frac{q}{2^p - 1}}$と書くことができるので、任意のnに応じたp,qを求めて、aをq回かけ平方根をp回とるという操作を何度も繰り返すことで、n乗根の近似値を求めることができるというものです。

改良点

大まかな方法は上の記事の通りですが、コードに実装するにあたっていくつかの点を改良してます。以下がその改良点です。

• 　ｎを素因数分解する　例えばaの15乗根を求めるときに、直接15乗根を求めるのではなく、$ a^{\frac{1}{15}} = (a^{\frac{1}{3}})^{\frac{1}{5}} $　としてaの3数が乗根を求めた後にそれの5乗根を求めるという風にしています。
• 　２のべき乗根を求める部分を別個に設ける　平方根を求める際に、上記のやり方だと対応しきれないのでその部分は別個に設けています。

なお、pとqを求める方法としては　$ 2^p - 1 $を小さい順から求めて、nで割り切れたときのpと$ \frac{2^p - 1}{n} $を採用してます。そして、今回の計算はどこかで計算を打ち切る必要があるのですが、回数は100回繰り返してのちにループを終了するようにしています。

実装コード

from sympy import factorint
from decimal import Decimal, getcontext

def sqrt_by_hand(n: float, digits: int = 100) -> float:
    """
    開平法によって √n を digits 桁（小数点以下）まで正確に求める。
    小数や大きな数値にも安全に対応。
    """
    getcontext().prec = digits * 2  # 精度を十分確保
    n = Decimal(str(n))  # Decimal化で指数表記を防ぐ
    s = format(n, 'f')   # 指数表記を防いで通常文字列化

    if '.' in s:
        int_part, frac_part = s.split('.')
    else:
        int_part, frac_part = s, ''

    if len(int_part) % 2 == 1:
        int_part = '0' + int_part
    if len(frac_part) % 2 == 1:
        frac_part += '0'

    int_pairs = [int(int_part[i:i + 2]) for i in range(0, len(int_part), 2)]
    frac_pairs = [int(frac_part[i:i + 2]) for i in range(0, len(frac_part), 2)]
    frac_pairs += [0] * digits

    pairs = int_pairs + frac_pairs
    result = ""
    remainder = 0
    divisor = 0

    for p in pairs:
        remainder = remainder * 100 + p
        x = 0
        while (divisor * 10 + (x + 1)) * (x + 1) <= remainder:
            x += 1
        remainder -= (divisor * 10 + x) * x
        divisor = divisor * 10 + 2 * x
        result += str(x)

    int_len = len(int_pairs)
    decimal_str = result[:int_len] + "." + result[int_len:int_len + digits]
    return float(decimal_str)
    
def p_q_detect(num):
    p = 0
    q = 0
    while True:
        p += 1
        tmp = 2**p - 1
        if tmp%num == 0:
            q = tmp//num
            return p,q

def calc(p,q,n,num):
    for _ in range(100):
        num *= n**q
        for _ in range(p):
            num = sqrt_by_hand(num)
    return num
    
n = 3
a = 2 # a^(1/n)
ans = 1 #最終的に出力する値

factor = factorint(n)
print(factor)
for k,v in factor.items():
    if k == 2:
        continue#平方根は後で計算するのでスキップ
    else:
        p,q = p_q_detect(k)
        # print(a,b)
        for _ in range(v):
            ans = calc(p,q,a,ans)
            print(ans)
            a = ans
            # print(n)
if ans == 1:#何も計算していない場合 2のべき乗根の時、ここを通る
    ans *= a
if 2 in factor.keys():
    for _ in range(factor[2]):
        ans = sqrt_by_hand(ans,100)
print(ans)

検証

今回は2のn乗（n=2~10)をとりあえず複数個求めてみることにしました。誤差はこちらの記事(https://qiita.com/Hisa_/items/57f976f218a24bbf373c )で使ったものを使用しています

n	ans	誤差
2	1.4142135623730951	2.220446049250313e-16
3	1.2599210498948732	0.0
4	1.189207115002721	1.1102230246251565e-16
5	1.1486983549970349	5.551115123125783e-16
6	1.122462048309373	2.220446049250313e-16
7	1.1040895136738123	5.551115123125783e-16
8	1.0905077326652577	2.220446049250313e-16
9	1.080059738892306	5.551115123125783e-16
10	1.0717734625362931	3.3306690738754696e-16

3の誤差が0.0になってますが確認してみたら誤差が小さくなりすぎて桁落ちを起こしているみたいです。とにかく、かなり精度の良いn乗根を求められているみたいです。

さいごに

やりたいことはできたので概ね満足していますが、このやり方だと一つ欠陥があって、nの素因数を$ \frac{q}{2^p - 1} $で表す必要があるのですが、この式ですべての素数にしかを表せるか分かってません。要はすべてのn乗根を表現できるとは限らないということです。なので、今後はある程度精度を落としても良いので、確実にすべてのnで扱えるようにソースを改良したいと思います。