# Python でビット列を解読する
 
import struct
 
def decode_float32(hex_str):
    """32ビット16進数文字列からIEEE 754の各要素を解読"""
    n = int(hex_str, 16)
    sign = (n >> 31) & 1
    exponent = (n >> 23) & 0xFF
    mantissa = n & 0x7FFFFF
 
    print(f"16進数:    {hex_str}")
    print(f"2進数:     {n:032b}")
    print(f"符号(S):   {sign} ({'負' if sign else '正'})")
    print(f"指数部(E): {exponent} (= {exponent} - 127 = {exponent - 127})")
    print(f"仮数部(M): {mantissa:023b}")
 
    if exponent == 0 and mantissa == 0:
        value = 0.0 * (-1 if sign else 1)
        print(f"分類: ゼロ ({'+' if not sign else '-'}0)")
    elif exponent == 0:
        value = (-1)**sign * (mantissa / 2**23) * 2**(-126)
        print(f"分類: 非正規化数")
    elif exponent == 255 and mantissa == 0:
        value = float('inf') * (-1 if sign else 1)
        print(f"分類: {'負' if sign else '正'}の無限大")
    elif exponent == 255:
        value = float('nan')
        print(f"分類: NaN")
    else:
        value = (-1)**sign * (1 + mantissa / 2**23) * 2**(exponent - 127)
        print(f"分類: 正規化数")
 
    print(f"値:        {value}")
    return value
 
# 使用例
decode_float32("40D00000")  # → 6.5
decode_float32("C1460000")  # → -12.375
decode_float32("3DCCCCCD")  # → 0.10000000149011612（0.1の近似）

/* C言語でのビット列解読 */
#include <stdio.h>
#include <stdint.h>
#include <string.h>
 
void decode_float32(float f) {
    uint32_t bits;
    memcpy(&bits, &f, sizeof(bits));  /* 型パンニング（安全な方法） */
 
    uint32_t sign     = (bits >> 31) & 1;
    uint32_t exponent = (bits >> 23) & 0xFF;
    uint32_t mantissa = bits & 0x7FFFFF;
 
    printf("値:       %.*g\n", 9, f);
    printf("ビット列: ");
    for (int i = 31; i >= 0; i--) {
        printf("%d", (bits >> i) & 1);
        if (i == 31 || i == 23) printf(" ");
    }
    printf("\n");
    printf("符号:     %u (%s)\n", sign, sign ? "負" : "正");
    printf("指数部:   %u (実際の指数 = %d)\n", exponent, (int)exponent - 127);
    printf("仮数部:   0x%06X\n", mantissa);
 
    if (exponent == 0 && mantissa == 0)
        printf("分類: %sゼロ\n", sign ? "負の" : "正の");
    else if (exponent == 0)
        printf("分類: 非正規化数\n");
    else if (exponent == 255 && mantissa == 0)
        printf("分類: %s無限大\n", sign ? "負の" : "正の");
    else if (exponent == 255)
        printf("分類: NaN\n");
    else
        printf("分類: 正規化数\n");
}
 
int main(void) {
    decode_float32(6.5f);
    decode_float32(-12.375f);
    decode_float32(0.1f);
    return 0;
}

# Python で符号付きゼロを確認
import math
 
pos_zero = +0.0
neg_zero = -0.0
 
print(pos_zero == neg_zero)         # True
print(math.copysign(1, pos_zero))   # 1.0
print(math.copysign(1, neg_zero))   # -1.0
 
# 符号付きゼロの検出
def is_negative_zero(x):
    return x == 0.0 and math.copysign(1, x) < 0
 
print(is_negative_zero(-0.0))  # True
print(is_negative_zero(+0.0))  # False

import math
 
inf = float('inf')
 
# 無限大の生成と演算
print(1.0 / 0.0)       # inf（Python ではデフォルトで例外なし...ではなく ZeroDivisionError）
# 注意: Python では 1.0/0.0 は ZeroDivisionError
# float('inf') で直接生成する
 
print(inf + inf)        # inf
print(inf + 1e308)      # inf
print(inf * -1)         # -inf
print(inf * 0)          # nan
print(inf - inf)        # nan
print(inf / inf)        # nan
print(1.0 / inf)        # 0.0
print(inf > 1e308)      # True
 
