Go/Hypot

頑強なHypot

package main

import (
	"fmt"
	"math"
)

// SimpleHypot は単純な実装
func simpleHypot(x, y float64) float64 {
	return math.Sqrt(x*x + y*y)
}

// 改善されたhypot実装
func hypot(x, y float64) float64 {
	getHighWord := func(f float64) int32 {
		return int32(math.Float64bits(f) >> 32)
	}

	getLowWord := func(f float64) int32 {
		return int32(math.Float64bits(f))
	}

	setHighWord := func(f *float64, hw int32) {
		bits := math.Float64bits(*f)
		*f = math.Float64frombits(uint64(hw)<<32 | uint64(bits&0xFFFFFFFF))
	}

	ha := getHighWord(x) & 0x7fffffff
	hb := getHighWord(y) & 0x7fffffff

	a, b := x, y
	if hb > ha {
		a, b = y, x
		ha, hb = hb, ha
	}

	a = math.Abs(a)
	b = math.Abs(b)

	if ha-hb > 0x3c00000 {
		return a + b
	}

	k := 0
	if ha > 0x5f300000 {
		if ha >= 0x7ff00000 {
			w := math.Abs(x+0.0) - math.Abs(y+0.0)
			low := getLowWord(a)
			if (ha&0xfffff | low) == 0 {
				w = a
			}
			low = getLowWord(b)
			if (hb ^ 0x7ff00000 | low) == 0 {
				w = b
			}
			return w
		}
		ha -= 0x25800000
		hb -= 0x25800000
		k += 600
		setHighWord(&a, ha)
		setHighWord(&b, hb)
	}

	if hb < 0x20b00000 {
		if hb <= 0x000fffff {
			low := getLowWord(b)
			if (hb | low) == 0 {
				return a
			}
			t1 := 0.0
			setHighWord(&t1, 0x7fd00000)
			b *= t1
			a *= t1
			k -= 1022
		} else {
			ha += 0x25800000
			hb += 0x25800000
			k -= 600
			setHighWord(&a, ha)
			setHighWord(&b, hb)
		}
	}

	w := a - b
	if w > b {
		t1 := 0.0
		setHighWord(&t1, ha)
		t2 := a - t1
		w = math.Sqrt(t1*t1 - (b*(-b) - t2*(a+t1)))
	} else {
		a = a + a
		y1 := 0.0
		setHighWord(&y1, hb)
		y2 := b - y1
		t1 := 0.0
		setHighWord(&t1, ha+0x00100000)
		t2 := a - t1
		w = math.Sqrt(t1*y1 - (w*(-w) - (t1*y2 + t2*b)))
	}

	if k != 0 {
		t1 := 0.0
		setHighWord(&t1, int32((1023+k)<<20))
		return t1 * w
	}
	return w
}

// 比較関数
func compareResults(x, y float64) {
	h1 := hypot(x, y)
	h2 := simpleHypot(x, y)
	h3 := math.Hypot(x, y)

	fmt.Printf("\nTest case: x=%g, y=%g\n", x, y)
	fmt.Printf("Custom hypot:     %g\n", h1)
	fmt.Printf("Simple hypot:     %g\n", h2)
	fmt.Printf("Standard hypot:   %g\n", h3)

	// NaNチェックを含む比較
	if math.IsNaN(h1) || math.IsNaN(h2) || math.IsNaN(h3) {
		fmt.Printf("NaN検出！\n")
		fmt.Printf("Custom IsNaN:   %v\n", math.IsNaN(h1))
		fmt.Printf("Simple IsNaN:   %v\n", math.IsNaN(h2))
		fmt.Printf("Standard IsNaN: %v\n", math.IsNaN(h3))
		return
	}

	// 結果が異なる場合の相対誤差計算
	if h1 != h2 || h1 != h3 {
		fmt.Printf("差異検出！\n")
		if !math.IsInf(h1, 0) && !math.IsInf(h2, 0) {
			relErr12 := math.Abs(h1-h2) / math.Max(math.Abs(h1), math.Abs(h2))
			fmt.Printf("Custom vs Simple 相対誤差: %g\n", relErr12)
		} else {
			fmt.Printf("Custom vs Simple: 少なくとも1つが無限大\n")
		}

		if !math.IsInf(h1, 0) && !math.IsInf(h3, 0) {
			relErr13 := math.Abs(h1-h3) / math.Max(math.Abs(h1), math.Abs(h3))
			fmt.Printf("Custom vs Standard 相対誤差: %g\n", relErr13)
		} else {
			fmt.Printf("Custom vs Standard: 少なくとも1つが無限大\n")
		}
	}
}

