利用者:Ef3/Hypot

hypot.go

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=%.22g, y=%.22g\n", x, y)
	fmt.Printf("Custom hypot:     %.22g\n", h1)
	fmt.Printf("Simple hypot:     %.22g\n", h2)
	fmt.Printf("Standard hypot:   %.22g\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)
}

hypot.js

// IEEE 754 Float64用のhypot実装
function hypot(x, y) {
    // Float64のビット操作ヘルパー
    function getHighWord(f) {
        // Float64ArrayとUint32Arrayで上位32ビットを取得
        const float64 = new Float64Array(1);
        float64[0] = f;
        const uint32 = new Uint32Array(float64.buffer);
        // JavaScriptはリトルエンディアンが一般的なので、インデックス1を使用
        return uint32[1];
    }

    function getLowWord(f) {
        const float64 = new Float64Array(1);
        float64[0] = f;
        const uint32 = new Uint32Array(float64.buffer);
        return uint32[0];
    }

    function setHighWord(f, hw) {
        const float64 = new Float64Array(1);
        float64[0] = f;
        const uint32 = new Uint32Array(float64.buffer);
        uint32[1] = hw;
        return float64[0];
    }

    let ha = getHighWord(x) & 0x7fffffff;
    let hb = getHighWord(y) & 0x7fffffff;

    let a, b;
    if (hb > ha) {
        [a, b] = [y, x];
        [ha, hb] = [hb, ha];
    } else {
        [a, b] = [x, y];
    }

    a = Math.abs(a);
    b = Math.abs(b);

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

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

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

    let w = a - b;
    if (w > b) {
        const t1 = setHighWord(0, ha);
        const t2 = a - t1;
        w = Math.sqrt(t1 * t1 - (b * (-b) - t2 * (a + t1)));
    } else {
        a = a + a;
        const y1 = setHighWord(0, hb);
        const y2 = b - y1;
        const t1 = setHighWord(0, ha + 0x00100000);
        const t2 = a - t1;
        w = Math.sqrt(t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)));
    }

    if (k !== 0) {
        const t1 = setHighWord(0, (1023 + k) << 20);
        return t1 * w;
    }
    return w;
}

// テストコード
function compareResults(x, y) {
    const h1 = hypot(x, y);
    const h2 = Math.sqrt(x * x + y * y);
    const h3 = Math.hypot(x, y);

    console.log(`\nTest case: x=${x.toPrecision(17)}, y=${y.toPrecision(17)}`);
    console.log(`Custom hypot:     ${h1.toPrecision(17)}`);
    console.log(`Simple hypot:     ${h2.toPrecision(17)}`);
    console.log(`Standard hypot:   ${h3.toPrecision(17)}`);

    // NaNチェック
    if (Number.isNaN(h1) || Number.isNaN(h2) || Number.isNaN(h3)) {
        console.log("NaN検出！");
        console.log(`Custom IsNaN:   ${Number.isNaN(h1)}`);
        console.log(`Simple IsNaN:   ${Number.isNaN(h2)}`);
        console.log(`Standard IsNaN: ${Number.isNaN(h3)}`);
        return;
    }

    // 結果比較
    if (h1 !== h2 || h1 !== h3) {
        console.log("差異検出！");
        if (!Number.isFinite(h1) || !Number.isFinite(h2)) {
            console.log("Custom vs Simple: 少なくとも1つが無限大");
        } else {
            const relErr12 = Math.abs(h1 - h2) / Math.max(Math.abs(h1), Math.abs(h2));
            console.log(`Custom vs Simple 相対誤差: ${relErr12.toPrecision(17)}`);
        }

        if (!Number.isFinite(h1) || !Number.isFinite(h3)) {
            console.log("Custom vs Standard: 少なくとも1つが無限大");
        } else {
            const relErr13 = Math.abs(h1 - h3) / Math.max(Math.abs(h1), Math.abs(h3));
            console.log(`Custom vs Standard 相対誤差: ${relErr13.toPrecision(17)}`);
        }
    }
}

