平方根を用いた立方根の近似と実装例

for HP 42S and C99
平方根を用いた立方根の近似と実装例

2003 年 3 月 4 日
2025 年 4 月 30 日改訂
細田隆之
English

ここでは、平方根のみが利用可能な環境を想定した、実数および複素数に対する立方根の近似手法を紹介します。実数版では、区間 x ∈ [0.125, 8] に最適化された 3 項、あるいは推奨される 4 項の Puiseux 的近似式を初期値として用い、Padé 型の 3次収束の漸化式によって高精度な結果を得ます。複素数版では、より広い値域に対応するため、複数の平方根を組み合わせて x^{m / 2ⁿ} の形での指数近似による初期値生成を行い、収束の安定性を確保しています。

平方根による立方根の近似

x の平方根は x^1/2 で与えられ、 x の立方根は x^1/3 で与えられます。立方根の指数 1/3 を m / 2ⁿ | m = round(2ⁿ / 3）で近似する場合、近似誤差は次式で与えられます。
m / 2ⁿ × 3 - 1

その誤差は ±1/2n（負号は n が偶数の場合）となる。具体的な値を次に示します。
  1 /   2 × 3 - 1 = +0.5
  1 /   4 × 3 - 1 = -0.25
  3 /   8 × 3 - 1 = +0.125
  5 /  16 × 3 - 1 = -0.0625
 11 /  32 × 3 - 1 = +0.03125
 21 /  64 × 3 - 1 = -0.015625
 43 / 128 × 3 - 1 = +0.0078125
 85 / 256 × 3 - 1 = -0.00390625
Fig.1: 立方根 x ^1/3 を x ^m/2ⁿ で近似した場合の近似誤差
これらのうち、[1e-10, 1e+10] の範囲において、初期値の近似誤差が半桁、相対誤差で言えば ±1/√10 ≈ ±0.316 を達成するのに x^11/32 が好都合なのがわかりました。

$x_0 &=& \sqrt{ \sqrt{a \sqrt{ \sqrt{a \sqrt{a} } } } } = a^{11/32} = a^{0.34375}$

初期近似式 (1) の最大近似誤差は x ∈ [1e-12:1e+12] において、x = 1e+12 のときに約 +0.366、すなわち約 0.437 桁の精度でしたので、後述の 3次収束の漸化式 (6) の 4回の反復によって精度が 3⁴ 倍すなわち 81 倍向上し、 0.437 × 3⁴ ≈ 35.3 桁 > 16 桁となって IEEE 754 の倍精度に到達できることになります。

漸化式による立方根の計算

式 (1) で得られた近似値を初期値として、ニュートン法のような漸化式により立方根を求めます。
$x = \sqrt[3]{a}$
相対誤差 h の次のように定義します。
$h = 1 - \frac{a}{x^3}$
立方根関数 x^1/3 は x ≠ 0 において無限級数に展開できるので、 (1 + h)^1/3 は次の級数展開で表せます。
$(1 + h)^{1/3} \approx 1 + \frac{1}{3}h - \frac{1}{9} h^2 + \frac{5}{81} h^3 - \frac{10}{243} h^4 + \frac{22}{729} h^5 + O(h^6)$
より効率的な反復式を得るために、この級数を (1, 1) パデ近似 (Padé approximation) に置き換え、1次と 2次の係数が級数展開の係数と一致するようにします。
$(1 + h)^{1/3} \approx \frac{2 h + 3}{h + 3}$
式(3)と式(5)から、以下の漸化式が導かれます。
$x_{i+1} &\leftarrow& x_i\left(\displaystyle\frac{x^3_i + 2a}{2x^3_i + a}\right)$

解の近傍では、ニュートン・ラフソン法は 2次の収束を示しますが、 (1, 1) パデ近似による漸化式は 2次の項まで級数展開と一致するため、この漸化式は 3次の収束を示します。
また、ニュートン・ラフソン法では、初期推定値が真値から離れている場合に不安定になることがありますが、このパデ近似による漸化式では、そのような場合にも優れた数値的安定性を示します。導出過程は異なりますが、この漸化式は Halley's method によるものと等価です。

この方法は a の n乗根 (ただし a ≠ 0) にも拡張でき、次の 3次収束の漸化式が導かれます。

また、パデ近似の次数を (2, 2) または (3, 3) にした場合にも、式 (6) と同様に、以下の漸化式 (8)、(9) が導かれます。それぞれの漸化式は、5次および 7次の収束を示します。
$x_{i+1} &\leftarrow& x_{i}\left(\frac{5x_i^6 + 35ax_i^3 + 14a^2}{14x_i^6 + 35 a x_i^3 + 5a^2}\right)$

$x_{i+1} &\leftarrow& x_{i}\left(\frac{2x^9 + 30ax^6 + 42a^2x^3 + 7a^3}{ 7x^9 + 42ax^6 + 30a^2x^3 + 2a^3}\right)$

初期近似誤差と漸近計算の収束性

Fig.2: cbrt 初期近似誤差とその収束挙動 (IEEE 754倍精度環境での計算）

複素数の立方根の実装例 (HP 42S 及び C99 ソースコード)

使用計算式

初期近似：
$x_0 &=& \sqrt{ \sqrt{a \sqrt{ \sqrt{a \sqrt{a} } } } } = a^{11/32} = a^{0.34375}$

3次収束の漸化式：
$x_\mathrm{n+1} &\leftarrow& x_\mathrm{n} \cdot \frac{2a + x_\mathrm{n}^3}{a + 2x_\mathrm{n}^3}$

cbrt - for HP 42S
Rev.1.01 (March 22, 2025) ccbrt - complex cube root function for C99

00 { 66-Byte Prgm }
01▸LBL "cbrt"
02 ENTER
03 X=0?
04 RTN
05 ENTER
06 ENTER
07 SQRT
08 ×
09 SQRT
10 SQRT
11 ×
12 SQRT
13 SQRT
14▸LBL 00
15 ENTER
16 ENTER
17 X↑2
18 ×
19 ENTER
20 +
21 LASTX
22 RCL+ ST T
23 RCL+ ST T
24 X<>Y
25 RCL+ ST T
26 ÷
27 X<>Y
28 ×
29 LASTX
30 RCL÷ ST Y
31 1
32 -
33 ABS
34 1E-11
35 X>Y?
36 GTO 01
37 R↓
38 R↓
39 GTO 00
40▸LBL 01
41 R↓
42 R↓
43 RTN
44 .END.

/* REFERNCE ONLY:
   FOR EXPLANATION OF THE ALGORITHM AND NOT INTENDED FOR ACTUAL USE.*/

#include <math.h>
#include <float.h>
#include <complex.h>

double complex ccbrt(double complex x);
/* The ccbrt(x) function compute the cube root of x. */

/* ccbrt Rev.1.43 (April 10, 2025) (c) 2003 Takayuki HOSODA */
double complex ccbrt(double complex a) {
    double complex x, xp, c;
    double         r, rp, s;
    if (a == 0) return 0;
    s = 1.0;
    //Simple range reduction
    while (cabs(a) < 0x1.0p-195) { a *= 0x1.0p+195; s *= 0x1.0p-65; }
    while (cabs(a) > 0x1.0p+195) { a *= 0x1.0p-195; s *= 0x1.0p+65; }
    while (cabs(a) < 0x1.0p-39)  { a *= 0x1.0p+39;  s *= 0x1.0p-13; }
    while (cabs(a) > 0x1.0p+39)  { a *= 0x1.0p-39;  s *= 0x1.0p+13; }
    //Initial approximation by a^(11/32)
    x = csqrt(csqrt(a*csqrt(csqrt(a*csqrt(a)))));
    r = 0x1.0p+39;
    do { //Asymptotic formula
        xp = x; rp = r;
        c = x * x * x;
        #if defined( USE_SEVENTH )
        x += x * (a - c) * (5 * (a * a  + c * c) + 17 * a * c )
            / (a * a * (a + a + 30 * c) + c * c * (42 * a + 7 * c));
        #elif defined( USE_FIFTH )
        x += 9 * x * (a + c) * (a - c) 
            / (14 * c * c + 35 * a * c + 5 * a * a);
        #else
        x += x * (a - c) / (c + c + a); //Cubic convergence
        #endif
        r = cabs(x - xp); //Correction radious
    } while (r != 0 && rp > r); //Convergence check by Urabe's method
    return s * x;
}

注：実用には例外処理とレンジ・リダクションが必要かも知れません。

ダウンロード : cbrt.raw — (Free42 用 raw プログラムファイル)

実数の立方根の実装例 (ANSI C ソースコード)

実数の範囲において sqrt 関数の他に、実数の指数部を直接扱う関数 ldexp, frexp が使える環境の場合には、レンジ・リダクションが簡単になります。
平方根、すなわち 1/2 のべきを用いた ピュイズー級数 (Puiseux series) 的な近似を用い、レンジ・リダクションによる [0.125, 8] の範囲において最適化した 3 項の初期近似式 (10) を下に示します。

$x_0 &=& 0.20493 + 0.88852 \sqrt{a} - 0.09345 a, \quad \ x \in [0.125, 8]$

初期近似式 (10) の最大近似誤差は x ≈ 3.065 のときに約 0.0148、すなわち約 1.83 桁の精度でしたので、前述の 3次収束の漸化式 (6) を 2回適用すれば 1.83 × 3² ≈ 16.46 桁 となって IEEE 754 の倍精度に到達できることになります。

Fig.3: 立方根を平方根を含む 3～5 項の級数で近似した場合の近似誤差 (x ∈ [0.125, 8])

もし、式 (10) を x87 の 80-bit 拡張倍精度での計算に適用しようとすると、3次収束の漸化式 2回では拡張倍精度の 19 桁の精度には到達できません。その場合には次の [0.125, 8] の範囲において最適化した 4 項の初期近似式 (11) では、最大誤差が約 0.0049 すなわち約 2.31 桁の精度でしたので、同様に 2.31 × 3² = 20.79 桁 となって拡張倍精度に到達できることになります。

$x_0 &=& c_1 + c_2 a^{1/2} + c_3 a + c_4 a^{3/2}, \quad x \in [0.125, 8] &(11)\\ \text{where:}&& c_1= 0.166466,\ c_2 = 1.025738, c_3 = - 0.22500,\ c_4 = 0.032796$

_cbrt - cube root function for C

#include <math.h> #include <float.h> double _cbrt(double d); /* The _cbrt(d) function compute the cube root of d.*/ /* _cbrt Rev.2.0 (April 23, 2025) (c) 2003 Takayuki HOSODA */ double _cbrt(double d) { double r, x, a, c; int e; if (d == 0) return 0; /* Range reduction */ r = frexp(fabs(d), &e); x = ldexp(r, e % 3); /* Initial approximation for 0.125 <= x <= 8, |err| < 0.0049 */ a = 0.1664647 + (1.025745 + 0.032798 * x) * sqrt(x) - 0.2250075 * x; /* Apply the asymptotic formula of cubic convergence twice */ c = a * a * a; a += a * (x - c) / (c + c + x); c = a * a * a; a += a * (x - c) / (c + c + x); /* A 1-ulp correction to ensure result consistency. */ r = fabs(x); if (a * a * a * (1.0 - DBL_EPSILON) > r) a *= 1.0 - 0.5 * DBL_EPSILON ; if (a * a * a * (1.0 + 2.0 * DBL_EPSILON) < r) a *= 1.0 + DBL_EPSILON ; /* Range expansion */ return ldexp(d > 0.0 ? a : -a, e / 3); }

_cbrtl - cube root function for C

#include <math.h> #include <float.h> long double _cbrtl(long double d); /* The _cbrtl(d) function compute the cube root of d.*/ /* _cbrtl Rev.2.0 (April 23, 2025) (c) 2003 Takayuki HOSODA */ long double _cbrtl(long double d) { long double r, x, a, c; int e; if (d == 0) return 0; /* Range reduction */ r = frexpl(fabsl(d), &e); x = ldexpl(r, e % 3); /* Initial approximation for 0.125 <= x <= 8, |err| < 0.0049 */ a = 0.1664647 + (1.025745 + 0.032798 * x) * sqrtl(x) - 0.2250075 * x; /* Apply the asymptotic formula of cubic convergence twice */ c = a * a * a; a += a * (x - c) / (c + c + x); c = a * a * a; a += a * (x - c) / (c + c + x); /* A 1-ulp correction to ensure result consistency. */ r = fabsl(x); if (a * a * a * (1.0 - LDBL_EPSILON) > r) a *= 1.0 - 0.5 * LDBL_EPSILON ; if (a * a * a * (1.0 + 2.0 * LDBL_EPSILON) < r) a *= 1.0 + LDBL_EPSILON ; /* Range expansion */ return ldexpl(d > 0.0 ? a : -a, e / 3); }

注：

GCC 7.5 を用いた実験では、_cbrt および _cbrtl 関数に対して、1,000,000 個のランダムな倍精度・拡張倍精度の値を用いて評価を行いました。 相対誤差 _cbrt(d)³ / d - 1 および _cbrtl(d)³ / d - 1 は、それぞれ ±4.44089209850063e-16 および ±2.16840434497101e-19 以内に収まり、いずれも IEEE 754 におけるマシン・イプシロンの 2 倍以内でした。

いずれの関数も負の引数に対応しており、d < 0 の場合には、それぞれ -_cbrt(-d), -_cbrtl(-d) を返します。これは、立方根が奇数乗根であるため、実数全体で連続に定義できるという性質に基づいています。

実数の立方根の初期近似式

実数立方根の実装例では、初期近似として平方根の項を含む近似式 (10) を用いました。この近似式は以下の考えに基づいて設計されています。
x ≈ 1 付近における x² の立方根の対して次数 (2, 0) の Padé 近似を行います。
その結果得られた式において x を √x で置き換えることで、 x 自身の立方根に対する近似式が得られます。
このとき、近似対象となる範囲も二乗されるため、立方根に対する近似誤差が少ない領域が拡がる効果が生じます。
得られた近似式をもとに、x = 1 で 1 になる条件を維持しつつ、 x ∈ [0.125, 8] における近似誤差の絶対値が最小になるよう、係数の最適化を行いました。
この結果が近似式 (10) です（Fig.4 参照）。

Fig.4: 立方根の二乗領域での近似 (x ∈ [0.125, 8])

その他の立方根初期近似式例

参考として、同様の考えに基づいて設計した、より高次の近似式とその比較グラフを以下に示します。これらは、平方根項を含む 4～5 項程度の式によって立方根を近似した例です。
近似式 :

Fig.5: 立方根を平方根を含む 4～5 項の展開式で近似した場合の近似誤差 (x ∈ [0.125, 8])

まとめ

本稿では、平方根のみを用いて立方根を近似する手法と、その実装例について考察しました。初期近似には、実数版では Puiseux 級数に類似した形に基づく近似式を、複素数版では x ^m/2ⁿ 型の指数近似用い、続いて Padé 型の反復法により精度の向上を図る構成を採りました。

提案する近似式および反復法は、比較的簡便な形でありながら、実数および複素数のいずれにおいても良好な収束性と精度を示すことが確認されました。また、平方根の計算コストは除算と同程度であるため、本稿で紹介したアプローチは、広く一般的なものではないかもしれませんが、立方根計算における実用的な選択肢の一つになり得ると考えています。

今後、より高次の初期近似や他の反復法との比較検討により、さらに高精度または高速な手法が得られる可能性も期待されます。

※本稿で述べる、平方根を用いた立方根の初期近似に関する手法 — すなわち、 x ^2/3 の Padé 近似を x = 1 近傍で構成し、その後 x を √x に置き換えることで x ^1/3 の関数に変換し、結果として有効な近似範囲を拡大するという方法 — は、 Puiseux 近似の着想に基づいて著者が独自に構成した手法であり、執筆時点においては同様の例を文献上に確認できておりません。

参考文献

"Convergence of Numerical Iteration in Solution of Equations.", Urabe Minoru, J. Sci. Hiroshima Univ. Ser. A 19 (1956), no. 3, 479--489. doi:10.32917/hmj/1556071264.
"cbrt.c", GNU C Library — (accessed April 11, 2025, version 2.41) glibc-2.41
Jean-Michel Muller, "Elementary Functions: Algorithms and Implementation", Birkhauser, 3rd Edition, 2016.
Okumura, Haruhiko; Shudo, Kazuyuki; Sugiura, Masaki; Tsuchimura, Nobuyuki; Tsuru, Kazuo; Hosoda, Takayuki; Matsui, Yoshimitsu; & Mitsunari, Shigeo. (2003). Java ni yoru Arugorizumu Jiten [Encyclopedia of Algorithms with Java]. Gijutsu Hyoron Sha. ISBN4-7741-1729-3 [in Japanese]

追補

cbrt の使用法 (HP 42S)

操作

[XEQ] "cbrt"   — cube root function

入力

REG X      : argument value a

出力

REG y      : a
REG X      : cbrt(a)

cbrt の計算例 (HP 42S)

 2 "cbrt" → 1.25992104989
-3 "cbrt" → 0.72112478515 + 1.24902476648 i
0 ENTER 1 COMPLEX "cbrt" → 0.866025403784 + 0.5 i
0 ENTER 10 COMPLEX "cbrt" → 1.86579517236 + 1.07721734502 i
-16 ENTER 16 COMPLEX "cbrt" → 2 + 2 i
1e9 ENTER "cbrt" → 1000

外部サイト

"Cube root" — Wikipedia
"Padé approximant" — Wikipedia
"Puiseux series" — Wikipedia