Numerical Utilities Example — hypot_safe, prime_factors, cpu_time¶
This example exercises three miscellaneous IMSL utility functions:
utilities.hypot_safe() for overflow-safe hypotenuse computation
(including extreme values such as 10²⁰⁰),
utilities.prime_factors() for integer factorisation, and
utilities.cpu_time() (IMSL CPSEC) for CPU-time measurement.
Machine constants (epsilon, float max, float min) are read from
numpy and floating-point precision is visualised as relative
rounding error vs magnitude, alongside a float32-vs-float64 comparison.
Example Code¶
"""IMSL numerical utilities example: hypot_safe, prime_factors, cpu_time,
and floating-point precision analysis.
Outputs:
- Machine constants and results printed to stdout
- SVG floating-point precision visualisation saved to
test_output/example_imsl_numerical_utils.svg
"""
from __future__ import annotations
import sys
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from utilities import hypot_safe, prime_factors, cpu_time, cpu_time as cpsec
def run_demo_imsl_numerical_utils() -> Dict[str, object]:
"""Demonstrate numerical utility functions and floating-point properties.
Uses :func:`utilities.hypot_safe` for overflow-safe hypotenuse
computation, :func:`utilities.prime_factors` for integer factorisation,
and :func:`utilities.cpu_time` (IMSL CPSEC) for timing. Computes
machine epsilon and float limits from numpy and visualises relative
rounding error as a function of magnitude.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``machine_eps`` (float),
``float_max`` (float), ``float_min`` (float),
``hypot_results`` (list), ``factored`` (dict),
and ``plot_path`` (str).
"""
# ── Machine constants ─────────────────────────────────────────────────────
finfo = np.finfo(np.float64)
machine_eps = float(finfo.eps)
float_max = float(finfo.max)
float_min = float(finfo.tiny)
print("\nIMSL Numerical Utilities Example")
print("=" * 55)
print("\n Machine Constants (float64 / IEEE 754 double)")
print(f" {'Machine epsilon':<28} {machine_eps:.6e}")
print(f" {'Float max':<28} {float_max:.6e}")
print(f" {'Float min (tiny)':<28} {float_min:.6e}")
print(f" {'Mantissa bits':<28} {finfo.nmant}")
print(f" {'Decimal digits':<28} {finfo.precision}")
print(f" {'Python sys.float_info.epsilon':<28} {sys.float_info.epsilon:.6e}")
# ── hypot_safe examples ───────────────────────────────────────────────────
hypot_cases = [
(3.0, 4.0),
(5.0, 12.0),
(1e200, 1e200), # would overflow with naive a^2+b^2
(1e-200, 1e-200), # would underflow with naive approach
(0.0, 7.3),
]
hypot_results: list[dict] = []
print("\n hypot_safe — overflow/underflow-safe hypotenuse")
print(f" {'a':>12} {'b':>12} {'hypot_safe':>16} {'naive ok?':>10}")
print(" " + "-" * 56)
for a, b in hypot_cases:
h = hypot_safe(a, b)
try:
naive = (a**2 + b**2) ** 0.5
naive_ok = "yes" if np.isfinite(naive) else "overflow"
except OverflowError:
naive_ok = "overflow"
hypot_results.append({"a": a, "b": b, "hypot": h})
print(f" {a:>12.3e} {b:>12.3e} {h:>16.6e} {naive_ok:>10}")
# ── prime_factors examples ────────────────────────────────────────────────
test_numbers = [12, 60, 100, 360, 1024, 9999, 97, 1000003]
factored: dict[int, list[int]] = {}
print("\n prime_factors — integer factorisation")
print(f" {'n':>10} {'factors'}")
print(" " + "-" * 40)
for n in test_numbers:
factors = prime_factors(n)
factored[n] = factors
expr = " × ".join(str(f) for f in factors)
print(f" {n:>10d} {expr}")
# ── cpu_time / cpsec timing ───────────────────────────────────────────────
t0 = cpu_time()
_ = np.linalg.eig(np.random.default_rng(0).random((200, 200)))
elapsed = cpu_time() - t0
print(f"\n cpu_time (CPSEC): eig(200×200) took {elapsed*1000:.3f} ms CPU time")
# ── Floating-point precision visualisation ────────────────────────────────
magnitudes = np.logspace(-15, 15, 300)
rel_errors = np.array([
abs((x + machine_eps * x) - x) / x if x != 0 else machine_eps
for x in magnitudes
])
# Round-trip through float32 to show precision loss
mag32 = np.logspace(-6, 6, 200)
x64 = mag32.astype(np.float64)
x32 = mag32.astype(np.float32).astype(np.float64)
rel_err32 = np.abs(x64 - x32) / (np.abs(x64) + 1e-300)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_numerical_utils.svg"
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Panel 1: relative rounding error vs magnitude
ax1 = axes[0]
ax1.loglog(magnitudes, rel_errors, color="#b45309", linewidth=1.8,
label=f"float64 ε={machine_eps:.1e}")
ax1.axhline(machine_eps, color="gray", linewidth=1.0, linestyle="--",
label="Machine epsilon")
ax1.set_xlabel("Magnitude")
ax1.set_ylabel("Relative rounding error")
ax1.set_title("Float64 precision vs magnitude")
ax1.legend(fontsize=8)
ax1.grid(True, which="both", alpha=0.3)
# Panel 2: float32 vs float64 relative error
ax2 = axes[1]
ax2.semilogy(mag32, rel_err32 + 1e-40, color="#d97706", linewidth=1.8,
label="float32 round-trip error")
ax2.axhline(np.finfo(np.float32).eps, color="#92400e", linewidth=1.0,
linestyle="--", label=f"float32 ε={np.finfo(np.float32).eps:.1e}")
ax2.set_xlabel("Magnitude")
ax2.set_ylabel("Relative error (log)")
ax2.set_title("float32 vs float64 round-trip")
ax2.legend(fontsize=8)
ax2.grid(True, which="both", alpha=0.3)
# Panel 3: hypot_safe magnitudes
ax3 = axes[2]
labels = [f"({a:.0e},{b:.0e})" for a, b in hypot_cases]
hypot_vals = [r["hypot"] for r in hypot_results]
finite_mask = np.isfinite(hypot_vals)
bar_vals = [v if np.isfinite(v) else 0.0 for v in hypot_vals]
colors_h = ["#b45309" if m else "#ef4444" for m in finite_mask]
ax3.bar(range(len(labels)), bar_vals, color=colors_h, alpha=0.85)
ax3.set_yscale("symlog", linthresh=1.0)
ax3.set_xticks(range(len(labels)))
ax3.set_xticklabels(labels, rotation=20, ha="right", fontsize=7)
ax3.set_ylabel("hypot_safe result (symlog)")
ax3.set_title("hypot_safe — extreme values")
ax3.grid(True, alpha=0.3, axis="y")
fig.suptitle(
"IMSL Numerical Utils: hypot_safe / prime_factors / cpu_time / float64 precision",
fontweight="bold",
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved to: {plot_path}")
return {
"machine_eps": machine_eps,
"float_max": float_max,
"float_min": float_min,
"hypot_results": hypot_results,
"factored": factored,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_numerical_utils()
Plot Output¶
Left: relative rounding error of float64 vs magnitude on a log–log scale.
Centre: float32 round-trip relative error showing ~7 decimal digits of
precision. Right: utilities.hypot_safe() results for extreme
(overflow/underflow) inputs on a symlog y-axis.¶
Console Output¶
IMSL Numerical Utilities Example
=======================================================
Machine Constants (float64 / IEEE 754 double)
Machine epsilon 2.220446e-16
Float max 1.797693e+308
Float min (tiny) 2.225074e-308
Mantissa bits 52
Decimal digits 15
Python sys.float_info.epsilon 2.220446e-16
hypot_safe ΓÇö overflow/underflow-safe hypotenuse
a b hypot_safe naive ok?
--------------------------------------------------------
3.000e+00 4.000e+00 5.000000e+00 yes
5.000e+00 1.200e+01 1.300000e+01 yes
1.000e+200 1.000e+200 1.414214e+200 overflow
1.000e-200 1.000e-200 1.414214e-200 yes
0.000e+00 7.300e+00 7.300000e+00 yes
prime_factors ΓÇö integer factorisation
n factors
----------------------------------------
12 2 × 2 × 3
60 2 × 2 × 3 × 5
100 2 × 2 × 5 × 5
360 2 × 2 × 2 × 3 × 3 × 5
1024 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2
9999 3 × 3 × 11 × 101
97 97
1000003 1000003
cpu_time (CPSEC): eig(200×200) took 15.625 ms CPU time
Plot saved to: test_output\example_imsl_numerical_utils.svg