// テストの実行
console.log("=== 通常のケース ===");
compareResults(3, 4);
compareResults(1e-8, 1e-8);

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

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

console.log("\n=== 特殊な値 ===");
compareResults(Infinity, 1.0);
compareResults(NaN, 1.0);
compareResults(0.0, 0.0);

console.log("\n=== オーバーフローが起こりやすいケース ===");
compareResults(1e308, 1e308);
compareResults(Number.MAX_VALUE/2, Number.MAX_VALUE/2);

console.log("\n=== アンダーフローが起こりやすいケース ===");
compareResults(Number.MIN_VALUE, Number.MIN_VALUE);
compareResults(1e-308, 1e-308);

hypot.rs

use std::f64::{self, INFINITY, NAN};

// 単純な実装
fn simple_hypot(x: f64, y: f64) -> f64 {
    (x * x + y * y).sqrt()
}

// 改善された実装
fn hypot(mut x: f64, mut y: f64) -> f64 {
    // ビット操作のためのヘルパー関数
    fn get_high_word(f: f64) -> i32 {
        (f.to_bits() >> 32) as i32
    }

    fn get_low_word(f: f64) -> i32 {
        f.to_bits() as i32
    }

    fn set_high_word(f: &mut f64, hw: i32) {
        let bits = f.to_bits();
        *f = f64::from_bits((hw as u64) << 32 | (bits & 0xFFFFFFFF));
    }

    let mut ha = get_high_word(x) & 0x7fffffff;
    let mut hb = get_high_word(y) & 0x7fffffff;
    
    let mut a = x;
    let mut b = y;

    if hb > ha {
        std::mem::swap(&mut a, &mut b);
        std::mem::swap(&mut ha, &mut hb);
    }

    a = a.abs();
    b = b.abs();

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

    let mut k = 0;
    if ha > 0x5f300000 {
        if ha >= 0x7ff00000 {
            let w = (x + 0.0).abs() - (y + 0.0).abs();
            let low = get_low_word(a);
            if (ha & 0xfffff | low) == 0 {
                return a;
            }
            let low = get_low_word(b);
            if (hb ^ 0x7ff00000 | low) == 0 {
                return b;
            }
            return w;
        }
        ha -= 0x25800000;
        hb -= 0x25800000;
        k += 600;
        set_high_word(&mut a, ha);
        set_high_word(&mut b, hb);
    }

    if hb < 0x20b00000 {
        if hb <= 0x000fffff {
            let low = get_low_word(b);
            if (hb | low) == 0 {
                return a;
            }
            let mut t1 = 0.0;
            set_high_word(&mut t1, 0x7fd00000);
            b *= t1;
            a *= t1;
            k -= 1022;
        } else {
            ha += 0x25800000;
            hb += 0x25800000;
            k -= 600;
            set_high_word(&mut a, ha);
            set_high_word(&mut b, hb);
        }
    }

    let w = a - b;
    if w > b {
        let mut t1 = 0.0;
        set_high_word(&mut t1, ha);
        let t2 = a - t1;
        let w = ((t1 * t1) - (b * (-b) - t2 * (a + t1))).sqrt();
        if k != 0 {
            let mut t1 = 0.0;
            set_high_word(&mut t1, ((1023 + k) << 20) as i32);
            return t1 * w;
        }
        return w;
    } else {
        a = a + a;
        let mut y1 = 0.0;
        set_high_word(&mut y1, hb);
        let y2 = b - y1;
        let mut t1 = 0.0;
        set_high_word(&mut t1, ha + 0x00100000);
        let t2 = a - t1;
        let w = ((t1 * y1) - (w * (-w) - (t1 * y2 + t2 * b))).sqrt();
        if k != 0 {
            let mut t1 = 0.0;
            set_high_word(&mut t1, ((1023 + k) << 20) as i32);
            return t1 * w;
        }
        return w;
    }
}

