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)
}