Wikibooks

利用者:Ef3/Hypot2

< 利用者:Ef3
#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <float.h>

#define SPLIT (0x1.0p27 + 1)  // 16進浮動小数点リテラルを使用

// sq関数: 高精度の平方計算
static void sq(double x, double *hi, double *lo) {
    double xc = x * SPLIT;
    double xh = x - xc + xc;
    double xl = x - xh;
    
    *hi = x * x;
    *lo = xh * xh - *hi + 2.0 * xh * xl + xl * xl;
}

// hypot関数: 高精度の2つの値のユークリッド距離
double hypot_custom(double x, double y) {
    // x と y を符号なし64ビット整数に変換
    uint64_t ux = *(uint64_t*)&x, uy  = *(uint64_t*)&y;

    // 絶対値を取得
    ux &= ~(1ULL << 63);  // 符号ビットをクリア
    uy &= ~(1ULL << 63);  // 符号ビットをクリア

    // |x| >= |y| を確認
    if (ux < uy) {
        uint64_t tmp = ux;
        ux = uy;
        uy = tmp;
    }

    // 指数部を取り出す
    int ex = ux >> 52;
    int ey = uy >> 52;

    // 浮動小数点数に戻す
    x = *(double*)&ux;
    y = *(double*)&uy;

    // 特殊ケース処理
    if (ey == 0x7FF) {
        return y;  // y が NaN または Inf の場合
    }
    if (ex == 0x7FF || uy == 0.0) {
        return x;  // x が NaN または Inf または y が 0 の場合
    }
    if (ex - ey > 64) {
        return x + y;  // 大きな差がある場合、単純に x + y を返す
    }

    // 精密な平方根計算の準備
    double z = 1.0;
    if (ex > 0x3FF + 510) {
        z = 0x1.0p700;  // 2^700
        x *= 0x1.0p-700;  // x と y をスケーリング
        y *= 0x1.0p-700;
    } else if (ey < 0x3FF - 450) {
        z = 0x1.0p-700;  // 2^-700
        x *= 0x1.0p700;  // x と y をスケーリング
        y *= 0x1.0p700;
    }

    // 高精度の平方計算
    double hx, lx, hy, ly;
    sq(x, &hx, &lx);
    sq(y, &hy, &ly);

    // 最終的な平方根の計算
    return z * sqrt(lx + ly + hx + hy);
}

int main() {
    printf("Hypotenuse: %1$g(%1$a) %2$g(%2$a)\n", hypot_custom(3.0, 4.0), hypot(3,4));
    printf("Hypotenuse: %1$g(%1$a) %2$g(%2$a)\n", hypot_custom(1e300, 1e300), hypot(1e300, 1e300));
    printf("Hypotenuse: %1$g(%1$a) %2$g(%2$a)\n", hypot_custom(1e-300, 1e-300), hypot(1e-300, 1e-300));
    printf("Hypotenuse: %1$g(%1$a) %2$g(%2$a)\n", hypot_custom(1e-300, 1e300), hypot(1e-300, 1e300));
    printf("Hypotenuse: %1$g(%1$a) %2$g(%2$a)\n", hypot_custom(1e300, 1e-300), hypot(1e300, 1e-300));

    return 0;
}

const SPLIT = Math.pow(2, 27) + 1; // FLT_EVAL_METHODに基づくSPLITの選択

function sq(x) {
  const xc = x * SPLIT;
  const xh = x - xc + xc;
  const xl = x - xh;
  const hi = x * x;
  const lo = xh * xh - hi + 2 * xh * xl + xl * xl;
  return { hi, lo };
}

function hypot(x, y) {
  // 絶対値を取得
  x = Math.abs(x);
  y = Math.abs(y);

  // |x| >= |y| となるように並べ替え
  if (x < y) {
    [x, y] = [y, x];
  }

  const ex = Math.floor(Math.log2(x));
  const ey = Math.floor(Math.log2(y));

  // 特殊ケース処理
  if (!isFinite(y)) return y;
  if (!isFinite(x) || y === 0) return x;
  if (ex - ey > 64) return x + (y * y) / (2 * x);

  // 精密な平方根計算の準備
  let z = 1;
  if (ex > 510) {
    z = Math.pow(2, 700);
    x *= Math.pow(2, -700);
    y *= Math.pow(2, -700);
  } else if (ey < -450) {
    z = Math.pow(2, -700);
    x *= Math.pow(2, 700);
    y *= Math.pow(2, 700);
  }

  const { hi: hx, lo: lx } = sq(x);
  const { hi: hy, lo: ly } = sq(y);
  return z * Math.sqrt(ly + lx + hy + hx);
}

// テスト
console.log(hypot(3, 4)); // 5.0
console.log(hypot(1e-300, 1e-300));
console.log(hypot(1e300, 1e300));

// テストコード
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);

const SPLIT: f64 = (1u64 << 27) as f64 + 1.0; // SPLIT の定義

fn sq(x: f64) -> (f64, f64) {
    let xc = x * SPLIT;
    let xh = x - xc + xc;
    let xl = x - xh;
    let hi = x * x;
    let lo = xh * xh - hi + 2.0 * xh * xl + xl * xl;
    (hi, lo)
}