fn compare_results(x: f64, y: f64) {
    let h1 = hypot(x, y);
    let h2 = simple_hypot(x, y);
    let h3 = f64::hypot(x, y);

    println!("\nTest case: x={:e}, y={:e}", x, y);
    println!("Custom hypot:     {:e}", h1);
    println!("Simple hypot:     {:e}", h2);
    println!("Standard hypot:   {:e}", h3);

    // NaNチェック
    if h1.is_nan() || h2.is_nan() || h3.is_nan() {
        println!("NaN検出！");
        println!("Custom IsNaN:   {}", h1.is_nan());
        println!("Simple IsNaN:   {}", h2.is_nan());
        println!("Standard IsNaN: {}", h3.is_nan());
        return;
    }

    // 結果が異なる場合の相対誤差計算
    if h1 != h2 || h1 != h3 {
        println!("差異検出！");
        if !h1.is_infinite() && !h2.is_infinite() {
            let rel_err12 = (h1 - h2).abs() / h1.abs().max(h2.abs());
            println!("Custom vs Simple 相対誤差: {:e}", rel_err12);
        } else {
            println!("Custom vs Simple: 少なくとも1つが無限大");
        }
        if !h1.is_infinite() && !h3.is_infinite() {
            let rel_err13 = (h1 - h3).abs() / h1.abs().max(h3.abs());
            println!("Custom vs Standard 相対誤差: {:e}", rel_err13);
        } else {
            println!("Custom vs Standard: 少なくとも1つが無限大");
        }
    }
}

fn main() {
    println!("=== 通常のケース ===");
    compare_results(3.0, 4.0);
    compare_results(1e-8, 1e-8);

    println!("\n=== 極端な値のケース ===");
    compare_results(1e300, 1e300);
    compare_results(1e-300, 1e-300);

    println!("\n=== 非常に異なる大きさの値 ===");
    compare_results(1e308, 1e-308);
    compare_results(1e50, 1e-50);

    println!("\n=== 特殊な値 ===");
    compare_results(INFINITY, 1.0);
    compare_results(NAN, 1.0);
    compare_results(0.0, 0.0);

    println!("\n=== オーバーフローが起こりやすいケース ===");
    compare_results(1e308, 1e308);
    compare_results(f64::MAX/2.0, f64::MAX/2.0);

    println!("\n=== アンダーフローが起こりやすいケース ===");
    compare_results(f64::MIN_POSITIVE, f64::MIN_POSITIVE);
    compare_results(1e-308, 1e-308);
}

hypot.c++

#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstdint>
#include <iomanip>

// ビット操作のためのユーティリティ関数
union DoubleBits {
    double d;
    uint64_t bits;
};

int32_t getHighWord(double f) {
    DoubleBits db{f};
    return static_cast<int32_t>(db.bits >> 32);
}

int32_t getLowWord(double f) {
    DoubleBits db{f};
    return static_cast<int32_t>(db.bits & 0xFFFFFFFF);
}

void setHighWord(double& f, int32_t hw) {
    DoubleBits db{f};
    db.bits = (static_cast<uint64_t>(hw) << 32) | (db.bits & 0xFFFFFFFF);
    f = db.d;
}

// 単純な実装
double simpleHypot(double x, double y) {
    return std::sqrt(x * x + y * y);
}

