べき乗を実装して既存のモノと比較してみる

はじめに

こんにちは、麻菜結と申します。この記事は実装したべき乗関数と既存のモノを比較してみる記事です。
注意事項を以下に記載します。

計算結果の正確性は保証しません、既存の機能を用いましょう。
エラー処理やコーナーケースの対応は考慮していません、既存の機能を用いましょう。
すでにあるものは自分で再実装する必要はありません、既存の機能を用いましょう。

それではLet's Go!

アルゴリズム

べき乗関数とは実数$x$と$y$があったとき、
$$
f(x, y) = x ^ y
$$
と計算するような関数です。プログラミング言語ではpowなどの名前で実装されていますね。
例えば$3 ^ 2$だったら$9$、$8 ^ {-1}$なら$0.125$、$2 ^ {\frac{1}{2}}$だったら$1.41421356..$、となります。ここでは実数範囲だけを考え、複素数は考えないことにします。
この右肩につく$y$を任意の数を扱いたい場合、単純な複数回の掛け算や逆数だけでは求めることはできません。上にあげたひとよひとよにひとみごろの例はそれで、$\frac{1}{2} = 0.5$のような値が入ってしまうともう少し工夫が必要になります。
まず$y$をこのように分解します。

$$
y = s\frac{n}{d}
$$

ただし、$s$は$1$か$-1$の値をとり、$n$と$d$は整数値を取ります。つまり、実数$y$の分数表現を行います。
このようにすると、べき乗関数がこのように記述することができます。

$$
f(x, y) = ((\sqrt[d]{x})^{n})^{s}
$$

日本語で説明すると、$x$の$\sqrt[d]{x}$を計算し、それを$n$乗して、最後に$y$が負の値だった時に逆数にします。

分数表現はこのサイトを参考にいたしました。
FloatToFraction (JavaScript版)

$\sqrt[d]{x}$はニュートン法によって求めます。良さそうな参考文献をあげられないので、各自ニュートン法　解説で調べると良いと思います。

実装例

プログラムの実装の中ではexponentが該当部分です。powじゃないのは特に意味はないのでお気になさらず。

C言語

#include <stdio.h>
#include <math.h>

#define ABS(x) ((x) < 0 ? -(x) : (x))
#define ERR_CPCY 0.0001
#define MAX_DENOM 100000

typedef struct {
    int numr;
    int denom;
} frac;

double myPow(double x, int k){
    int i;
    double r = 1.0;
    if(k == 0){
        return 1.0;
    }
    for(i = 0; i < k; i++){
        r = r * x;
    }
    return r;
}

frac double2frac(double r){
    int d, n1, n2;
    double v1, v2;
    frac f;
    d = 1;
    while(d < MAX_DENOM){
        n1 = (int)(r * d);
        n2 = n1 + 1;
        v1 = (double)n1 / (double)d;
        v2 = (double)n2 / (double)d;
        if(ABS(v1 - r) < ERR_CPCY){
            f.numr = n1;
            f.denom = d;
            return f;
        }else if(ABS(v2 - r) < ERR_CPCY){
            f.numr = n2;
            f.denom = d;
        }
        d++;
    }
    return f;
}

double root(double k, int r){
    double x, nx;
    if(k == 0.0 || r < 1.0){
        return 0.0;
    }
    if(r == 1){
        return k;
    }
    x = k;
    nx = x - (myPow(x, r) - k) / (r * (myPow(x, r - 1)));
    while(!(ABS(x - nx) < ERR_CPCY)){
        x = nx;
        nx = x - (myPow(x, r) - k) / (r * (myPow(x, r - 1)));
    }
    return nx;
}

double exponent(double x, double y){
    int isMinus;
    frac f;
    isMinus = (y < 0);
    y = ABS(y);
    f = double2frac(y);
    if(x < 1.0 && f.denom != 1){
        return 0.0;
    }
    x = root(x, f.denom);
    x = myPow(x, f.numr);
    x = isMinus ? (1 / x) : x;
    return x;
}

int main(){
    double x = 2.0;
    double y = 0.5;
    
    printf("%.01lf ** %.01lf = %.020lf\n", x, y, exponent(x, y));
    printf("%.01lf ** %.01lf = %.020lf\n", x, y, pow(x, y));
    
    return 0;
}

2.0 ** 0.5 = 1.41421356237468986983
2.0 ** 0.5 = 1.41421356237309514547

12桁目から誤差が生まれています。まあそこそこですね。

実行はこちらからどうぞ
[C] gcc 12.1.0 - Wandbox

Python

ERR_CPCY = 0.0001
MAX_DENOM = 100000


def myPow(x, k):
    r = 1.0
    if(k == 0):
        return 1.0
    for i in range(k):
        r = r * x
    return r

def double2frac(r):
    d, n1, n2 = 1, 0, 0
    v1, v2 = 0.0, 0.0
    numr, denom = 0, 0
    while(d < MAX_DENOM):
        n1 = int(r * d)
        n2 = n1 + 1
        v1 = n1 / d
        v2 = n2 / d
        if(abs(v1 - r) < ERR_CPCY):
            numr = n1;
            denom = d
            return (numr, denom)
        elif(abs(v2 - r) < ERR_CPCY):
            numr = n2
            denom = d
        d = d + 1
    return (numr, denom)

def root(k, r):
    x, nx = 0.0, 0.0
    if(k == 0.0 or r < 1.0):
        return 0.0
    if(r == 1):
        return k
    x = k
    nx = x - (myPow(x, r) - k) / (r * (myPow(x, r - 1)))
    while(not (abs(x - nx) < ERR_CPCY)):
        x = nx
        nx = x - (myPow(x, r) - k) / (r * (myPow(x, r - 1)))
    return nx

def exponent(x, y):
    isMinus = (y < 0)
    y = abs(y)
    numr, denom = double2frac(y)
    if(x < 1.0 and denom != 1):
        return 0.0
    x = root(x, denom)
    x = myPow(x, numr)
    x = 1 / x if isMinus else x
    return x

if __name__ == "__main__":
    x = 2.0
    y = 0.5
    print(f"{x} ** {y} = {exponent(x, y):.020}")
    print(f"{x} ** {y} = {x ** y:.020}")

2.0 ** 0.5 = 1.4142135623746898698
2.0 ** 0.5 = 1.4142135623730951455

こちらも同じく12桁目から誤差が出ています。
内部ではすべて倍精度浮動小数点だと思うのですが、ここまですべての実行結果が違うのはすこし面白いですね。（原因究明とかはしない）

[Python] CPython 3.10.2 - Wandbox

おわりに

お疲れさまでした。数学ライブラリとかが無いようなマイナーな言語でこれが実装出来たので、その時は何の作業も手がつかなくなるくらい興奮しました。
「Cで実装しなおして記事にしてまとめよ♪」とその日はるんるんだったのですが、あとで調べたらほとんどどの処理系にも付属しているmathにpowあることがわかり「私は何をしているんだろう…」と落ち込みましたが、比較記事として書くことにして没ネタを回避しました。
みなさんも、もしべき乗が実装されていない処理系でべき乗が必要になった時には以上のプログラムを参照してみてください （もっと良いアルゴリズムが確実にあるだろうし、そもそも数学ライブラリが無いような環境で難しい算数をするな）
ここまで読んでいただき、ありがとうございました。