fn hypot(x: f64, y: f64) -> f64 {
    // 絶対値を取得
    let mut x = x.abs();
    let mut y = y.abs();

    // |x| >= |y| となるように並べ替え
    if x < y {
        std::mem::swap(&mut x, &mut y);
    }

    let ex = x.to_bits() >> 52 & 0x7FF; // x の指数部
    let ey = y.to_bits() >> 52 & 0x7FF; // y の指数部

    // 特殊ケース処理
    if ey == 0x7FF {
        return y; // y が NaN または inf の場合
    }
    if ex == 0x7FF || y == 0.0 {
        return x; // x が inf または y が 0 の場合
    }
    if ex as i32 - ey as i32 > 64 {
        return x + y * y / (2.0 * x); // |x| >> |y| の場合
    }

    // 精密な平方根計算の準備
    let mut z = 1.0;
    if ex > 1023 + 510 {
        z = 2f64.powi(700);
        x *= 2f64.powi(-700);
        y *= 2f64.powi(-700);
    } else if ey < 1023 - 450 {
        z = 2f64.powi(-700);
        x *= 2f64.powi(700);
        y *= 2f64.powi(700);
    }

    let (hx, lx) = sq(x);
    let (hy, ly) = sq(y);
    z * (lx + ly + hx + hy).sqrt()
}

fn main() {
    let x = 3.0;
    let y = 4.0;
    println!("{}", hypot(x, y)); // 5.0
}

#include <iostream>
#include <cmath>
#include <cstdint>
#include <cstring> // 追加: memcpyを使うために必要

// FLT_EVAL_METHOD に基づく SPLIT の定義
#if FLT_EVAL_METHOD > 1 && LDBL_MANT_DIG == 64
#define SPLIT (0x1p32 + 1)
#else
#define SPLIT (0x1p27 + 1)
#endif

// sq関数
void sq(double x, double& hi, double& lo) {
    double xc = x * SPLIT;
    double xh = x - xc + xc;
    double xl = x - xh;
    hi = x * x;
    lo = xh * xh - hi + 2.0 * xh * xl + xl * xl;
}

// hypot関数
double hypot(double x, double y) {
    uint64_t ux, uy;
    std::memcpy(&ux, &x, sizeof(x)); // 修正: memcpyはstd::ではなくグローバル名前空間にある
    std::memcpy(&uy, &y, sizeof(y)); // 修正: 同様

    // 絶対値を取得
    ux &= ~((uint64_t)1 << 63);
    uy &= ~((uint64_t)1 << 63);

    if (ux < uy) {
        std::swap(ux, uy);
    }

    int ex = (ux >> 52) & 0x7FF;
    int ey = (uy >> 52) & 0x7FF;

    x = *reinterpret_cast<double*>(&ux);
    y = *reinterpret_cast<double*>(&uy);

    // 特殊ケース処理
    if (ey == 0x7FF) {
        return y;
    }
    if (ex == 0x7FF || y == 0.0) {
        return x;
    }
    if (ex - ey > 64) {
        return x + y;
    }

    // 精密な平方根計算の準備
    double z = 1.0;
    if (ex > 1023 + 510) {
        z = std::pow(2.0, 700);
        x *= std::pow(2.0, -700);
        y *= std::pow(2.0, -700);
    } else if (ey < 1023 - 450) {
        z = std::pow(2.0, -700);
        x *= std::pow(2.0, 700);
        y *= std::pow(2.0, 700);
    }

    double hx, lx, hy, ly;
    sq(x, hx, lx);
    sq(y, hy, ly);
    return z * std::sqrt(lx + ly + hx + hy);
}

int main() {
    double x = 3.0;
    double y = 4.0;
    std::cout << hypot(x, y) << std::endl; // 5.0
    return 0;
}

import struct
import math

# FLT_EVAL_METHOD に基づく SPLIT の定義
SPLIT = (2 ** 27 + 1)  # それに相当する値を定義

# sq関数
def sq(x):
    xc = x * SPLIT
    xh = x - xc + xc
    xl = x - xh
    hi = x * x
    lo = xh * xh - hi + 2.0 * xh * xl + xl * xl
    return hi, lo

# hypot関数
def hypot(x, y):
    # x と y のビットパターンを整数に変換する
    ux = struct.unpack('>Q', struct.pack('>d', x))[0]
    uy = struct.unpack('>Q', struct.pack('>d', y))[0]

    # 絶対値を取得
    ux &= ~((1 << 63))  # 符号ビットをクリア
    uy &= ~((1 << 63))  # 同様

    # |x| >= |y| を確認
    if ux < uy:
        ux, uy = uy, ux

    # 指数部の取得
    ex = (ux >> 52) & 0x7FF
    ey = (uy >> 52) & 0x7FF

    x = struct.unpack('>d', struct.pack('>Q', ux))[0]
    y = struct.unpack('>d', struct.pack('>Q', uy))[0]

    # 特殊ケース処理
    if ey == 0x7FF:
        return y
    if ex == 0x7FF or y == 0.0:
        return x
    if ex - ey > 64:
        return x + y

    # 精密な平方根計算の準備
    z = 1.0
    if ex > 1023 + 510:
        z = 2 ** 700
        x *= 2 ** -700
        y *= 2 ** -700
    elif ey < 1023 - 450:
        z = 2 ** -700
        x *= 2 ** 700
        y *= 2 ** 700

    # 平方根の計算
    hx, lx = sq(x)
    hy, ly = sq(y)

    return z * math.sqrt(lx + ly + hx + hy)

# 実行例
x = 3.0
y = 4.0
print(hypot(x, y))  # 結果: 5.0
"https://ja.wikibooks.org/w/index.php?title=利用者:Ef3/Hypot2&oldid=263972" より作成