integer sqrt and hypot

整数の平方根と√(x² + y²) と立方根のアルゴリズム

Ｃ言語で書いてあるが、簡単にアセンブリ言語や他言語に移植出来るように書いてあるつもりである。
32bit / 16bit -> 16bit の除算器と 16bit * 16bit -> 32bit の乗算器を持つプロセッサでの利用を想定している。

32bit 整数版 sqrt (平方根)

0x8000 未満の数は 16bit 上位にスケーリングしてから計算している。ニュートン法を行う前に初期値を近似で求め、ニュートン法を２回行い 1 LSB 誤差が得られた後に、誤差を 1/2 LSB に収めるために LSB を調整している。
平方根をニュートン法で求める場合、誤差は必ず負になるので、1 LSB の誤差が得られた時点で数値を増やす方向に調整すれば良い。誤差を 1/2 LSB に納めるためには、与えられた整数を x 、計算結果の整数を y として

y² < x < (y + 1)² … (1)

となる y が存在し、y に対して誤差が最大の 1/2 LSB となるときの x は

x == (y + 1/2)²

であるが、 x は整数なので

x == y² + y

のときである。従って、 y² + y < x の場合に y + 1 を答えにすればよい。
誤差のヒストグラムは 0 ～ 2³¹ の全数検査で次のようになる。

err(+1/4..+1/2 LSB) = 536872070
err(-1/4..+1/4 LSB) = 1073739508
err(-1/2..-1/4 LSB) = 536872070

初期値の近似式は次のように１次近似で計算が簡単になるように係数を調整した。

g(x) = (22 * x + 11) / 32 … (2)

(2) 式で約 3% の精度が得られる。

Fig.1 sqrt(x) の近似誤差 (初期値とニュートン法１回目)
sqrt(x)近似誤差

/* NAME
 *    sqrtl - square root function
 * SYNOPSIS
 *    long
 *    sqrtl(long x)
 * DISCRIPTIONS
 *    The sqrtl() function compute the non-negative square root of x.
 * ERROR
 *    Below 1/2 LSB.
 * SEE ALSO
 *    sqrt(3), http://www.finetune.co.jp/~lyuka/fract/sqrt_hypot.html
 * COPYRIGHT
 *    Copyright 2002, Takayuki HOSODA. All rights reserved.
 */
long
_sqrtl(a)
    long a;
{
    register long x;
    register long t;
    register long s;
    int scale;

    x = a;
    if (x > 0) {
        scale = 0;
        if (x < 0x8000) {
            x <<= 16;
            scale = 8;
            a = x;
        }
        x >>= 8;
        s = 8;
        for (t = 0x400000L; x < t; t >>= 2)
            s--;
        t = 88;
        t <<= s;
        x *= 22;
        s += 5;
        x >>= s;                      /* -3.1e-2 < err < +2.9e-2 */
        s = a;
        t += x;
        x = s;
        s /= t;
        s += t;
        s >>= 1;                      /* -4.8e-4 < err <= 0      */
        t = x;
        x /= s;
        x += s;
        x >>= 1;                      /* -1.2e-7 < err <= 0      */
        s = x;
        s++;
        s *= x;
        if (t > s)                    /* adjust LSB              */
            x++;
        if (scale) {
            x += 127;
            x >>= 8;
        }
    }
    return x;
}

16bit 整数版 hypot

ニュートン法を行う前に初期値を近似で求め、ニュートン法を２回行い 1 LSB 誤差が得られた後に、誤差を 1/2 LSB に収めるために LSB を調整している。
初期値の近似式は、１次近似で 1% 以下の精度が得られる次の式 (3) が適している。

g(x) = (53 * x + 129) / 128 … (3)

しかし、ここでは、得られる精度は約 2% だが計算のより簡単な次の式 (4) を採用した。

h(x) = 13 * x / 32 + 1 … (4)

どちらの場合でも 1 LSB 以内の誤差を得るためには２回のニュートン法が必要であるので大差ないからである。
誤差の様子を Fig.2 に示す。

Fig.2 sqrt(1+x) の近似誤差 (初期値とニュートン法１回目)

sqrt のプログラム例と同様に誤差を 1/2 LSB に納めるために LSB を調整している。

/* NAME
 *    hypots  - euclidean distance function
 * SYNOPSIS
 *    unsigned short
 *    hypots(short x, short y)
 * DESCRIPTION
 *     The hypots() function compute the sqrt(x * x + y * y).
 * ERROR
 *    Below 1/2 LSB
 * SEE ALSO
 *    hypot(3),  http://www.finetune.co.jp/~lyuka/fract/hypot.html
 * COPYRIGHT
 *    Copyright 2002, Takayuki HOSODA. All rights reserved.
 */

