浮動小数の除算

除算は遅いとか、いろいろ御託は省略する

組み込み系でリアルタイム処理したいからレイテンシが云々とか、そういう御託は省略する。
知りたいことは、 Newton-Raphson 法で、 C 言語で、どうやって浮動小数の逆数を計算するか、だ。
自分もそれが知りたかったし、意外と Google 先生でパッと出てこなかった。 Copilot 君も微妙な回答だった。

基本的な考え方

$x$ の逆数を求めるとして、初期値 $y_0$ を与えて、次の漸化式をループで回せばできる。

y_{t+1}=y_t*(2-x*y_t)

初期値は、できるだけ逆数に近くとる必要がある。変な値だとそもそも収束しない。
初期値を与える為に、複雑な計算や条件分岐を入れるのは意味が全くないのだけれど、初期値として何が適切か書かれた記事を探すのも途中で面倒になった。さっと出てくるなら今回の記事も書かなかった。

初期値を浮動小数のビット分解から作る

どうやってビット分解しても構わないが、折角 C 言語なので、共用体を使う。どうせ Python などを使う人は逆数の計算で悩んだりしないわけだし、 C 言語に特化していたっていいじゃん。

typedef union{
    float data;
    struct{
        unsigned long man:23;
        unsigned long exp:8;
        unsigned long sign:1;
    }bit;
}U_FLOAT_32;

初期値は、次のように与える。

U_FLOAT_32 uX = {fX};
U_FLOAT_32 uY = {1.0f};
uY.bit.man = ~(uX.man);
uY.bit.exp = ~(uX.exp) - 2;

これで何故よいのかは、紙とペンで書いてみれば分かる。私はそうした。
例えば、 $1.9$ の逆数は、 $ca. 0.526$ だけど、仮数部は $1.X$ で正規化するから、仮数部だけ見たら、 $1.9 → ca. 1.052$ となる。同様に、 $1.1 → ca. 1.818$ 、 $1.5 → ca. 1.333$ 。
つまり、逆数をとると $\sqrt{2} = ca. 1.414$ を挟んで反対側へ写る。それをざっくりやるなら、ビット反転で十分だろう、ということである。
指数部については、 7Fh $=127$ をゼロとするわけだから、ビット反転して $-1$ する。後は、仮数部の正規化があるから、更に $-1$ する。

反復の必要回数

確り検証したわけではないけれど、 float 型の場合で $3$ 回、 double 型なら $4$ 回必要。
用途によっては $1$ 回減らしても精度十分。何桁の精度で十分か考えた方がよい。組み込み系で使うならそこまで精度要らなかったりしないかい？
先ほど、初期値の仮数部をざっくりで与えたので、 $2$ 回で十分だったり不十分だったりしているけれど、精度判定の処理を入れるくらいなら多く回した方が早い。

実装例

本当は $0$ の逆数とか、 $1$ の逆数とか、負の数とかあるので、その対応が必要。
そういった入力をしない前提とすれば、コードは短くできる。高速にしたいなら、入力されるはずがないデータを除く処理すら無駄だ。
とりあえず負値だけ例として入れておく。

float calc_inv_float(float fX){
    U_FLOAT_32 uX = {fX};
    U_FLOAT_32 uY;

    // 負値対応の例
    if(uX.bit.sign){
        fX = -fX;
    }

    // 初期値
    uY.bit.man = ~(uX.bit.man);
    uY.bit.exp = ~(uX.bit.exp) - 2;

    // 反復計算
    uY.data = uY.data * ( 2.0f - fX * uY.data );
    uY.data = uY.data * ( 2.0f - fX * uY.data );
    uY.data = uY.data * ( 2.0f - fX * uY.data );

    // 負値対応の例
    if(uX.bit.sign){
        uY.bit.sign = 1;
    }

    return uY.data;
}

double 型の場合

共用体のビット数を以下のように変更する。

typedef union{
    double data;
    struct{
        unsigned long long man:52;
        unsigned long long exp:11;
        unsigned long long sign:1;
    }bit;
}U_FLOAT_64

また、反復を $1$ 回増やす。

double calc_inv_double(double dX){
    U_FLOAT_64 uX = {dX};
    U_FLOAT_64 uY;

    // 負値対応の例
    if(uX.bit.sign){
        dX = -dX;
    }

    // 初期値
    uY.bit.man = ~(uX.bit.man);
    uY.bit.exp = ~(uX.bit.exp) - 2;

    // 反復計算
    uY.data = uY.data * ( 2.0 - dX * uY.data );
    uY.data = uY.data * ( 2.0 - dX * uY.data );
    uY.data = uY.data * ( 2.0 - dX * uY.data );
    uY.data = uY.data * ( 2.0 - dX * uY.data );

    // 負値対応の例
    if(uX.bit.sign){
        uY.bit.sign = 1;
    }

    return uY.data;
}

こういうのを書くと、 C++ だったらもっと汎用性を持って書けるなあ、とか思ってしまうけれど、除算でサブルーチン組むような組み込み系で C++ 使うわけがないね。

ちなみに、最近の組み込み系で使うマイコン、例えば Cortex-M33 だと、 STMicro の仕様書に依れば 32bit 浮動小数の加減乗算 $1$ クロック、除算 $14$ クロックらしい。先のコードはどう見ても $14$ クロック以上使ってしまうから、もっとショボいマイコンじゃないとわざわざ実装する意味がない。