271 lines
12 KiB
Python
271 lines
12 KiB
Python
"""Independent stdlib references for the five new exponential-family functions.
|
|
|
|
Run: python3 path/to/generate_exponential_functions.py [--check]
|
|
Exact binary64 high+low sums are formed at 2200 digits and checked with Fraction.
|
|
Decimal ln/exp (or a cancellation-safe expm1 series) runs at 450 and 650 digits;
|
|
all inputs, flags, and 120-significant-digit references must agree. Integer base-2
|
|
and base-10 powers use exact rational arithmetic, also checked against ln/exp.
|
|
No production arithmetic, binary64 transcendental functions, or rounded DD
|
|
constants are used. Zeros, nonfinite inputs, and enormous inputs belong in the
|
|
main suite. These rows specify an error bound, not correctly rounded DD results.
|
|
"""
|
|
from collections import Counter
|
|
from decimal import Decimal, localcontext
|
|
from fractions import Fraction
|
|
from pathlib import Path
|
|
import argparse
|
|
import math
|
|
import random
|
|
|
|
|
|
OUTPUT = Path(__file__).with_name('ExponentialFunctionsReferenceData.cs')
|
|
OPERATIONS = {'Exp2': 2, 'Exp10': 10, 'ExpM1': None, 'Exp2M1': 2, 'Exp10M1': 10}
|
|
EPSILON = math.ulp(0.0)
|
|
|
|
|
|
def exact_sum(high, low):
|
|
with localcontext() as context:
|
|
context.prec = 2200
|
|
value = Decimal.from_float(high) + Decimal.from_float(low)
|
|
assert Fraction(value) == Fraction(high) + Fraction(low), (high, low)
|
|
return value
|
|
|
|
|
|
def exact_decimal(value):
|
|
with localcontext() as context:
|
|
context.prec = 2200
|
|
result = Decimal(value.numerator) / Decimal(value.denominator)
|
|
assert Fraction(result) == value
|
|
return result
|
|
|
|
|
|
OVERFLOW = exact_decimal(Fraction(2) ** 1024 - Fraction(2) ** 970)
|
|
HALF_EPSILON = exact_decimal(Fraction(2) ** -1075)
|
|
|
|
|
|
def expm1(value):
|
|
# exp(value)-1 would erase tiny inputs even with hundreds of digits. For
|
|
# |value| <= 0.5, sum x + x^2/2! + ... until the context stops changing.
|
|
if abs(value) > Decimal('0.5'):
|
|
return value.exp() - 1
|
|
total = term = value
|
|
for denominator in range(2, 10000):
|
|
term = term * value / denominator
|
|
updated = total + term
|
|
if updated == total:
|
|
return updated
|
|
total = updated
|
|
raise AssertionError('expm1 series did not converge')
|
|
|
|
|
|
def neighbors(boundary):
|
|
high = float(boundary)
|
|
low = float(boundary - Decimal.from_float(high))
|
|
adjacent = [math.nextafter(low, -math.inf), low, math.nextafter(low, math.inf)]
|
|
assert exact_sum(high, adjacent[0]) < boundary < exact_sum(high, adjacent[-1])
|
|
assert math.nextafter(adjacent[0], math.inf) == adjacent[1]
|
|
assert math.nextafter(adjacent[1], math.inf) == adjacent[2]
|
|
return [(high, residual) for residual in adjacent]
|
|
|
|
|
|
def inputs(operation, log_base):
|
|
base = OPERATIONS[operation]
|
|
minus_one = operation.endswith('M1')
|
|
for high in [-10.25, -2.0, -1.0, -0.125, 0.125, 0.75, 1.0, 3.25]:
|
|
yield 'ordinary', high, 0.0
|
|
for high in [-1.0, 1.0]:
|
|
for low in [math.ulp(high) / 4, -math.ulp(high) / 4, EPSILON, -EPSILON]:
|
|
yield 'dense-and-sparse-low', high, low
|
|
for magnitude in [1e-20, 1e-100, 1e-300, math.ldexp(1.0, -1022), 2 * EPSILON, EPSILON]:
|
|
for sign in [-1, 1]:
|
|
yield 'tiny', sign * magnitude, 0.0
|
|
for sign in [-1, 1]:
|
|
for high, low in neighbors(Decimal(sign) / (2 * log_base)):
|
|
yield 'absolute-log-half-switch', high, low
|
|
if minus_one:
|
|
# Concrete dispatch uses the input high (0.5 for e/2, 0.125 for 10).
|
|
# Include both neighboring highs and lows on either side of the anchor.
|
|
switch = 0.125 if base == 10 else 0.5
|
|
for sign in [-1, 1]:
|
|
high = sign * switch
|
|
for adjacent in [math.nextafter(high, -math.inf), high, math.nextafter(high, math.inf)]:
|
|
yield 'input-high-series-switch', adjacent, 0.0
|
|
for low in [-math.ulp(high) / 4, math.ulp(high) / 4]:
|
|
yield 'input-high-series-switch', high, low
|
|
for high, low in neighbors(Decimal(sign) * Decimal(2) ** -54 / log_base):
|
|
yield 'quadratic-series-switch', high, low
|
|
if operation == 'Exp2':
|
|
for sign in [-1, 1]:
|
|
low = sign * math.ldexp(1.0, -500)
|
|
for adjacent in [math.nextafter(low, -math.inf), low, math.nextafter(low, math.inf)]:
|
|
yield 'sparse-correction-switch', 1000.0, adjacent
|
|
rng = random.Random(20260916 + (base or 1))
|
|
# Both local/moderate arguments and the complete useful exponential output
|
|
# range, expressed in natural-log coordinates, with reproducible dense lows.
|
|
for label, lower, upper in [('random-moderate', -8, 8), ('random-range', -748, 712)]:
|
|
for _ in range(6):
|
|
high = float(Decimal.from_float(rng.uniform(lower, upper)) / log_base)
|
|
yield label, high, rng.choice([-1, 1]) * math.ulp(high) / 4
|
|
if base == 2:
|
|
for exponent in [-1076, -1075, -1074, -1022, -100, -1, 1, 10, 53, 100, 1023, 1024]:
|
|
yield 'integer-power', float(exponent), 0.0
|
|
elif base == 10:
|
|
for exponent in [-324, -323, -308, -1, 1, 22, 308, 309]:
|
|
yield 'integer-power', float(exponent), 0.0
|
|
if minus_one:
|
|
for high in ([-40.0, -100.0, -745.0, -746.0] if base is None else
|
|
[-54.0, -100.0, -1075.0] if base == 2 else [-20.0, -100.0, -324.0]):
|
|
yield 'negative-saturation', high, 0.0
|
|
if base == 2:
|
|
# The exponential correction is exactly half epsilon at -1075,
|
|
# but exp2m1 is near -1, NOT an underflowing result.
|
|
for low in [-EPSILON, EPSILON]:
|
|
yield 'negative-saturation-sparse-low', -1075.0, low
|
|
# For M1, exp(x ln b) must reach OVERFLOW + 1, not OVERFLOW.
|
|
overflow_log = (OVERFLOW + int(minus_one)).ln() / log_base
|
|
for high, low in neighbors(overflow_log):
|
|
yield 'overflow-adjacent-low', high, low
|
|
if not minus_one:
|
|
# Base 2 has the exactly attainable threshold x=-1075; do not obtain
|
|
# it by an inexact quotient of two rounded logarithms.
|
|
underflow_log = Decimal(-1075) if base == 2 else HALF_EPSILON.ln() / log_base
|
|
for high, low in neighbors(underflow_log):
|
|
yield 'underflow-adjacent-low', high, low
|
|
normal_log = Decimal(-1022) if base == 2 else Decimal(2).ln() * -1022 / log_base
|
|
for high, low in neighbors(normal_log):
|
|
yield 'min-normal-adjacent-low', high, low
|
|
|
|
|
|
def evaluate(operation, high, low, log_base):
|
|
value = exact_sum(high, low)
|
|
base = OPERATIONS[operation]
|
|
minus_one = operation.endswith('M1')
|
|
argument = value * log_base
|
|
result = expm1(argument) if minus_one else argument.exp()
|
|
if base is not None and value == value.to_integral_value():
|
|
rational = Fraction(base) ** int(value) - int(minus_one)
|
|
exact = exact_decimal(rational)
|
|
# Cross-check integer fixtures by an independent mathematical identity.
|
|
assert abs(result - exact) <= abs(exact) * Decimal('1e-400')
|
|
result = exact
|
|
# Use exact input versus the independently computed thresholds: comparing a
|
|
# rounded exponential with half epsilon can misclassify an exact midpoint.
|
|
overflow = value >= (OVERFLOW + int(minus_one)).ln() / log_base
|
|
if minus_one:
|
|
underflow = abs(result) <= HALF_EPSILON
|
|
assert (result < 0) == (value < 0)
|
|
else:
|
|
threshold = Decimal(-1075) if base == 2 else HALF_EPSILON.ln() / log_base
|
|
underflow = value <= threshold
|
|
with localcontext() as output_context:
|
|
output_context.prec = 120
|
|
rounded = +result
|
|
assert abs(result - rounded) <= abs(result) * Decimal(2) ** -350
|
|
return format(rounded, 'e'), overflow, underflow
|
|
|
|
|
|
def generate_at_precision(precision):
|
|
rows = []
|
|
counts = Counter()
|
|
with localcontext() as context:
|
|
context.prec = precision
|
|
for operation, base in OPERATIONS.items():
|
|
log_base = Decimal(base).ln() if base else Decimal(1)
|
|
seen = set()
|
|
boundary_flags = {}
|
|
for label, high, low in inputs(operation, log_base):
|
|
assert math.isfinite(high) and math.isfinite(low)
|
|
assert high != 0 and abs(low) <= math.ulp(high) / 2
|
|
assert float(Fraction(high) + Fraction(low)) == high, (high, low)
|
|
result = evaluate(operation, high, low, log_base)
|
|
if label in ['overflow-adjacent-low', 'underflow-adjacent-low']:
|
|
flag_index = 1 if label.startswith('overflow') else 2
|
|
boundary_flags.setdefault(label, []).append(result[flag_index])
|
|
key = (high, low)
|
|
if key in seen:
|
|
continue
|
|
seen.add(key)
|
|
rows.append((operation, label, high, low, *result))
|
|
counts[operation] += 1
|
|
assert boundary_flags['overflow-adjacent-low'][0] is False
|
|
assert boundary_flags['overflow-adjacent-low'][-1] is True
|
|
if not operation.endswith('M1'):
|
|
assert boundary_flags['underflow-adjacent-low'][0] is True
|
|
assert boundary_flags['underflow-adjacent-low'][-1] is False
|
|
return rows, counts
|
|
|
|
|
|
def verify_sparse_splits(precision):
|
|
# Component-retention sentinels in DoubleDoubleExponentialFunctionsTests:
|
|
# unlike the 120-digit rows, these retain the tiny correction beside 2^n.
|
|
with localcontext() as context:
|
|
context.prec = precision
|
|
ln2 = Decimal(2).ln()
|
|
for exponent, expected_low in [(1, EPSILON), (500, 1.1210060331144859e-173),
|
|
(1000, 3.6694906201918696e-23)]:
|
|
for sign in [-1, 1]:
|
|
value = exact_sum(float(exponent), sign * EPSILON)
|
|
result = (value * ln2).exp()
|
|
high = math.ldexp(1.0, exponent)
|
|
assert float(result) == high
|
|
assert float(result - Decimal.from_float(high)) == sign * expected_low
|
|
|
|
|
|
def generate():
|
|
first, counts = generate_at_precision(450)
|
|
second, second_counts = generate_at_precision(650)
|
|
assert first == second and counts == second_counts, 'Precision stability check failed'
|
|
verify_sparse_splits(450)
|
|
verify_sparse_splits(650)
|
|
assert 200 <= len(first) <= 400
|
|
lines = []
|
|
previous_group = None
|
|
for operation, label, high, low, reference, overflow, underflow in first:
|
|
group = operation, label
|
|
if group != previous_group:
|
|
lines.append(f' // {operation}: {label}')
|
|
previous_group = group
|
|
lines.append(f' yield return new("{operation}", {high!r}, {low!r}, '
|
|
f'"{reference}", {str(overflow).lower()}, {str(underflow).lower()});')
|
|
content = '''// Generated by generate_exponential_functions.py; do not hand-edit.
|
|
// Exact binary64 sums verified with Fraction at 2200 digits; Decimal ln/exp and
|
|
// cancellation-safe expm1 at 450/650 digits; 120-digit references agree.
|
|
// Integer base-2/base-10 cases additionally use exact rational powers.
|
|
// Row: operation, high, low, reference, overflow, underflow.
|
|
// Overflow: exact result >= 2^1024 - 2^970 (binary64 overflow midpoint).
|
|
// Underflow: |exact result| <= 2^-1075, NOT exp(x) <= 2^-1075 for M1.
|
|
// Compare exact component sums with relative tolerance 2^-100 + double.Epsilon
|
|
// absolute; reference rounding uncertainty is bounded by 2^-350 relative.
|
|
// Tiny corrections to 1/-1 below 120 digits are not component-retention oracles.
|
|
// Signed zeros, infinities, NaNs, and enormous arguments are tested separately.
|
|
using Xunit;
|
|
|
|
namespace Just.PreciseMath.Tests.ReferenceData;
|
|
|
|
internal static class ExponentialFunctionsReferenceData
|
|
{
|
|
internal static IEnumerable<TheoryDataRow<string, double, double, string, bool, bool>> Cases()
|
|
{
|
|
''' + '\n'.join(lines) + '\n }\n}\n'
|
|
return content, counts
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument('--check', action='store_true', help='Verify fixtures without writing')
|
|
args = parser.parse_args()
|
|
content, counts = generate()
|
|
if args.check:
|
|
assert OUTPUT.read_text() == content, 'Fixture is stale; regenerate it'
|
|
action = 'Verified'
|
|
else:
|
|
OUTPUT.write_text(content)
|
|
action = 'Generated'
|
|
print(f'{action} {sum(counts.values())} references at 450/650 digits: {dict(counts)}')
|
|
print('Exact 2200-digit sums, rational integer powers, adjacent-low boundary flags, '
|
|
'normalized pairs, and reference rounding bounds verified.')
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|