# 無限大の判定
print(math.isinf(inf))          # True
print(math.isinf(-inf))         # True
print(math.isinf(1e308))        # False
print(math.isinf(1e308 * 10))   # True（オーバーフロー）

import math
import numpy as np
 
x = float('nan')
 
# NaN の根本的性質: 自分自身と等しくない
print(x == x)     # False  ← IEEE 754 で規定された動作
print(x != x)     # True
print(x > 0)      # False
print(x < 0)      # False
print(x >= 0)     # False
print(x <= 0)     # False
# → NaN との比較は != 以外すべて False
 
# NaN の判定方法
print(math.isnan(x))            # True  ← 推奨
print(x != x)                    # True  ← 伝統的イディオム（非推奨）
print(np.isnan(x))              # True  ← NumPy
 
# NaN の伝播（「毒」のように広がる）
print(x + 1)          # nan
print(x * 0)          # nan
print(x ** 0)         # 1.0  ← 例外! IEEE 754 で規定
print(0 * float('inf'))  # nan
print(max(x, 5))      # nan（Python 標準）
print(min(x, 5))      # nan
 
# NaN を含むリストの集約
values = [1.0, 2.0, float('nan'), 4.0]
print(sum(values))     # nan（1つでも NaN があると結果が NaN）
print(max(values))     # nan
print(min(values))     # nan
 
# NaN 安全な集約（NumPy）
arr = np.array(values)
print(np.nansum(arr))    # 7.0（NaN を無視して合計）
print(np.nanmean(arr))   # 2.333...（NaN を除外して平均）
print(np.nanmax(arr))    # 4.0

# 0.1 + 0.2 問題の詳細な確認
 
from decimal import Decimal
 
# float の正確な値を Decimal で確認
print(f"0.1 の格納値: {Decimal(0.1)}")
print(f"0.2 の格納値: {Decimal(0.2)}")
print(f"0.3 の格納値: {Decimal(0.3)}")
print(f"0.1+0.2 の値: {Decimal(0.1) + Decimal(0.2)}")
print()
print(f"差分: {(Decimal(0.1) + Decimal(0.2)) - Decimal(0.3)}")
print(f"   = {float(0.1) + float(0.2) - float(0.3)}")
 
# 出力:
# 0.1 の格納値: 0.1000000000000000055511151231257827021181583404541015625
# 0.2 の格納値: 0.200000000000000011102230246251565404236316680908203125
# 0.3 の格納値: 0.299999999999999988897769753748434595763683319091796875
# 0.1+0.2 の値: 0.3000000000000000166533453693773481063544750213623046875
# 差分: 1.77635683940025046E-17
 
# 主要言語での結果
# Python:     0.1 + 0.2 == 0.30000000000000004
# JavaScript: 0.1 + 0.2 === 0.30000000000000004
# C/C++:      0.1 + 0.2 == 0.30000000000000004
# Java:       0.1 + 0.2 == 0.30000000000000004
# Ruby:       0.1 + 0.2 == 0.30000000000000004
# → 全ての言語で同じ結果（IEEE 754 準拠のため）

# Python での丸めモードの確認
from decimal import Decimal, ROUND_HALF_UP, ROUND_HALF_EVEN
 
# Python の組み込み round() は偶数丸め
print(round(0.5))   # 0  ← 偶数丸め
print(round(1.5))   # 2  ← 偶数丸め
print(round(2.5))   # 2  ← 偶数丸め
print(round(3.5))   # 4  ← 偶数丸め
 
# Decimal で明示的に丸めモードを指定
d = Decimal('2.5')
print(d.quantize(Decimal('1'), rounding=ROUND_HALF_UP))    # 3（四捨五入）
print(d.quantize(Decimal('1'), rounding=ROUND_HALF_EVEN))  # 2（偶数丸め）

# 桁落ちの典型例: 2次方程式の解の公式
 
import math
 
def quadratic_naive(a, b, c):
    """素朴な解の公式（桁落ちが発生しうる）"""
    discriminant = b*b - 4*a*c
    sqrt_d = math.sqrt(discriminant)
    x1 = (-b + sqrt_d) / (2 * a)
    x2 = (-b - sqrt_d) / (2 * a)
    return x1, x2
 
