152 lines
6.9 KiB
Python
152 lines
6.9 KiB
Python
"""Independent finite Pow fixtures; Python standard library only.
|
|
|
|
Run from any directory: python3 path/to/generate_real_pow.py [--check].
|
|
Build each binary64 component sum EXACTLY at 2200 decimal digits and verify
|
|
against Fraction. Decimal.ln/exp evaluate at 450 and 650 digits; the rounded
|
|
120-significant-digit references must agree. Even 1 +/- 2^-1074 is retained
|
|
before logarithm evaluation, so huge exponents can amplify sparse corrections.
|
|
No library-under-test operations or binary64 Pow supply expected values.
|
|
Each scalar case discards exponent low BEFORE independently evaluating its oracle.
|
|
"""
|
|
from decimal import Decimal, localcontext
|
|
from fractions import Fraction
|
|
from pathlib import Path
|
|
import argparse
|
|
import math
|
|
import random
|
|
import sys
|
|
|
|
|
|
OUTPUT = Path(__file__).with_name('RealPowReferenceData.cs')
|
|
|
|
|
|
def exact_sum(high, low):
|
|
with localcontext() as context:
|
|
context.prec = 2200
|
|
result = Decimal.from_float(high) + Decimal.from_float(low)
|
|
assert Fraction(result) == Fraction(high) + Fraction(low), (high, low)
|
|
return result
|
|
|
|
|
|
def inputs():
|
|
# (label, base high, base low, exponent high, exponent low).
|
|
rows = []
|
|
for high in [0.125, 0.5, 1.25, 2.0, 10.0, 1e-308, 1e308,
|
|
sys.float_info.max, sys.float_info.min, math.ulp(0.0),
|
|
3 * math.ulp(0.0), math.nextafter(sys.float_info.min, 0.0)]:
|
|
for power in [0.5, -0.5, 0.25, -0.25]:
|
|
rows.append(('magnitude-fractional', high, 0.0, power, 0.0))
|
|
for high in [0.75, 1.25, 2.0, 10.0, 1e-300, 1e300]:
|
|
for sign in [-1, 1]:
|
|
rows.append(('both-residuals', high, sign * math.ulp(high) / 4,
|
|
-0.75, -sign * math.ulp(0.75) / 4))
|
|
for low in [2.0 ** -54, -2.0 ** -54, 1e-300, -1e-300, 1e-308, -1e-308,
|
|
math.ulp(0.0), -math.ulp(0.0)]:
|
|
powers = [1e16, -1e16] if abs(low) > 1e-100 else [1e300, -1e300]
|
|
if abs(low) < 1e-307:
|
|
powers.extend([sys.float_info.max, -sys.float_info.max])
|
|
for power in powers:
|
|
for exponent_low in [0.0, math.ulp(power) / 4]:
|
|
rows.append(('near-one-amplification', 1.0, low, power, exponent_low))
|
|
# Tiny deltas on both sides of likely sparse-kernel cutoffs.
|
|
for scale in [-26, -27, -40, -53, -100, -500, -1000, -1022, -1073]:
|
|
for sign in [-1, 1]:
|
|
rows.append(('sparse-scale-sweep', 1.0, sign * 2.0 ** scale,
|
|
math.ldexp(0.75, min(-scale, 1023)), 0.125))
|
|
# Nonzero exponent lows distinguish DD from scalar, including high integers.
|
|
for high in [2.0, 10.0, 1e308, 1e-308]:
|
|
for power, low in [(1.0, 2.0 ** -54), (-1.0, -2.0 ** -54),
|
|
(0.0, math.ulp(0.0)), (0.5, 2.0 ** -55)]:
|
|
rows.append(('exponent-low-decisive', high, 0.0, power, low))
|
|
for power in [float(-(2 ** 31)), float(2 ** 31 - 1), float(2 ** 31),
|
|
float(-(2 ** 31) - 1), float(2 ** 53), float(2 ** 54)]:
|
|
for low in [-0.5, 0.0, 0.5]:
|
|
rows.append(('integer-boundary', 1.0, 2.0 ** -54, power, low))
|
|
# Exact DD integer parity includes low; scalar has its own (usually even) sign.
|
|
for power, low in [(2.0 ** 53, 1.0), (2.0 ** 54, -1.0),
|
|
(-2.0 ** 53, -1.0), (2.0 ** 54, 2.0)]:
|
|
rows.append(('negative-integral-low-parity', -1.0, -2.0 ** -54, power, low))
|
|
for power in [-63.0, -7.0, 7.0, 63.0]:
|
|
rows.append(('negative-int-dispatch', -1.25, 2.0 ** -55, power, 0.0))
|
|
# Finite range shoulders are far from ambiguous overflow/zero thresholds.
|
|
for power in [1023.75, -1021.75, -1073.5, -1074.25]:
|
|
rows.append(('finite-range-shoulder', 2.0, 2.0 ** -54, power, 2.0 ** -45))
|
|
rng = random.Random(1618033)
|
|
for index in range(96):
|
|
high = math.ldexp(rng.uniform(1.125, 1.875), rng.randint(-1000, 1000))
|
|
low = rng.choice([-1, 1]) * math.ulp(high) / 4
|
|
# Binary64 log selects well-separated finite cases, NEVER their references.
|
|
target_log = rng.uniform(-740.0, 708.0)
|
|
power = target_log / math.log(high)
|
|
rows.append((f'random-{index:03}', high, low, power,
|
|
rng.choice([-1, 1]) * math.ulp(power) / 4))
|
|
return rows
|
|
|
|
|
|
def reference(high, low, power_high, power_low, precision):
|
|
base = exact_sum(high, low)
|
|
power = exact_sum(power_high, power_low)
|
|
sign = 1
|
|
if base < 0:
|
|
assert power == power.to_integral_value(), (base, power)
|
|
sign = -1 if int(power) % 2 else 1
|
|
base = base.copy_abs()
|
|
with localcontext() as context:
|
|
context.prec = precision
|
|
result = (base.ln() * power).exp() * sign
|
|
# No disputed overflow boundary expectations: leave a large margin.
|
|
assert result.copy_abs() < Decimal.from_float(sys.float_info.max) * Decimal('0.999999999999')
|
|
with localcontext() as output_context:
|
|
output_context.prec = 120
|
|
rounded = +result
|
|
# Verify a generous 2^-350 reference-relative allowance explicitly.
|
|
assert abs(result - rounded) <= abs(rounded) * Decimal(2) ** -350
|
|
return format(rounded, 'e')
|
|
|
|
|
|
def generate():
|
|
rows = []
|
|
references = 0
|
|
for label, high, low, power, power_low in inputs():
|
|
for scalar in [False, True]:
|
|
oracle_low = 0.0 if scalar else power_low
|
|
first = reference(high, low, power, oracle_low, 450)
|
|
assert first == reference(high, low, power, oracle_low, 650), (label, scalar)
|
|
references += 1
|
|
rows.append(f' yield return new("{label}", {high!r}, {low!r}, '
|
|
f'{power!r}, {power_low!r}, {str(scalar).lower()}, "{first}");')
|
|
content = '''// Generated by generate_real_pow.py; do not hand-edit reference literals.
|
|
// Exact binary64 sums; Decimal ln/exp at 450/650 digits, rounded to 120 digits.
|
|
// Reference uncertainty is explicitly allowed as 2^-350 relative in the tests.
|
|
using Xunit;
|
|
|
|
namespace Just.PreciseMath.Tests.ReferenceData;
|
|
|
|
internal static class RealPowReferenceData
|
|
{
|
|
public static IEnumerable<TheoryDataRow<string, double, double, double, double, bool, string>> Cases()
|
|
{
|
|
''' + '\n'.join(rows) + '\n }\n}\n'
|
|
return content, references
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument('--check', action='store_true', help='Verify checked-in fixture without writing it')
|
|
args = parser.parse_args()
|
|
# Sparse input construction sentinels include more than 1000 decimal places.
|
|
for high, low in [(1.0, math.ulp(0.0)), (1.0, -math.ulp(0.0)),
|
|
(sys.float_info.max, math.ulp(0.0))]:
|
|
exact_sum(high, low)
|
|
content, count = generate()
|
|
if args.check:
|
|
assert OUTPUT.read_text() == content, 'Fixture is stale; regenerate it'
|
|
print(f'Verified {count} finite references at 450/650 digits and exact component sums.')
|
|
else:
|
|
OUTPUT.write_text(content)
|
|
print(f'Generated {count} finite references in {OUTPUT.name}.')
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|