// 改善された実装
double hypot(double x, double y) {
    int32_t ha = getHighWord(x) & 0x7fffffff;
    int32_t hb = getHighWord(y) & 0x7fffffff;
    
    double a = x;
    double b = y;

    if (hb > ha) {
        std::swap(a, b);
        std::swap(ha, hb);
    }

    a = std::abs(a);
    b = std::abs(b);

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

    int k = 0;
    if (ha > 0x5f300000) {
        if (ha >= 0x7ff00000) {
            double w = std::abs(x + 0.0) - std::abs(y + 0.0);
            int32_t low = getLowWord(a);
            if ((ha & 0xfffff | low) == 0) {
                return a;
            }
            low = getLowWord(b);
            if ((hb ^ 0x7ff00000 | low) == 0) {
                return b;
            }
            return w;
        }
        ha -= 0x25800000;
        hb -= 0x25800000;
        k += 600;
        setHighWord(a, ha);
        setHighWord(b, hb);
    }

    if (hb < 0x20b00000) {
        if (hb <= 0x000fffff) {
            int32_t low = getLowWord(b);
            if ((hb | low) == 0) {
                return a;
            }
            double 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);
        }
    }

    double w = a - b;
    if (w > b) {
        double t1 = 0.0;
        setHighWord(t1, ha);
        double t2 = a - t1;
        w = std::sqrt(t1 * t1 - (b * (-b) - t2 * (a + t1)));
        if (k != 0) {
            double t1 = 0.0;
            setHighWord(t1, (1023 + k) << 20);
            return t1 * w;
        }
        return w;
    } else {
        a = a + a;
        double y1 = 0.0;
        setHighWord(y1, hb);
        double y2 = b - y1;
        double t1 = 0.0;
        setHighWord(t1, ha + 0x00100000);
        double t2 = a - t1;
        w = std::sqrt(t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)));
        if (k != 0) {
            double t1 = 0.0;
            setHighWord(t1, (1023 + k) << 20);
            return t1 * w;
        }
        return w;
    }
}

void compareResults(double x, double y) {
    double h1 = hypot(x, y);
    double h2 = simpleHypot(x, y);
    double h3 = std::hypot(x, y);

    std::cout << std::scientific << std::setprecision(20);
    std::cout << "\nTest case: x=" << x << ", y=" << y << std::endl;
    std::cout << "Custom hypot:     " << h1 << std::endl;
    std::cout << "Simple hypot:     " << h2 << std::endl;
    std::cout << "Standard hypot:   " << h3 << std::endl;

    // NaNチェック
    if (std::isnan(h1) || std::isnan(h2) || std::isnan(h3)) {
        std::cout << "NaN検出！" << std::endl;
        std::cout << "Custom IsNaN:   " << std::isnan(h1) << std::endl;
        std::cout << "Simple IsNaN:   " << std::isnan(h2) << std::endl;
        std::cout << "Standard IsNaN: " << std::isnan(h3) << std::endl;
        return;
    }

    // 結果が異なる場合の相対誤差計算
    if (h1 != h2 || h1 != h3) {
        std::cout << "差異検出！" << std::endl;
        if (!std::isinf(h1) && !std::isinf(h2)) {
            double relErr12 = std::abs(h1 - h2) / std::max(std::abs(h1), std::abs(h2));
            std::cout << "Custom vs Simple 相対誤差: " << relErr12 << std::endl;
        } else {
            std::cout << "Custom vs Simple: 少なくとも1つが無限大" << std::endl;
        }
        if (!std::isinf(h1) && !std::isinf(h3)) {
            double relErr13 = std::abs(h1 - h3) / std::max(std::abs(h1), std::abs(h3));
            std::cout << "Custom vs Standard 相対誤差: " << relErr13 << std::endl;
        } else {
            std::cout << "Custom vs Standard: 少なくとも1つが無限大" << std::endl;
        }
    }
}