def quadratic_stable(a, b, c):
    """桁落ちを回避する安定版"""
    discriminant = b*b - 4*a*c
    sqrt_d = math.sqrt(discriminant)
 
    # b の符号に応じて桁落ちしない方の解を先に計算
    if b >= 0:
        q = -0.5 * (b + sqrt_d)
    else:
        q = -0.5 * (b - sqrt_d)
 
    x1 = q / a
    x2 = c / q  # ビエタの公式: x1 * x2 = c/a を利用
    return x1, x2
 
# テスト: a=1, b=10^8, c=1 → 真の解は x ≈ -10^(-8), x ≈ -10^8
a, b, c = 1, 1e8, 1
 
naive = quadratic_naive(a, b, c)
stable = quadratic_stable(a, b, c)
 
print(f"素朴な解:   x1 = {naive[0]:.15e}, x2 = {naive[1]:.15e}")
print(f"安定な解:   x1 = {stable[0]:.15e}, x2 = {stable[1]:.15e}")
print(f"理論値:     x1 ≈ -1e-08,             x2 = -1e+08")
 
# 素朴な解: x1 = -7.450580596923828e-09 ← 誤差大
# 安定な解: x1 = -1.000000000000000e-08 ← 正確

# 情報落ちの確認
 
# 演算順序による結果の違い
a = 1e16
b = 1.0
 
print(f"(a + b) - a = {(a + b) - a}")   # 0.0（b の情報が消失）
print(f"(a - a) + b = {(a - a) + b}")   # 1.0（正しい結果）
 
# より深刻な例: 大量の小さな値を大きな値に加算
big = 1e15
result_forward = big
for i in range(1000000):
    result_forward += 1.0
 
result_reverse = 0.0
for i in range(1000000):
    result_reverse += 1.0
result_reverse += big
 
print(f"前方加算: {result_forward}")      # 精度が低い
print(f"逆方向:   {result_reverse}")      # やや正確
print(f"理論値:   {big + 1000000.0}")

# 丸め誤差の蓄積: 0.1 を 1000 回加算
 
# 素朴な加算
total_naive = 0.0
for i in range(1000):
    total_naive += 0.1
print(f"素朴な加算:     {total_naive}")         # 99.99999999999986
print(f"誤差:           {total_naive - 100.0}")  # -1.4e-13
 
# math.fsum（内部で拡張精度を使用）
import math
total_fsum = math.fsum([0.1] * 1000)
print(f"math.fsum:      {total_fsum}")           # 100.00000000000007
print(f"誤差:           {total_fsum - 100.0}")    # 7.1e-14
 
# Kahan 補償加算アルゴリズム
def kahan_sum(values):
    """Kahan の補償加算: 丸め誤差を補正項で追跡"""
    total = 0.0
    compensation = 0.0  # 誤差の蓄積を追跡する補正項
    for value in values:
        y = value - compensation      # 補正を適用
        t = total + y                 # 加算（ここで丸め誤差が発生）
        compensation = (t - total) - y  # 丸め誤差を捕捉
        total = t
    return total
 
total_kahan = kahan_sum([0.1] * 1000)
print(f"Kahan加算:      {total_kahan}")          # 100.00000000000007
print(f"誤差:           {total_kahan - 100.0}")   # 非常に小さい
 
# Neumaier の改良版（Kahan より頑健）
def neumaier_sum(values):
    """Neumaier の補償加算: |total| < |value| の場合も正しく補正"""
    total = 0.0
    compensation = 0.0
    for value in values:
        t = total + value
        if abs(total) >= abs(value):
            compensation += (total - t) + value
        else:
            compensation += (value - t) + total
        total = t
    return total + compensation
 
total_neumaier = neumaier_sum([0.1] * 1000)
print(f"Neumaier加算:   {total_neumaier}")

# ULP と整数精度の限界を確認
 
import numpy as np
 
# binary32 の整数精度の限界
f32 = np.float32
 