unsigned short 
_hypots(x, y)
    short x;
    short y;
{
    register unsigned long a;
    register unsigned long b;
    register unsigned long c;
    a = x;
    b = y;
    if(x < 0) a = -x;
    if(y < 0) b = -y;
    c = b;
    if(a != 0) {
        c = a;
        if (b != 0) {
            if(a < b) {
                a = b; b = c;
            }
            b *= b;
            c = b;
            b /= a;
            b *= 13;
            b >>= 5;
            b += a;    /*       0  < error < +1.9e-2 */
            a *= a;
            c += a;
            a = c;     /* apply newton-raphson methode twice */
            a /= b;
            a += b;
            a >>= 1;   /* -1.7e-4 < error <  0       */
            b = c;
            c /= a;
            c += a;
            c >>= 1;   /* -1.3e-8 < error <  0       */
            a = c;
            a++;
            a *= c;
            if (b > a) /* adjust LSB */
                c++;
        }
    }
    return (unsigned short)c;
}

14bit 整数版 hypot

14bit 精度を得るために、初期値の近似式に式 (3) を使用してニュートン法を１回に減らしている。

/* NAME
 *    hypot14  - euclidean distance function
 * SYNOPSIS
 *    short
 *    hypot14(short x, short y)
 * DESCRIPTION
 *     The hypot14() functions compute the sqrt(x * x + y * y).
 * ERROR
 *     Below 1/2 LSB
 * BUGS
 *     Input x, y must be within range of -11,715 to 11,715 for 1/2 LSB error.
 * SEE ALSO
 *    hypot(3),  http://www.finetune.co.jp/~lyuka/fract/hypot.html
 * COPYRIGHT
 *    Copyright 2002, Takayuki HOSODA. All rights reserved.
 */

#include <math.h>

short 
hypot14(short x, short y) {
    short x;
    short y;
{
    register long a;
    register long b;
    register long c;

    a = (long)abs(x);
    b = (long)abs(y);
    c = b;
    if (a != 0) {
        c = a;
        if (b != 0) {
            if(a < b) {
                a = b; b = c;
            }
            b *= b;
            c = b;
            b /= a;
            b *= 53;
            b += a;
            b >>= 7;
            b += a;    /* |err| < 8.5e-3 */
            a *= a;
            c += a;    /* apply newton-raphson methode */
            a = c;
            c /= b;
            c += b;
            c >>= 1;   /* -7.2e-5 < error < 0, 1/sqrt(2^27) = 8.6e-5 */
            b = c;
            b++;
            b *= c;
            if (a > b) /* adjust LSB */
                c++;
        }
    }
    return c;
}

32bit 整数版 cbrt (立方根)

初期値を近似式(5)で求めたあと、漸化式(6)で計算している。

g(x) = (39 * x + 29) / 64 … (5)

近似式(5)は1次近似で計算が簡単になるように係数を調整した。これで約 6% の精度が得られる。
プログラム例の近似部分では x が 1 の時に近似値が 0 にならないように 1 を加えている。

漸化式(6)は 1次/1次の Pade 近似をもとにしたもので、3次収束を示す。^*Note1

x_n+1 = x_n * (x_n³ + 2 a) / (2 x_n³ + a) … (6)

近似と漸化式1回で 1.5e-4 以下の誤差になり、漸化式2回で 2.2e-12 以下の誤差になるが、
32bit 整数の立方根は高々±1,290 (3桁)程度なので漸化式の計算は1回で十分である。
誤差の様子を Fig.3 に示す。

Fig.3 cbrt(x) の近似誤差

sqrt のプログラム例と同様に誤差を 1/2 LSB に納めるために LSB を調整している。

/* NAME
 *    cbrtl  - cube root function.
 * SYNOPSIS
 *    long
 *    cbrtl(long x)
 * DISCRIPTIONS
 *    The cbrtl() function compute the cube root of x.
 * ERROR
 *    Below 1/2 LSB.
 * SEE ALSO
 *    cbrt(3), http://www.finetune.co.jp/~lyuka/fract/sqrt_hypot.html
 * COPYRIGHT
 *    Copyright 2003, Takayuki HOSODA. All rights reserved.
 */
long
cbrtl(a)
    long a;
{
    long x, t, b;

    if (a == 0)
        return 0;
    if (a == 0x80000000L)
        return -1290;
    x = b = abs(a);
    for (t = 1; x > t; t <<= 1)
        x >>= 2;
    x = ((x * 39 + t * 29) >> 6) + 1;       /* -0.06 < err < 0.06 */
    t = b / (x * x);
    x = (x * (t + t + x)) / (x + x + t);    /* |err| < 1.5e-4 */
    t = x * x;
    if ((6 * t + 3 * x) < (b - x * t) * 4)
        x++;                                /* adjust LSB */
    if (a < 0)
        return -x;
    return x;
}

Note1: 一般に a の n 乗根(但し n≠0) の 3次収束の漸化式は

で与えられる。