int main() {
    std::cout << "=== 通常のケース ===" << std::endl;
    compareResults(3.0, 4.0);
    compareResults(1e-8, 1e-8);

    std::cout << "\n=== 極端な値のケース ===" << std::endl;
    compareResults(1e300, 1e300);
    compareResults(1e-300, 1e-300);

    std::cout << "\n=== 非常に異なる大きさの値 ===" << std::endl;
    compareResults(1e308, 1e-308);
    compareResults(1e50, 1e-50);

    std::cout << "\n=== 特殊な値 ===" << std::endl;
    compareResults(INFINITY, 1.0);
    compareResults(NAN, 1.0);
    compareResults(0.0, 0.0);

    std::cout << "\n=== オーバーフローが起こりやすいケース ===" << std::endl;
    compareResults(1e308, 1e308);
    compareResults(std::numeric_limits<double>::max()/2.0, 
                  std::numeric_limits<double>::max()/2.0);

    std::cout << "\n=== アンダーフローが起こりやすいケース ===" << std::endl;
    compareResults(std::numeric_limits<double>::min(), 
                  std::numeric_limits<double>::min());
    compareResults(1e-308, 1e-308);

    return 0;
}

hypot.py

import math
import struct

# ビット操作のためのヘルパー関数
def get_high_word(f: float) -> int:
    bits = struct.unpack('Q', struct.pack('d', f))[0]
    return (bits >> 32) & 0xFFFFFFFF

def get_low_word(f: float) -> int:
    bits = struct.unpack('Q', struct.pack('d', f))[0]
    return bits & 0xFFFFFFFF

def set_high_word(f: float, hw: int) -> float:
    bits = struct.unpack('Q', struct.pack('d', f))[0]
    bits = (bits & 0xFFFFFFFF) | ((hw & 0xFFFFFFFF) << 32)
    return struct.unpack('d', struct.pack('Q', bits))[0]

# 単純な実装
def simple_hypot(x: float, y: float) -> float:
    return math.sqrt(x * x + y * y)

# 改善された実装
def hypot(x: float, y: float) -> float:
    ha = get_high_word(x) & 0x7fffffff
    hb = get_high_word(y) & 0x7fffffff
    
    if hb > ha:
        x, y = y, x
        ha, hb = hb, ha
    
    a = abs(x)
    b = abs(y)
    
    if ha - hb > 0x3c00000:
        return a + b
    
    k = 0
    if ha > 0x5f300000:
        if ha >= 0x7ff00000:
            w = abs(x + 0.0) - abs(y + 0.0)
            low = get_low_word(a)
            if (ha & 0xfffff | low) == 0:
                return a
            low = get_low_word(b)
            if (hb ^ 0x7ff00000 | low) == 0:
                return b
            return w
        ha -= 0x25800000
        hb -= 0x25800000
        k += 600
        a = set_high_word(a, ha)
        b = set_high_word(b, hb)
    
    if hb < 0x20b00000:
        if hb <= 0x000fffff:
            low = get_low_word(b)
            if (hb | low) == 0:
                return a
            t1 = set_high_word(0.0, 0x7fd00000)
            b *= t1
            a *= t1
            k -= 1022
        else:
            ha += 0x25800000
            hb += 0x25800000
            k -= 600
            a = set_high_word(a, ha)
            b = set_high_word(b, hb)
    
    w = a - b
    if w > b:
        t1 = set_high_word(0.0, ha)
        t2 = a - t1
        w = math.sqrt(t1 * t1 - (b * (-b) - t2 * (a + t1)))
        if k != 0:
            return set_high_word(0.0, ((1023 + k) << 20)) * w
        return w
    else:
        a = a + a
        y1 = set_high_word(0.0, hb)
        y2 = b - y1
        t1 = set_high_word(0.0, ha + 0x00100000)
        t2 = a - t1
        w = math.sqrt(t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)))
        if k != 0:
            return set_high_word(0.0, ((1023 + k) << 20)) * w
        return w