print(f"2^23     = {f32(2**23)}")              # 8388608.0
print(f"2^23 + 1 = {f32(2**23 + 1)}")          # 8388609.0（正確）
print(f"2^24     = {f32(2**24)}")              # 16777216.0
print(f"2^24 + 1 = {f32(2**24 + 1)}")          # 16777216.0 ← 同じ値!
print(f"2^24 + 2 = {f32(2**24 + 2)}")          # 16777218.0
print(f"2^24 + 3 = {f32(2**24 + 3)}")          # 16777220.0 ← +4 されている!
 
print()
 
# binary64 の整数精度の限界
print(f"2^53     = {float(2**53)}")             # 9007199254740992.0
print(f"2^53 + 1 = {float(2**53 + 1)}")         # 9007199254740992.0 ← 同じ!
 
# JavaScript での影響
# JSON.parse('{"id": 9007199254740993}')
# → {"id": 9007199254740992}  ← ID が変わってしまう!
# → Twitter が snowflake ID を文字列で返すのはこのため
 
# ULP の計算
def ulp(x):
    """与えられた値における ULP を計算"""
    return np.spacing(x)
 
for val in [0.5, 1.0, 1000.0, 1e6, 1e15]:
    print(f"ULP({val:>10}) = {ulp(val):.6e}")

# --- アンチパターン: == による浮動小数点比較 ---
 
# 危険なコード
total = 0.0
for _ in range(10):
    total += 0.1
 
if total == 1.0:       # ← 永遠に True にならない可能性!
    print("合計は1.0")
else:
    print(f"合計は {total}")  # 合計は 0.9999999999999999
 
# ループの終了条件での危険
x = 0.0
while x != 1.0:       # ← 無限ループの危険!
    x += 0.1
    if x > 2.0:       # 安全弁がないと本当に無限ループ
        break
 
# --- 正しいパターン ---
 
import math
 
# パターン1: math.isclose()（Python 3.5+）
print(math.isclose(0.1 + 0.2, 0.3))  # True
# デフォルト: rel_tol=1e-9, abs_tol=0.0
 
# パターン2: 相対誤差による比較
def nearly_equal(a, b, rel_tol=1e-9, abs_tol=1e-12):
    """相対誤差と絶対誤差の両方を考慮した比較"""
    if a == b:  # Inf == Inf, +0 == -0 を正しく扱う
        return True
    if math.isnan(a) or math.isnan(b):
        return False
    diff = abs(a - b)
    return diff <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
 
# パターン3: ループでは < や > を使う
x = 0.0
while x < 1.0:  # != ではなく < を使用
    x += 0.1

/* C言語での浮動小数点比較 */
#include <math.h>
#include <float.h>
#include <stdbool.h>
 
/* アンチパターン */
bool bad_compare(double a, double b) {
    return a == b;  /* 浮動小数点では信頼できない */
}
 
/* 正しいパターン: 相対・絶対誤差の組み合わせ */
bool nearly_equal(double a, double b, double rel_tol, double abs_tol) {
    if (a == b) return true;  /* Inf, 0 の扱い */
    if (isnan(a) || isnan(b)) return false;
 
    double diff = fabs(a - b);
    double larger = fmax(fabs(a), fabs(b));
 
    return diff <= fmax(rel_tol * larger, abs_tol);
}
 
/* 使用例 */
int main(void) {
    double x = 0.1 + 0.2;
    double y = 0.3;
 
    /* NG */
    if (x == y) { /* 到達しない */ }
 
    /* OK */
    if (nearly_equal(x, y, 1e-9, 1e-12)) {
        /* 正しく到達する */
    }
 
    return 0;
}

# --- アンチパターン: 金融計算に float を使用 ---
 
# 危険なコード
price = 19.99
tax_rate = 0.08
tax = price * tax_rate          # 1.5992000000000002
total = price + tax             # 21.5892
print(f"税込: ${total:.2f}")    # $21.59（表示上は正しく見えるが...）
 
# 大量の取引で誤差が蓄積
daily_amounts = [0.01] * 1000000  # 100万件の1セント取引
total = sum(daily_amounts)
print(f"合計: ${total:.2f}")  # $10000.00 にならない可能性
 