func main() {
	fmt.Println("=== 通常のケース ===")
	compareResults(3, 4)
	compareResults(1e-8, 1e-8)

	fmt.Println("\n=== 極端な値のケース ===")
	compareResults(1e300, 1e300)
	compareResults(1e-300, 1e-300)

	fmt.Println("\n=== 非常に異なる大きさの値 ===")
	compareResults(1e308, 1e-308)
	compareResults(1e50, 1e-50)

	fmt.Println("\n=== 特殊な値 ===")
	compareResults(math.Inf(1), 1.0)
	compareResults(math.NaN(), 1.0)
	compareResults(0.0, 0.0)

	fmt.Println("\n=== オーバーフローが起こりやすいケース ===")
	compareResults(1e308, 1e308)
	compareResults(math.MaxFloat64/2, math.MaxFloat64/2)

	fmt.Println("\n=== アンダーフローが起こりやすいケース ===")
	compareResults(math.SmallestNonzeroFloat64, math.SmallestNonzeroFloat64)
	compareResults(1e-308, 1e-308)
}

実行結果

=== 通常のケース ===

Test case: x=3, y=4
Custom hypot:     5
Simple hypot:     5
Standard hypot:   5

Test case: x=1e-08, y=1e-08
Custom hypot:     1.4142135623730952e-08
Simple hypot:     1.4142135623730952e-08
Standard hypot:   1.4142135623730952e-08

=== 極端な値のケース ===

Test case: x=1e+300, y=1e+300
Custom hypot:     1.4142135623730952e+300
Simple hypot:     +Inf
Standard hypot:   1.4142135623730952e+300
差異検出！
Custom vs Simple: 少なくとも1つが無限大
Custom vs Standard 相対誤差: 0

Test case: x=1e-300, y=1e-300
Custom hypot:     1.414213562373095e-300
Simple hypot:     0
Standard hypot:   1.4142135623730952e-300
差異検出！
Custom vs Simple 相対誤差: 1
Custom vs Standard 相対誤差: 1.1722481355006683e-16

=== 非常に異なる大きさの値 ===

Test case: x=1e+308, y=1e-308
Custom hypot:     1e+308
Simple hypot:     +Inf
Standard hypot:   1e+308
差異検出！
Custom vs Simple: 少なくとも1つが無限大
Custom vs Standard 相対誤差: 0

Test case: x=1e+50, y=1e-50
Custom hypot:     1e+50
Simple hypot:     1e+50
Standard hypot:   1e+50

=== 特殊な値 ===

Test case: x=+Inf, y=1
Custom hypot:     +Inf
Simple hypot:     +Inf
Standard hypot:   +Inf

Test case: x=NaN, y=1
Custom hypot:     NaN
Simple hypot:     NaN
Standard hypot:   NaN
NaN検出！
Custom IsNaN:   true
Simple IsNaN:   true
Standard IsNaN: true

Test case: x=0, y=0
Custom hypot:     0
Simple hypot:     0
Standard hypot:   0

=== オーバーフローが起こりやすいケース ===

Test case: x=1e+308, y=1e+308
Custom hypot:     1.4142135623730951e+308
Simple hypot:     +Inf
Standard hypot:   1.4142135623730951e+308
差異検出！
Custom vs Simple: 少なくとも1つが無限大
Custom vs Standard 相対誤差: 0

Test case: x=8.988465674311579e+307, y=8.988465674311579e+307
Custom hypot:     1.2711610061536462e+308
Simple hypot:     +Inf
Standard hypot:   1.2711610061536462e+308
差異検出！
Custom vs Simple: 少なくとも1つが無限大
Custom vs Standard 相対誤差: 0

=== アンダーフローが起こりやすいケース ===

Test case: x=5e-324, y=5e-324
Custom hypot:     5e-324
Simple hypot:     0
Standard hypot:   5e-324
差異検出！
Custom vs Simple 相対誤差: 1
Custom vs Standard 相対誤差: 0

Test case: x=1e-308, y=1e-308
Custom hypot:     1.414213562373095e-308
Simple hypot:     0
Standard hypot:   1.414213562373095e-308
差異検出！
Custom vs Simple 相対誤差: 1
Custom vs Standard 相対誤差: 0