logarithmic functions implementation

1-tests/Just.PreciseMath.Tests/ReferenceData/generate_logarithmic_functions.py

@@ -0,0 +1,137 @@
+"""Generate exact-input logarithmic-family references using only the standard library.
+
+Run this file with --check to verify the checked-in C# fixtures. Component sums
+and 1+x are formed exactly at 2200 digits and checked with Fraction. Decimal.ln
+and base conversion are evaluated independently at 450 and 650 digits, requiring
+identical 120-digit output and range flags. Tests do not require Python.
+"""
+from decimal import Decimal, localcontext
+from fractions import Fraction
+from pathlib import Path
+import argparse
+import math
+import random
+import sys
+
+from generate_log import exact_sum, inputs
+
+
+OUTPUT = Path(__file__).with_name('LogarithmicFunctionsReferenceData.cs')
+
+
+def plus_one_inputs():
+    epsilon = math.ulp(0.0)
+    yield -1.0, epsilon
+    yield -1.0, 2.0 ** -1000
+    yield -1.0, 2.0 ** -54
+    for high in [-0.999, -0.75, -0.5, -0.125, 0.125, 0.5, 1.0, 3.0, 9.0,
+                 sys.float_info.min, 1e100, sys.float_info.max]:
+        for low in [0.0, math.ulp(high) / 4, -math.ulp(high) / 4]:
+            yield high, low
+    for exponent in [-1074, -1073, -1022, -1000, -500, -100, -55, -54, -53]:
+        for sign in [-1, 1]:
+            high = sign * 2.0 ** exponent
+            for low in [0.0, math.ulp(high) / 4, -math.ulp(high) / 4]:
+                yield high, low
+    # Sparse lows and both sides of the actual high-component dispatch boundaries.
+    for high in [-2.0 ** -100, 2.0 ** -100]:
+        for low in [-2.0 ** -180, 2.0 ** -180]:
+            yield high, low
+    for high in [-0.5, 0.5, -2.0 ** -54, 2.0 ** -54]:
+        for neighbor in [math.nextafter(high, -math.inf), high, math.nextafter(high, math.inf)]:
+            for low in [-math.ulp(neighbor) / 4, 0.0, math.ulp(neighbor) / 4]:
+                yield neighbor, low
+    rng = random.Random(271828)
+    for _ in range(96):
+        high = rng.uniform(-0.99, 1.0)
+        yield high, rng.choice([-1, 1]) * math.ulp(high) / 4
+
+
+def base_inputs():
+    values = [(0.5, 0.0), (3.0, 2.0 ** -53), (10.0, -2.0 ** -51),
+              (math.ulp(0.0), 0.0), (sys.float_info.max, 0.0),
+              (1.0, 2.0 ** -54), (1.0, -2.0 ** -54),
+              (1.0, math.ulp(0.0)), (1.0, -math.ulp(0.0)),
+              (1.0, 2.0 ** -500), (1.0, -2.0 ** -500)]
+    for high, low in values:
+        for base_high, base_low in values:
+            yield high, low, base_high, base_low
+
+
+def cases():
+    for _, high, low in inputs():
+        for operation in ['Log2', 'Log10']:
+            yield operation, high, low, 0.0, 0.0
+    for high, low in plus_one_inputs():
+        for operation in ['LogP1', 'Log2P1', 'Log10P1']:
+            yield operation, high, low, 0.0, 0.0
+    for high, low, base_high, base_low in base_inputs():
+        yield 'LogBase', high, low, base_high, base_low
+
+
+def reference(operation, high, low, base_high, base_low, precision):
+    value = exact_sum(high, low)
+    if operation.endswith('P1'):
+        with localcontext() as context:
+            context.prec = 2200
+            value += 1
+        assert Fraction(value) == Fraction(high) + Fraction(low) + 1
+    assert value > 0
+    base = exact_sum(base_high, base_low) if operation == 'LogBase' else None
+    with localcontext() as context:
+        context.prec = precision
+        result = value.ln()
+        if base is not None:
+            result /= base.ln()
+        elif operation.startswith('Log2'):
+            result /= Decimal(2).ln()
+        elif operation.startswith('Log10'):
+            result /= Decimal(10).ln()
+        with localcontext() as output_context:
+            output_context.prec = 120
+            rounded = +result
+        assert abs(result - rounded) <= abs(rounded) * Decimal(2) ** -350
+        # Independently classify against binary64's exact range midpoints.
+        magnitude = abs(Fraction(result))
+        overflow = magnitude >= Fraction((2 ** 54 - 1) * 2 ** 970)
+        underflow = magnitude != 0 and magnitude <= Fraction(1, 2 ** 1075)
+        return format(rounded, 'e') if rounded else '0e+0', overflow, underflow
+
+
+def generate():
+    rows = []
+    for case in dict.fromkeys(cases()):
+        result = reference(*case, 450)
+        assert result == reference(*case, 650), case
+        operation, high, low, base_high, base_low = case
+        text, overflow, underflow = result
+        rows.append(f'        yield return new("{operation}", {high!r}, {low!r}, '
+                    f'{base_high!r}, {base_low!r}, "{text}", '
+                    f'{str(overflow).lower()}, {str(underflow).lower()});')
+    content = '''// Generated by generate_logarithmic_functions.py; do not hand-edit reference literals.
+// Exact component sums and 1+x; Decimal.ln at 450/650 digits, rounded to 120 digits.
+namespace Just.PreciseMath.Tests.ReferenceData;
+
+internal static class LogarithmicFunctionsReferenceData
+{
+    internal static IEnumerable<TheoryDataRow<string, double, double, double, double, string, bool, bool>> 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')
+    args = parser.parse_args()
+    content, count = generate()
+    if args.check:
+        assert OUTPUT.read_text() == content, 'Fixture is stale; regenerate it'
+        print(f'Verified {count} logarithmic-family references at 450/650 digits.')
+    else:
+        OUTPUT.write_text(content)
+        print(f'Generated {count} references in {OUTPUT.name}.')
+
+
+if __name__ == '__main__':
+    main()