# --- 正しいパターン1: Decimal 型 ---
from decimal import Decimal, ROUND_HALF_UP, getcontext
 
# 精度を設定
getcontext().prec = 28
 
price = Decimal('19.99')       # 文字列から生成（float 経由しない!）
tax_rate = Decimal('0.08')
tax = (price * tax_rate).quantize(
    Decimal('0.01'),
    rounding=ROUND_HALF_UP
)
total = price + tax
print(f"税込: ${total}")       # $21.59（正確）
 
# 大量の取引でも正確
daily_amounts = [Decimal('0.01')] * 1000000
total = sum(daily_amounts)
print(f"合計: ${total}")       # $10000.00（正確）
 
# --- 正しいパターン2: 整数演算（セント単位） ---
price_cents = 1999             # $19.99 = 1999 セント
tax_rate_bps = 800             # 8% = 800 ベーシスポイント
tax_cents = (price_cents * tax_rate_bps + 5000) // 10000  # 四捨五入
total_cents = price_cents + tax_cents
 
print(f"税込: ${total_cents / 100:.2f}")  # $21.59（正確）
 
# 注意: Decimal('0.1') と Decimal(0.1) は異なる!
print(Decimal('0.1'))   # 0.1（正確）
print(Decimal(0.1))     # 0.1000000000000000055511151231257827...（float経由の誤差）

# 落とし穴1: 演算の順序依存性（結合法則の不成立）
a, b, c = 1e20, -1e20, 1.0
print(f"(a + b) + c = {(a + b) + c}")   # 1.0
print(f"a + (b + c) = {a + (b + c)}")   # 0.0  ← 異なる結果!
 
# 落とし穴2: 分配法則の不成立
a, b, c = 1e15, 1.0, -1e15
print(f"a × (b + c) = {a * (b + c)}")     # 期待: a × 1.0 - a × 1e15... 複雑
# 一般に a*(b+c) ≠ a*b + a*c
 
# 落とし穴3: 比較の非推移性
# a < b かつ b < c でも a < c とは限らない（NaN がある場合）
a, b, c = 1.0, float('nan'), 2.0
print(f"a < b: {a < b}")  # False
print(f"b < c: {b < c}")  # False
print(f"a < c: {a < c}")  # True  ← NaN が比較を壊す
 
# 落とし穴4: ソートの不安定性
import random
values = [1.0, float('nan'), 2.0, float('nan'), 0.5]
# sorted(values) → NaN の位置が不定、ソートが壊れる可能性
 
# 落とし穴5: ハッシュの一貫性
# Python では hash(0) == hash(0.0) == hash(Decimal('0'))
# しかし hash(float('nan')) は呼び出すごとに同じだが、
# NaN == NaN が False なので dict のキーとして使うと問題が起こる
 
# 落とし穴6: 型変換の罠（JavaScript）
# JSON.parse('{"value": 9007199254740993}') → 9007199254740992
# → 大きな整数 ID が JSON パース時に変化する

import math
from typing import Optional
 
def robust_float_equal(
    a: float,
    b: float,
    rel_tol: float = 1e-9,
    abs_tol: float = 1e-12
) -> bool:
    """
    堅牢な浮動小数点比較関数。
 
    エッジケースを正しく処理:
    - NaN: NaN同士でも False を返す（IEEE 754 準拠）
    - Inf: 同符号の Inf 同士は True
    - -0 と +0: True（IEEE 754 準拠）
    - 非常に小さい値: abs_tol で判定
    - 通常の値: rel_tol で判定
 
    Parameters:
        a, b: 比較する浮動小数点数
        rel_tol: 相対許容誤差（デフォルト 1e-9）
        abs_tol: 絶対許容誤差（デフォルト 1e-12）
 
    Returns:
        a と b が十分に近いかどうか
    """
    # NaN の処理（NaN は何とも等しくない）
    if math.isnan(a) or math.isnan(b):
        return False
 
    # 完全一致（Inf == Inf, +0 == -0 を含む）
    if a == b:
        return True
 
    # 差分の計算
    diff = abs(a - b)
 
    # 無限大の場合（符号が異なる Inf 同士）
    if math.isinf(a) or math.isinf(b):
        return False
 
    # 相対誤差 or 絶対誤差で判定
    larger = max(abs(a), abs(b))
    return diff <= max(rel_tol * larger, abs_tol)
 
 
