"""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()