def compare_results(x: float, y: float) -> None:
    h1 = hypot(x, y)
    h2 = simple_hypot(x, y)
    h3 = math.hypot(x, y)
    
    print(f"\nTest case: x={x:e}, y={y:e}")
    print(f"Custom hypot:     {h1:e}")
    print(f"Simple hypot:     {h2:e}")
    print(f"Standard hypot:   {h3:e}")
    
    # NaNチェック
    if math.isnan(h1) or math.isnan(h2) or math.isnan(h3):
        print("NaN検出！")
        print(f"Custom IsNaN:   {math.isnan(h1)}")
        print(f"Simple IsNaN:   {math.isnan(h2)}")
        print(f"Standard IsNaN: {math.isnan(h3)}")
        return
    
    # 結果が異なる場合の相対誤差計算
    if h1 != h2 or h1 != h3:
        print("差異検出！")
        if not math.isinf(h1) and not math.isinf(h2):
            rel_err12 = abs(h1 - h2) / max(abs(h1), abs(h2))
            print(f"Custom vs Simple 相対誤差: {rel_err12:e}")
        else:
            print("Custom vs Simple: 少なくとも1つが無限大")
        if not math.isinf(h1) and not math.isinf(h3):
            rel_err13 = abs(h1 - h3) / max(abs(h1), abs(h3))
            print(f"Custom vs Standard 相対誤差: {rel_err13:e}")
        else:
            print("Custom vs Standard: 少なくとも1つが無限大")

def main():
    print("=== 通常のケース ===")
    compare_results(3.0, 4.0)
    compare_results(1e-8, 1e-8)
    
    print("\n=== 極端な値のケース ===")
    compare_results(1e300, 1e300)
    compare_results(1e-300, 1e-300)
    
    print("\n=== 非常に異なる大きさの値 ===")
    compare_results(1e308, 1e-308)
    compare_results(1e50, 1e-50)
    
    print("\n=== 特殊な値 ===")
    compare_results(float('inf'), 1.0)
    compare_results(float('nan'), 1.0)
    compare_results(0.0, 0.0)
    
    print("\n=== オーバーフローが起こりやすいケース ===")
    compare_results(1e308, 1e308)
    compare_results(float('inf')/2, float('inf')/2)
    
    print("\n=== アンダーフローが起こりやすいケース ===")
    compare_results(2.2250738585072014e-308, 2.2250738585072014e-308)  # min positive
    compare_results(1e-308, 1e-308)

if __name__ == "__main__":
    main()

hypot_time.py

import math
import struct
import timeit

# ビット操作のためのヘルパー関数
def get_high_word(f: float) -> int:
    bits = struct.unpack('Q', struct.pack('d', f))[0]
    return (bits >> 32) & 0xFFFFFFFF

def set_high_word(f: float, hw: int) -> float:
    bits = struct.unpack('Q', struct.pack('d', f))[0]
    bits = (bits & 0xFFFFFFFF) | ((hw & 0xFFFFFFFF) << 32)
    return struct.unpack('d', struct.pack('Q', bits))[0]

# 単純な実装
def simple_hypot(x: float, y: float) -> float:
    return math.sqrt(x * x + y * y)

# 改善された実装
def hypot(x: float, y: float) -> float:
    ha = get_high_word(x) & 0x7fffffff
    hb = get_high_word(y) & 0x7fffffff
    
    if hb > ha:
        x, y = y, x
        ha, hb = hb, ha
    
    a = abs(x)
    b = abs(y)
    
    if ha - hb > 0x3c00000:
        return a + b
    
    k = 0
    if ha > 0x5f300000:
        if ha >= 0x7ff00000:
            w = abs(x + 0.0) - abs(y + 0.0)
            if (ha & 0xfffff | get_high_word(a) & 0xffffffff) == 0:
                return a
            if (hb ^ 0x7ff00000 | get_high_word(b) & 0xffffffff) == 0:
                return b
            return w
        ha -= 0x25800000
        hb -= 0x25800000
        k += 600
        a = set_high_word(a, ha)
        b = set_high_word(b, hb)
    
    if hb < 0x20b00000:
        if hb <= 0x000fffff:
            if hb == 0 and get_high_word(b) == 0:
                return a
            t1 = set_high_word(0.0, 0x7fd00000)
            b *= t1
            a *= t1
            k -= 1022
        else:
            ha += 0x25800000
            hb += 0x25800000
            k -= 600
            a = set_high_word(a, ha)
            b = set_high_word(b, hb)
    
    w = a - b
    if w > b:
        t1 = set_high_word(0.0, ha)
        t2 = a - t1
        w = math.sqrt(t1 * t1 - (b * (-b) - t2 * (a + t1)))
        if k != 0:
            return set_high_word(0.0, ((1023 + k) << 20)) * w
        return w
    else:
        a = a + a
        y1 = set_high_word(0.0, hb)
        y2 = b - y1
        t1 = set_high_word(0.0, ha + 0x00100000)
        t2 = a - t1
        w = math.sqrt(t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)))
        if k != 0:
            return set_high_word(0.0, ((1023 + k) << 20)) * w
        return w