# テスト
assert robust_float_equal(0.1 + 0.2, 0.3)          # True
assert not robust_float_equal(1.0, 2.0)             # False
assert robust_float_equal(float('inf'), float('inf'))  # True
assert not robust_float_equal(float('inf'), float('-inf'))  # False
assert not robust_float_equal(float('nan'), float('nan'))   # False
assert robust_float_equal(0.0, -0.0)                # True
assert robust_float_equal(1e-15, 1.1e-15, rel_tol=0.1)  # True

# 条件数（Condition Number）による誤差の予測
 
import numpy as np
 
# 条件数の悪い連立方程式
A_bad = np.array([
    [1.0, 1.0],
    [1.0, 1.0001]
])
b = np.array([2.0, 2.0001])
 
cond = np.linalg.cond(A_bad)
print(f"条件数: {cond:.0f}")  # 約40000
 
# 条件数が大きい → 入力のわずかな変化で出力が大きく変動
x = np.linalg.solve(A_bad, b)
print(f"解: {x}")  # [1.0, 1.0]
 
# b をわずかに摂動
b_perturbed = b + np.array([0.0001, 0.0])
x_perturbed = np.linalg.solve(A_bad, b_perturbed)
print(f"摂動解: {x_perturbed}")  # 大きく異なる可能性
 
# 誤差の上界: ||δx||/||x|| ≤ cond(A) × ||δb||/||b||
# → 条件数が 10^4 なら、入力の 10^(-12) の誤差が 10^(-8) の結果誤差になりうる

# PyTorch での混合精度学習の例
import torch
from torch.cuda.amp import autocast, GradScaler
 
model = MyModel().cuda()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
scaler = GradScaler()  # Loss Scaling を自動管理
 
for data, target in dataloader:
    optimizer.zero_grad()
 
    # autocast: 演算ごとに最適な精度を自動選択
    with autocast():
        output = model(data)        # FP16 で順伝播
        loss = loss_fn(output, target)
 
    # Loss Scaling + FP16 逆伝播
    scaler.scale(loss).backward()
 
    # FP32 でパラメータ更新
    scaler.step(optimizer)
    scaler.update()

import struct
import math
 
def visualize_float64(value):
    """float64 のビット表現を詳細に可視化"""
    # float → bytes → int
    raw_bytes = struct.pack('>d', value)
    bits = int.from_bytes(raw_bytes, 'big')
 
    # フィールド抽出
    sign = (bits >> 63) & 1
    exponent = (bits >> 52) & 0x7FF
    mantissa = bits & 0xFFFFFFFFFFFFF
 
    # ビット文字列
    bit_str = f"{bits:064b}"
    formatted = f"{bit_str[0]} {bit_str[1:12]} {bit_str[12:]}"
 
    print(f"値:         {value}")
    print(f"16進数:     {raw_bytes.hex()}")
    print(f"ビット列:   {formatted}")
    print(f"符号(S):    {sign} ({'負' if sign else '正'})")
    print(f"指数部(E):  {exponent} (バイアス除去後: {exponent - 1023})")
    print(f"仮数部(M):  {mantissa:052b}")
 
    # 分類
    if exponent == 0 and mantissa == 0:
        print(f"分類:       {'負' if sign else '正'}のゼロ")
    elif exponent == 0:
        actual_exp = 1 - 1023
        value_calc = (-1)**sign * (mantissa / 2**52) * 2**actual_exp
        print(f"分類:       非正規化数 (実効指数: {actual_exp})")
    elif exponent == 2047 and mantissa == 0:
        print(f"分類:       {'負' if sign else '正'}の無限大")
    elif exponent == 2047:
        snan = (mantissa >> 51) & 1 == 0
        print(f"分類:       {'Signaling' if snan else 'Quiet'} NaN")
    else:
        actual_exp = exponent - 1023
        print(f"分類:       正規化数 (実効指数: {actual_exp})")
        print(f"有効値:     1.{mantissa:052b} × 2^{actual_exp}")
 
    # 隣接値
    if not (math.isnan(value) or math.isinf(value)):
        next_val = struct.unpack('>d', (bits + 1).to_bytes(8, 'big'))[0]
        prev_val = struct.unpack('>d', (bits - 1).to_bytes(8, 'big'))[0]
        print(f"次の値:     {next_val} (差: {next_val - value})")
        print(f"前の値:     {prev_val} (差: {value - prev_val})")
 
    print()
 
