111 lines
4.5 KiB
Python
111 lines
4.5 KiB
Python
"""Generate Log fixtures independently with Python's standard library.
|
|
|
|
Run from any directory: python3 path/to/generate_log.py [--check].
|
|
Construct exact binary64 component sums at 2200 decimal digits, verify against
|
|
Fraction, and require identical 120-digit Decimal.ln references at 450/650 digits.
|
|
The C# suite also multiplies the independent ln(2) reference by integer exponents
|
|
using BigInteger to check powers of two across every binary64 input exponent.
|
|
"""
|
|
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('LogReferenceData.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 reference(high, low, precision):
|
|
value = exact_sum(high, low)
|
|
assert value > 0
|
|
with localcontext() as context:
|
|
context.prec = precision
|
|
result = value.ln()
|
|
with localcontext() as output_context:
|
|
output_context.prec = 120
|
|
rounded = +result
|
|
assert abs(result - rounded) <= abs(rounded) * Decimal(2) ** -350
|
|
return format(rounded, 'e')
|
|
|
|
|
|
def inputs():
|
|
for high in [0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 3.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),
|
|
math.nextafter(sys.float_info.min, math.inf),
|
|
math.nextafter(1.0, 0.0), math.nextafter(1.0, math.inf)]:
|
|
yield 'scalar-boundary', high, 0.0
|
|
for scale in [-1074, -1073, -1022, -1000, -500, -1, 0, 1, 500, 1000, 1023]:
|
|
high = math.ldexp(1.5, scale)
|
|
for low in [math.ulp(high) / 4, -math.ulp(high) / 4]:
|
|
yield 'binade-residual', high, low
|
|
for high in [1.0, math.sqrt(2), 1 / math.sqrt(2), 1e-308, 1e308, sys.float_info.max]:
|
|
for low in [math.ulp(high) / 4, -math.ulp(high) / 4, math.ulp(0.0), -math.ulp(0.0)]:
|
|
yield 'dense-and-sparse-low', high, low
|
|
# Same high with adjacent lows crosses the tiny-delta kernel cutoff.
|
|
for low in [2.0 ** -54, -2.0 ** -54]:
|
|
for neighbor in [math.nextafter(low, -math.inf), low, math.nextafter(low, math.inf)]:
|
|
yield 'tiny-delta-transition', 1.0, neighbor
|
|
for scale in [-53, -54, -55, -100, -500, -1000, -1022, -1073, -1074]:
|
|
for sign in [-1, 1]:
|
|
yield 'near-one-log', 1.0, sign * 2.0 ** scale
|
|
# Centering transitions in both high and low, far from x=1 cancellation.
|
|
for high in [math.sqrt(2), 1 / math.sqrt(2)]:
|
|
for neighbor in [math.nextafter(high, 0.0), high, math.nextafter(high, math.inf)]:
|
|
for low in [math.ulp(neighbor) / 4, -math.ulp(neighbor) / 4]:
|
|
yield 'centering-transition', neighbor, low
|
|
rng = random.Random(3141592)
|
|
for index in range(256):
|
|
high = math.ldexp(rng.uniform(1.0, 1.999), rng.randint(-1074, 1023))
|
|
yield f'random-{index:03}', high, rng.choice([-1, 0, 1]) * math.ulp(high) / 4
|
|
|
|
|
|
def generate():
|
|
ln2 = reference(2.0, 0.0, 450)
|
|
assert ln2 == reference(2.0, 0.0, 650)
|
|
rows = []
|
|
for label, high, low in inputs():
|
|
first = reference(high, low, 450)
|
|
assert first == reference(high, low, 650), (label, high, low)
|
|
rows.append(f' yield return new("{label}", {high!r}, {low!r}, "{first}");')
|
|
content = '''// Generated by generate_log.py; do not hand-edit reference literals.
|
|
// Exact binary64 sums; Decimal.ln at 450/650 digits, rounded to 120 digits.
|
|
using Xunit;
|
|
|
|
namespace Just.PreciseMath.Tests.ReferenceData;
|
|
|
|
internal static class LogReferenceData
|
|
{
|
|
''' + f' internal const string Ln2 = "{ln2}";\n\n' + ''' internal static IEnumerable<TheoryDataRow<string, double, double, string>> Cases()
|
|
{
|
|
''' + '\n'.join(rows) + '\n }\n}\n'
|
|
return content, len(rows)
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument('--check', action='store_true', help='Verify fixtures without writing')
|
|
args = parser.parse_args()
|
|
content, count = generate()
|
|
if args.check:
|
|
assert OUTPUT.read_text() == content, 'Fixture is stale; regenerate it'
|
|
print(f'Verified {count} Log references and ln(2) at 450/650 digits with exact input sums.')
|
|
else:
|
|
OUTPUT.write_text(content)
|
|
print(f'Generated {count} Log references in {OUTPUT.name}.')
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|