# パフォーマンス比較関数
def performance_comparison():
    test_cases = [
        (3.0, 4.0),
        (1e-8, 1e-8),
        (1e300, 1e300),
        (1e-300, 1e-300),
        (1e308, 1e-308),
        (float('inf'), 1.0),
        (0.0, 0.0)
    ]

    print("=== パフォーマンス測定結果 ===")
    for x, y in test_cases:
        print(f"\nTest case: x={x:e}, y={y:e}")
        
        # 時間計測
        time_custom = timeit.timeit(lambda: hypot(x, y), number=100000)
        time_simple = timeit.timeit(lambda: simple_hypot(x, y), number=100000)
        time_standard = timeit.timeit(lambda: math.hypot(x, y), number=100000)
        
        print(f"Custom hypot:   {time_custom:.6f} 秒")
        print(f"Simple hypot:   {time_simple:.6f} 秒")
        print(f"Standard hypot: {time_standard:.6f} 秒")

# 実行
if __name__ == "__main__":
    performance_comparison()

hypot2.go

package main

import (
	"fmt"
	"math"
	"math/big"
	"testing"
	"time"
)

// fdlibm_hypot の実装（ここに実際の実装を入れてください）
func fdlibm_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
}

// SPLIT定数の定義
const SPLIT = 0x1p27 + 1

// sq は二乗計算を行い、高精度部分と低精度部分を返します
func sq(x float64) (hi, lo float64) {
	// xを分割して計算
	xc := x * SPLIT
	xh := x - xc + xc
	xl := x - xh
	
	// 高精度部分と低精度部分を計算
	hi = x * x
	lo = xh*xh - hi + 2*xh*xl + xl*xl
	
	return hi, lo
}

// Hypot は2つの数値の平方和の平方根を計算します
func musl_hypot(x, y float64) float64 {
	// ビット操作のために、float64をuint64として扱う
	ux := math.Float64bits(x)
	uy := math.Float64bits(y)
	
	// 絶対値を取る（符号ビットをクリア）
	ux &= (1<<63 - 1)
	uy &= (1<<63 - 1)
	
	// |x| >= |y| となるように並び替え
	if ux < uy {
		ux, uy = uy, ux
	}
	
	// 指数部を取得
	ex := int(ux >> 52)
	ey := int(uy >> 52)
	
	// float64に戻す
	x = math.Float64frombits(ux)
	y = math.Float64frombits(uy)
	
	// 特殊なケースの処理
	switch {
	case ey == 0x7ff: // y が inf/nan
		return y
	case ex == 0x7ff || uy == 0: // x が inf/nan または y が 0
		return x
	case ex-ey > 64: // y が x に比べて非常に小さい
		return x + y
	}
	
	// オーバーフロー/アンダーフロー対策のスケーリング
	var z float64 = 1
	switch {
	case ex > 0x3ff+510:
		z = 0x1p700
		x *= 0x1p-700
		y *= 0x1p-700
	case ey < 0x3ff-450:
		z = 0x1p-700
		x *= 0x1p700
		y *= 0x1p700
	}
	
	// 二乗計算
	hx, lx := sq(x)
	hy, ly := sq(y)
	
	// 最終結果の計算
	return z * math.Sqrt(ly+lx+hy+hx)
}