# デバッグ例
visualize_float64(0.1)
visualize_float64(0.2)
visualize_float64(0.1 + 0.2)
visualize_float64(0.3)

class Interval:
    """単純な区間演算クラス（デバッグ用）"""
 
    def __init__(self, lo, hi=None):
        if hi is None:
            # 浮動小数点数から区間を生成（丸め誤差を考慮）
            import sys
            eps = sys.float_info.epsilon
            if lo == 0.0:
                self.lo, self.hi = -5e-324, 5e-324
            else:
                ulp = abs(lo) * eps
                self.lo = lo - ulp
                self.hi = lo + ulp
        else:
            self.lo = lo
            self.hi = hi
 
    def __add__(self, other):
        return Interval(self.lo + other.lo, self.hi + other.hi)
 
    def __sub__(self, other):
        return Interval(self.lo - other.hi, self.hi - other.lo)
 
    def __mul__(self, other):
        products = [
            self.lo * other.lo, self.lo * other.hi,
            self.hi * other.lo, self.hi * other.hi
        ]
        return Interval(min(products), max(products))
 
    def __repr__(self):
        mid = (self.lo + self.hi) / 2
        radius = (self.hi - self.lo) / 2
        return f"[{self.lo:.17g}, {self.hi:.17g}] (幅: {self.hi - self.lo:.3e})"
 
    def contains(self, value):
        return self.lo <= value <= self.hi
 
# 0.1 + 0.2 の区間演算
a = Interval(0.1)
b = Interval(0.2)
c = a + b
print(f"0.1 の区間: {a}")
print(f"0.2 の区間: {b}")
print(f"0.1+0.2:    {c}")
print(f"0.3 を含む: {c.contains(0.3)}")

import math
import warnings
 
def check_numerical_stability(func, x, delta=1e-8):
    """関数の数値安定性を簡易チェック"""
    y = func(x)
    y_plus = func(x + delta)
    y_minus = func(x - delta)
 
    # 条件数の近似: |x * f'(x) / f(x)|
    if abs(y) > 0:
        deriv_approx = (y_plus - y_minus) / (2 * delta)
        cond_approx = abs(x * deriv_approx / y)
    else:
        cond_approx = float('inf')
 
    # 対称性チェック
    forward_diff = y_plus - y
    backward_diff = y - y_minus
    if abs(forward_diff) > 0:
        symmetry = abs(forward_diff - backward_diff) / abs(forward_diff)
    else:
        symmetry = 0.0
 
    print(f"f({x}) = {y}")
    print(f"条件数の推定: {cond_approx:.2e}")
    if cond_approx > 1e10:
        warnings.warn(f"条件数が非常に大きい（{cond_approx:.2e}）: 数値不安定の可能性")
    print(f"差分の対称性: {symmetry:.2e}")
    if symmetry > 0.01:
        warnings.warn("差分の対称性が低い: 丸め誤差の影響が大きい可能性")
 
# テスト: 桁落ちが発生する関数
def unstable_func(x):
    """x が 0 に近いとき桁落ちが発生"""
    return (1 - math.cos(x)) / (x * x)  # 理論値は x→0 で 0.5 に収束
 
def stable_func(x):
    """数学的に等価だが数値的に安定"""
    return 2 * (math.sin(x/2) / x) ** 2  # 半角公式を利用
 
print("=== 不安定な実装 ===")
check_numerical_stability(unstable_func, 1e-8)
print()
print("=== 安定な実装 ===")
check_numerical_stability(stable_func, 1e-8)

# Python で浮動小数点例外を検出する
 
import numpy as np
 
# NumPy の浮動小数点例外設定
# デフォルトでは警告のみ。'raise' で例外を発生させられる
old_settings = np.seterr(all='raise')  # 全ての浮動小数点例外で例外発生
 
