ちょっとむかしを思い出しつつ、AIに書いてもらった。1
典型的な√aを計算する。
二分法は、二分検索と同じように、lowとhighの中央を詰めていく。
ニュートン法は、$y=x^2-a$を微分して、$y'=2x$を使って、$x-\frac{y}{y'}=\frac{1}{2}(\frac{a}{x}+x)$を次の値にする。
初期値が、a>1のとき二分法は[1, a]、ニュートン法はaからスタートすればよいが、a<1のときは累乗すると小さくなるので二分法は[a, 1]、ニュートン法は1からスタートになる。
public class SquareRootUniversal {
private static final double EPSILON = 1e-10;
public static void main(String[] args) {
// a < 1 のケース (0.01) と a > 1 のケース (1e20) 両方でテスト
double[] testValues = {0.0001, 0.5, 2.0, 1e20};
for (double a : testValues) {
System.out.println("=== Target: sqrt(" + a + ") ===");
solveNewton(a);
solveBisection(a);
System.out.println();
}
}
/**
* 汎用ニュートン法
*/
public static void solveNewton(double a) {
if (a < 0) return;
if (a == 0) { System.out.println("Newton: 0.0"); return; }
// ポイント: a < 1 のときは 1.0 を初期値にすることで、必ず「解より大きい側」からスタートできる
double x = Math.max(1.0, a);
int steps = 0;
while (true) {
steps++;
double nextX = 0.5 * (x + a / x);
System.out.printf("Step %d: x = %.10e, nextX = %.10e\n", steps, x, nextX);
if (Math.abs(nextX - x) < EPSILON) {
x = nextX;
break;
}
x = nextX;
}
System.out.printf("Newton Result: %.15f (Steps: %d)\n", x, steps);
}
/**
* 汎用二分法
*/
public static void solveBisection(double a) {
if (a < 0) return;
if (a == 0) { System.out.println("Bisection: 0.0"); return; }
// a < 1 のときは範囲を [a, 1]、a >= 1 のときは [1, a] に設定する
double low = Math.min(a, 1.0);
double high = Math.max(a, 1.0);
int steps = 0;
while ((high - low) / high > EPSILON) {
steps++;
double mid = (low + high) / 2.0;
System.out.printf("Step %d: mid = %.10e, low = %.10e, high = %.10e\n", steps, mid, low, high);
if (mid * mid > a) {
high = mid;
} else {
low = mid;
}
}
System.out.printf("Bisection Result: %.15f (Steps: %d)\n", (low + high) / 2.0, steps);
}
}
ニュートン法の方が収束が早い。
ただし、1e20のとき、半分づつしか詰めていけない。
=== Target: sqrt(1.0E-4) ===
Step 1: x = 1.0000000000e+00, nextX = 5.0005000000e-01
Step 2: x = 5.0005000000e-01, nextX = 2.5012499000e-01
Step 3: x = 2.5012499000e-01, nextX = 1.2526239506e-01
Step 4: x = 1.2526239506e-01, nextX = 6.3030359624e-02
Step 5: x = 6.3030359624e-02, nextX = 3.2308448330e-02
Step 6: x = 3.2308448330e-02, nextX = 1.7701806998e-02
Step 7: x = 1.7701806998e-02, nextX = 1.1675473895e-02
Step 8: x = 1.1675473895e-02, nextX = 1.0120218365e-02
Step 9: x = 1.0120218365e-02, nextX = 1.0000714039e-02
Step 10: x = 1.0000714039e-02, nextX = 1.0000000025e-02
Step 11: x = 1.0000000025e-02, nextX = 1.0000000000e-02
Newton Result: 0.010000000000000 (Steps: 11)
Step 1: mid = 5.0005000000e-01, low = 1.0000000000e-04, high = 1.0000000000e+00
Step 2: mid = 2.5007500000e-01, low = 1.0000000000e-04, high = 5.0005000000e-01
Step 3: mid = 1.2508750000e-01, low = 1.0000000000e-04, high = 2.5007500000e-01
Step 4: mid = 6.2593750000e-02, low = 1.0000000000e-04, high = 1.2508750000e-01
Step 5: mid = 3.1346875000e-02, low = 1.0000000000e-04, high = 6.2593750000e-02
：
Step 36: mid = 1.0000000003e-02, low = 9.9999999888e-03, high = 1.0000000018e-02
Step 37: mid = 9.9999999960e-03, low = 9.9999999888e-03, high = 1.0000000003e-02
Step 38: mid = 9.9999999997e-03, low = 9.9999999960e-03, high = 1.0000000003e-02
Step 39: mid = 1.0000000001e-02, low = 9.9999999997e-03, high = 1.0000000003e-02
Step 40: mid = 1.0000000001e-02, low = 9.9999999997e-03, high = 1.0000000001e-02
Bisection Result: 0.010000000000131 (Steps: 40)
=== Target: sqrt(0.5) ===
Step 1: x = 1.0000000000e+00, nextX = 7.5000000000e-01
Step 2: x = 7.5000000000e-01, nextX = 7.0833333333e-01
Step 3: x = 7.0833333333e-01, nextX = 7.0710784314e-01
Step 4: x = 7.0710784314e-01, nextX = 7.0710678119e-01
Step 5: x = 7.0710678119e-01, nextX = 7.0710678119e-01
Newton Result: 0.707106781186548 (Steps: 5)
Bisection Result: 0.707106781163020 (Steps: 33)
=== Target: sqrt(2.0) ===
Step 1: x = 2.0000000000e+00, nextX = 1.5000000000e+00
Step 2: x = 1.5000000000e+00, nextX = 1.4166666667e+00
Step 3: x = 1.4166666667e+00, nextX = 1.4142156863e+00
Step 4: x = 1.4142156863e+00, nextX = 1.4142135624e+00
Step 5: x = 1.4142135624e+00, nextX = 1.4142135624e+00
Newton Result: 1.414213562373095 (Steps: 5)
Bisection Result: 1.414213562326040 (Steps: 33)
=== Target: sqrt(1.0E20) ===
Step 1: x = 1.0000000000e+20, nextX = 5.0000000000e+19
Step 2: x = 5.0000000000e+19, nextX = 2.5000000000e+19
Step 3: x = 2.5000000000e+19, nextX = 1.2500000000e+19
Step 4: x = 1.2500000000e+19, nextX = 6.2500000000e+18
Step 5: x = 6.2500000000e+18, nextX = 3.1250000000e+18
Step 6: x = 3.1250000000e+18, nextX = 1.5625000000e+18
Step 7: x = 1.5625000000e+18, nextX = 7.8125000000e+17
Step 8: x = 7.8125000000e+17, nextX = 3.9062500000e+17
Step 9: x = 3.9062500000e+17, nextX = 1.9531250000e+17
Step 10: x = 1.9531250000e+17, nextX = 9.7656250000e+16
Step 11: x = 9.7656250000e+16, nextX = 4.8828125000e+16
Step 12: x = 4.8828125000e+16, nextX = 2.4414062500e+16
Step 13: x = 2.4414062500e+16, nextX = 1.2207031250e+16
Step 14: x = 1.2207031250e+16, nextX = 6.1035156250e+15
Step 15: x = 6.1035156250e+15, nextX = 3.0517578125e+15
Step 16: x = 3.0517578125e+15, nextX = 1.5258789063e+15
Step 17: x = 1.5258789063e+15, nextX = 7.6293945317e+14
Step 18: x = 7.6293945317e+14, nextX = 3.8146972665e+14
Step 19: x = 3.8146972665e+14, nextX = 1.9073486346e+14
Step 20: x = 1.9073486346e+14, nextX = 9.5367431990e+13
Step 21: x = 9.5367431990e+13, nextX = 4.7683716519e+13
Step 22: x = 4.7683716519e+13, nextX = 2.3841859308e+13
Step 23: x = 2.3841859308e+13, nextX = 1.1920931751e+13
Step 24: x = 1.1920931751e+13, nextX = 5.9604700699e+12
Step 25: x = 5.9604700699e+12, nextX = 2.9802434236e+12
Step 26: x = 2.9802434236e+12, nextX = 1.4901384889e+12
Step 27: x = 1.4901384889e+12, nextX = 7.4510279840e+11
Step 28: x = 7.4510279840e+11, nextX = 3.7261850403e+11
Step 29: x = 3.7261850403e+11, nextX = 1.8644343752e+11
Step 30: x = 1.8644343752e+11, nextX = 9.3489896607e+10
Step 31: x = 9.3489896607e+10, nextX = 4.7279765453e+10
Step 32: x = 4.7279765453e+10, nextX = 2.4697417583e+10
Step 33: x = 2.4697417583e+10, nextX = 1.4373211954e+10
Step 34: x = 1.4373211954e+10, nextX = 1.0665299547e+10
Step 35: x = 1.0665299547e+10, nextX = 1.0020750636e+10
Step 36: x = 1.0020750636e+10, nextX = 1.0000021485e+10
Step 37: x = 1.0000021485e+10, nextX = 1.0000000000e+10
Step 38: x = 1.0000000000e+10, nextX = 1.0000000000e+10
Step 39: x = 1.0000000000e+10, nextX = 1.0000000000e+10
Newton Result: 10000000000.000000000000000 (Steps: 39)
Bisection Result: 10000000000.005482000000000 (Steps: 67)
指数部を半分にする
そこで出してきたのが、double型の指数部を半分にした初期値にしよう。
ソース見ても分からないので、デバッグを入れまくった。
if (a <= 0) return;
// --- 1. ビット演算による初期値の推定 ---
// IEEE 754 形式のビットパターンを取得
long bits = Double.doubleToLongBits(a);
System.out.printf("%d %d %.10e\n", (bits >> 52)-1023, (bits & 0x000FFFFFFFFFFFFFL), a);
// 指数部を半分にする魔法の操作
// 指数部はバイアス(1023)がかかっているので、それを考慮してビットをずらす
long initialBits;
double x;
initialBits = (bits >> 1) + (511L << 52);
x = Double.longBitsToDouble(initialBits);
System.out.printf("%d %d %.10e\n", (initialBits >> 52)-1023, (initialBits & 0x000FFFFFFFFFFFFFL), x);
initialBits = (bits >> 1) + (512L << 52);
x = Double.longBitsToDouble(initialBits);
System.out.printf("%d %d %.10e\n", (initialBits >> 52)-1023, (initialBits & 0x000FFFFFFFFFFFFFL), x);
initialBits = (((bits >> 52) - 1023) / 2 + 1023) << 52 | (bits & 0x000FFFFFFFFFFFFFL);
x = Double.longBitsToDouble(initialBits);
System.out.printf("%d %d %.10e\n", (initialBits >> 52)-1023, (initialBits & 0x000FFFFFFFFFFFFFL), x);
System.out.printf("Estimated Initial x0: %.10e\n", x);
初期値を1.16e10からスタートすると、39ステップが5ステップで収束する。
Target: sqrt(1.0E20) = 1.0E10
66 1599915997629504 1.0000000000e+20
32 3051757812500000 7.2053503417e+09
33 3051757812500000 1.4410700683e+10
33 1599915997629504 1.1641532183e+10
Estimated Initial x0: 1.1641532183e+10
Step 1: x = 1.1641532183e+10, nextX = 1.0115733387e+10
Step 2: x = 1.0115733387e+10, nextX = 1.0000662049e+10
Step 3: x = 1.0000662049e+10, nextX = 1.0000000022e+10
Step 4: x = 1.0000000022e+10, nextX = 1.0000000000e+10
Step 5: x = 1.0000000000e+10, nextX = 1.0000000000e+10
Fast Newton Result: 10000000000.000000000000000 (Steps: 5)
66 1599915997629504 1.0000000000e+20
求めたい数aの、バイアスを戻した指数（$2^{exponent}$）、仮数（$1+fraction×2^{-52}$）、double値（$2^{66}×(1+\frac{1599915997629504}{2^{52}})=10^{20}$）
initialBits = (bits >> 1) + (511L << 52);
32 3051757812500000 7.2053503417e+09
AIが出してきたやつ。まるごと1ビット右シフトして、指数部のバイアスが1023なので、半分の511を（指数部の52ビット目に合わせて）足してやれ。
仮数部も変わっている。指数部が切り捨てた511なので、平方根の値よりも小さい数になっている。
initialBits = (bits >> 1) + (512L << 52);
33 3051757812500000 1.4410700683e+10
そこで、511でなく、512の方が、大きくなるのではと、指数部がさっきより1大きい。2
initialBits = (((bits >> 52) - 1023) / 2 + 1023) << 52 | (bits & 0x000FFFFFFFFFFFFFL);
33 1599915997629504 1.1641532183e+10
仮数部を変えずに、指数部のみ半分にしてみた。
2乗根以外のk乗根のときにkで割ればよい。
なお、符号ビットもシフトされるのでは？と思ったら、冒頭で弾いている。
ちなみに三乗根の初期値はこれだと出してきた。
もう知らんがな。
public static double initialCubeRoot(double a) {
long bits = Double.doubleToLongBits(a);
// 指数部をだいたい1/3にするマジックナンバー
// 0x4000000000000000L 付近の調整値が必要
long initialBits = bits / 3 + 0x2A51000000000000L;
return Double.longBitsToDouble(initialBits);
}
ちなみにk乗根は、大きく離れているとざっくり$\frac{k-1}{k}$になるので、4乗根ならば、2乗根を2回の方が早いとな。
まあ100乗根なんて計算する日は来ないだろう。
セカント法
直近の2つの近似値 $(x_{n-1}, x_n)$ を結ぶ直線を引き、その直線が $x$ 軸と交わる点を次の近似値 $x_{n+1}$ にする。
public static void solveSecant(double x0, double x1, java.util.function.DoubleUnaryOperator f) {
double xPrev = x0;
double xCurr = x1;
for (int i = 0; i < 100; i++) {
double fPrev = f.applyAsDouble(xPrev);
double fCurr = f.applyAsDouble(xCurr);
// 分母が0（または極小）にならないかチェック
if (Math.abs(fCurr - fPrev) < 1e-18) break;
// セカント法の更新式
double xNext = xCurr - fCurr * (xCurr - xPrev) / (fCurr - fPrev);
if (Math.abs(xNext - xCurr) < Math.abs(xCurr) * 1e-15) {
System.out.println("Converged at: " + xNext);
return;
}
xPrev = xCurr;
xCurr = xNext;
}
}
逆2次補間法
直近の3点を通る放物線（2次式）を使って解を探す。
$x$ と $y$ の役割を入れ替えて、$x = g(y)$ という「横向きの放物線」として近似するのが「逆」2次補間法の特徴です。
- 通常の2次補間: $y = ax^2 + bx + c$
- 逆2次補間: $x = ay^2 + by + c$
3つの点 $(x_{n-2}, y_{n-2}), (x_{n-1}, y_{n-1}), (x_n, y_n)$ があるとき、次の近似値 $x_{n+1}$ は以下の式で計算される。
$$x_{n+1} = \frac{y_{n-1}y_n}{(y_{n-2}-y_{n-1})(y_{n-2}-y_n)}x_{n-2} + \frac{y_{n-2}y_n}{(y_{n-1}-y_{n-2})(y_{n-1}-y_n)}x_{n-1} + \frac{y_{n-2}y_{n-1}}{(y_n-y_{n-2})(y_n-y_{n-1})}x_n$$
ブレント法（Brent's method）
3つのハイブリッド
- 逆2次補間法 (Inverse Quadratic Interpolation)
- セカント法 (Secant Method)
- 二分法 (Bisection Method)
Commons Mathにあるのか。
// Apache Commons Math を使った例
UnivariateSolver solver = new BrentSolver();
double result = solver.solve(100, x -> Math.cos(x) - x, 0, 1);
ハレー法
ハレー彗星でお馴染みのおっさんが提案した方法。って、1700年ころ？
ニュートン法の式に、$f(x)$ の2階微分（$f''(x)$）の情報を加えることで、グラフの「曲がり具合」を計算に取り入れる。
$$x_{n+1} = x_n - \frac{2 f(x_n) f'(x_n)}{2 (f'(x_n))^2 - f(x_n) f''(x_n)}$$
高速逆平方根（Fast Inverse Square Root）
$1/\sqrt{x}$ を計算するための関数。
threehalfsって、3次元ベクトルの長さの逆数を計算したいためか。
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // floatをlongとして解釈（ビットの強引な取り出し）
i = 0x5f3759df - ( i >> 1 ); // ★謎のマジックナンバー★
y = * ( float * ) &i; // 再びfloatに戻す
y = y * ( threehalfs - ( x2 * y * y ) ); // ニュートン法（1回だけ）
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2回目は不要なのでコメントアウトされている
return y;
}
現在のCPU（x86やARM）には、この「逆平方根」をハードウェアレベルで高速に行う専用命令（RSQRTSSなど）が搭載されています。そのため、現代のプログラミングでこのトリッキーなコードを自前で書く必要はなくなりました。
って、ライブラリに関数がないと使えないわな。