// referenceHypot は高精度なhypot計算の参照実装です
func referenceHypot(x, y float64, prec uint) float64 {
    // 高精度の計算のためにbig.Floatを設定
    xBig := new(big.Float).SetPrec(prec).SetFloat64(x)
    yBig := new(big.Float).SetPrec(prec).SetFloat64(y)
    
    // x^2を計算
    x2 := new(big.Float).Mul(xBig, xBig)
    
    // y^2を計算
    y2 := new(big.Float).Mul(yBig, yBig)
    
    // x^2 + y^2を計算
    sum := new(big.Float).Add(x2, y2)
    
    // 平方根を計算
    result := new(big.Float).SetPrec(prec)
    result.Sqrt(sum)
    
    // float64に変換して返す
    f64, _ := result.Float64()
    return f64
}

// テストケースの構造体
type TestCase struct {
	x, y, expected float64
}

// 精度評価用の関数を修正
func evaluateAccuracyWithReference(impl func(x, y float64) float64, testCases []TestCase) (maxError float64, avgError float64) {
    const prec = 256 // 256ビットの精度を使用
    
    var sumError float64
    maxError = 0

    for _, tc := range testCases {
        result := impl(tc.x, tc.y)
        expected := referenceHypot(tc.x, tc.y, prec)
        
        error := math.Abs((result - expected) / expected)
        
        if error > maxError {
            maxError = error
        }
        sumError += error
    }

    avgError = sumError / float64(len(testCases))
    return maxError, avgError
}

// ベンチマーク用の関数
func BenchmarkHypot(b *testing.B, impl func(x, y float64) float64, testCases []TestCase) time.Duration {
	start := time.Now()
	for i := 0; i < b.N; i++ {
		tc := testCases[i%len(testCases)]
		impl(tc.x, tc.y)
	}
	return time.Since(start)
}

func main() {
    // テストケースの準備（既存のものを使用）
    testCases := []TestCase{
        {3, 4, 5},                              // 基本的な三角形
        {1e-308, 1e-308, math.Sqrt(2) * 1e-308}, // 極小値
        {1e308, 1e308, math.Sqrt(2) * 1e308},    // 極大値
        {1e100, 1e-100, 1e100},                  // 大きな差の値
        {123.456, 789.012, 798.6239742181658},   // 一般的な値
        {0, 1, 1},                               // ゼロを含む値
        {math.MaxFloat64 / 2, math.MaxFloat64 / 2, math.MaxFloat64 / math.Sqrt2}, // 大きな値
    }

    // 高精度参照実装を使用した精度評価
    fdlibmMaxErr, fdlibmAvgErr := evaluateAccuracyWithReference(fdlibm_hypot, testCases)
    muslMaxErr, muslAvgErr := evaluateAccuracyWithReference(musl_hypot, testCases)

    fmt.Println("高精度参照実装との比較による精度評価結果:")
    fmt.Printf("fdlibm_hypot - 最大相対誤差: %e, 平均相対誤差: %e\n", fdlibmMaxErr, fdlibmAvgErr)
    fmt.Printf("musl_hypot   - 最大相対誤差: %e, 平均相対誤差: %e\n", muslMaxErr, muslAvgErr)

    // 既存のベンチマーク評価も実行
    b := &testing.B{N: 1000000}
    fdlibmDuration := BenchmarkHypot(b, fdlibm_hypot, testCases)
    muslDuration := BenchmarkHypot(b, musl_hypot, testCases)

    fmt.Println("\n速度評価結果:")
    fmt.Printf("fdlibm_hypot - 実行時間: %v, 1回あたり: %v\n", 
        fdlibmDuration, 
        fdlibmDuration/time.Duration(b.N))
    fmt.Printf("musl_hypot   - 実行時間: %v, 1回あたり: %v\n", 
        muslDuration, 
        muslDuration/time.Duration(b.N))
}

hypot.js