try:
    result = np.float64(1e308) * np.float64(10)  # オーバーフロー
except FloatingPointError as e:
    print(f"捕捉: {e}")
 
try:
    result = np.float64(0.0) / np.float64(0.0)  # 無効演算
except FloatingPointError as e:
    print(f"捕捉: {e}")
 
np.seterr(**old_settings)  # 設定を元に戻す
 
# warnings モジュールでの検出
import warnings
warnings.filterwarnings('error', category=RuntimeWarning)
 
try:
    result = np.float64(1.0) / np.float64(0.0)
except RuntimeWarning as e:
    print(f"警告を捕捉: {e}")
 
warnings.resetwarnings()

/* C言語での浮動小数点例外のトラップ */
#define _GNU_SOURCE
#include <stdio.h>
#include <fenv.h>
#include <math.h>
#include <signal.h>
 
/* 浮動小数点例外ハンドラ */
void fpe_handler(int sig) {
    printf("浮動小数点例外が発生!\n");
 
    /* 例外フラグを確認 */
    if (fetestexcept(FE_DIVBYZERO))
        printf("  - ゼロ除算\n");
    if (fetestexcept(FE_OVERFLOW))
        printf("  - オーバーフロー\n");
    if (fetestexcept(FE_UNDERFLOW))
        printf("  - アンダーフロー\n");
    if (fetestexcept(FE_INVALID))
        printf("  - 無効演算\n");
    if (fetestexcept(FE_INEXACT))
        printf("  - 不正確（丸めが発生）\n");
 
    feclearexcept(FE_ALL_EXCEPT);
}
 
int main(void) {
    /* 例外フラグを使った事後チェック（推奨） */
    feclearexcept(FE_ALL_EXCEPT);
 
    volatile double a = 1.0;
    volatile double b = 0.0;
    volatile double result = a / b;
 
    if (fetestexcept(FE_DIVBYZERO)) {
        printf("ゼロ除算が発生: result = %f\n", result);
    }
 
    feclearexcept(FE_ALL_EXCEPT);
    result = sqrt(-1.0);
 
    if (fetestexcept(FE_INVALID)) {
        printf("無効演算が発生: result = %f\n", result);
    }
 
    return 0;
}

x = 1e16
print(x + 1 == x)
print(x + 2 == x)

# safe_equal
assert safe_equal(0.1 + 0.2, 0.3)
assert safe_equal(float('inf'), float('inf'))
assert not safe_equal(float('nan'), float('nan'))
assert safe_equal(+0.0, -0.0)
 
# safe_sum
assert abs(safe_sum([0.1] * 10) - 1.0) < 1e-15
 
# safe_mean
assert safe_mean([1.0, float('nan'), 3.0]) == 2.0
assert math.isnan(safe_mean([]))
 
# float_info
info = float_info(6.5)
assert info['sign'] == 0
assert info['exponent'] == 2
assert info['category'] == 'normal'

# 正確に計算される演算の例
 
# 1. 同符号の値の減算（Sterbenz の定理）
a = 1.5
b = 1.0
# a/2 ≤ b ≤ 2a を満たすので a - b は正確
print(a - b)  # 0.5（正確）
 
# 2. 2のべき乗の乗除算
x = 3.14159
print(x * 2.0)     # 6.28318（正確: 指数部の調整のみ）
print(x * 0.5)     # 1.570795（正確: 指数部の調整のみ）
print(x * 4.0)     # 12.56636（正確）
 
# 3. FMA（融合積和演算）
import math
# math.fma(a, b, c) = a*b + c を1回の丸めで計算（Python 3.13+）
# → a*b の中間結果が丸められないため、通常の a*b + c より正確
 
# 4. 二重倍精度（Double-Double）演算
def two_sum(a, b):
    """a + b を高精度に計算。s + e = a + b が正確に成り立つ"""
    s = a + b
    v = s - a
    e = (a - (s - v)) + (b - v)
    return s, e  # s は丸められた和、e は誤差
 
s, e = two_sum(1e16, 1.0)
print(f"和: {s}, 誤差: {e}")  # 和: 1e16, 誤差: 1.0
# → 失われた情報が e に保存されている