"""Trusted reference oracle for Bernoulli numbers.

This module is the **test oracle** for the verification lab. Everything the
learner generates is judged against it, so it is deliberately small, exact, and
independently checkable.

Conventions fixed for this module
---------------------------------
1. **First Bernoulli numbers** (also called ``B_n^-``, the convention of NIST's
   DLMF section 24.2): ``B1 = -1/2``.  The competing *second* convention has
   ``B1 = +1/2``; every other value is identical.  This is the single most common
   source of disagreement between two correct-looking programs, so the course
   teaches it explicitly rather than hiding it.
2. **Modern indexing**: ``bernoulli(n)`` returns ``B_n`` for ``n >= 0``, covering
   every index, including the odd ones that are zero.  Ada Lovelace's Note G uses
   the 19th-century convention in which only the non-zero terms are numbered, so
   her ``B1, B3, B5, B7`` are modern ``B2, B4, B6, B8``.  Use
   :func:`ada_to_modern_index` to move between the two; never mix them silently.
3. **Exact arithmetic**: all values are :class:`fractions.Fraction`.  Bernoulli
   numbers have large numerators and denominators and are not representable in
   binary floating point, so floats are treated as an error, not a rounding
   detail.

Algorithm
---------
The recurrence used here follows directly from the definition and needs no
special functions::

    sum_{j=0}^{n} C(n+1, j) * B_j = 0   for n >= 1

Rearranged for the highest term::

    B_n = -1/(n+1) * sum_{j=0}^{n-1} C(n+1, j) * B_j

This is the "second-convention" recurrence's sibling; with ``B0 = 1`` it produces
``B1 = -1/2`` directly, which is the convention this module fixes.

Independent verification
------------------------
:data:`PUBLISHED_TABLE` holds values transcribed from published tables so the
oracle can be checked against a source it did not compute.  ``B30`` is included
deliberately: small indices hide recurrence bugs that only appear once the
denominators grow.

Run this file directly to print the table the course quotes::

    python3 reference_bernoulli.py
"""

from __future__ import annotations

from fractions import Fraction
from functools import lru_cache
from math import comb

__all__ = [
    "CONVENTION",
    "PUBLISHED_TABLE",
    "ADA_NOTE_G_INDEX_MAP",
    "bernoulli",
    "bernoulli_sequence",
    "ada_bernoulli",
    "ada_to_modern_index",
    "modern_to_ada_index",
    "verify_against_published_table",
]

#: Human-readable name of the convention this oracle implements.
CONVENTION = "first Bernoulli numbers (B_n^-), B1 = -1/2"


@lru_cache(maxsize=None)
def bernoulli(n: int) -> Fraction:
    """Return ``B_n`` as an exact :class:`~fractions.Fraction`.

    This is the interface every submitted implementation must match:
    ``bernoulli(n: int) -> fractions.Fraction`` defined for ``n >= 0``, using
    ``B1 = -1/2``.

    >>> bernoulli(0)
    Fraction(1, 1)
    >>> bernoulli(1)
    Fraction(-1, 2)
    >>> bernoulli(4)
    Fraction(-1, 30)
    >>> bernoulli(3)
    Fraction(0, 1)
    """
    if not isinstance(n, int) or isinstance(n, bool):
        raise TypeError(f"n must be an int, got {type(n).__name__}")
    if n < 0:
        raise ValueError(f"n must be >= 0, got {n}")

    if n == 0:
        return Fraction(1)
    # Every odd index above 1 is zero. Returning early keeps the recurrence
    # cheap and states the property in code where a reader can see it.
    if n > 1 and n % 2 == 1:
        return Fraction(0)

    total = Fraction(0)
    for j in range(n):
        total += comb(n + 1, j) * bernoulli(j)
    return -total / (n + 1)


def bernoulli_sequence(start: int = 0, stop: int = 13) -> list[tuple[int, Fraction]]:
    """Return ``[(n, B_n), ...]`` for ``start <= n < stop``.

    Provided so learners can explore a range without writing a loop.
    """
    if start < 0:
        raise ValueError(f"start must be >= 0, got {start}")
    if stop < start:
        raise ValueError(f"stop ({stop}) must be >= start ({start})")
    return [(n, bernoulli(n)) for n in range(start, stop)]


def ada_to_modern_index(ada_index: int) -> int:
    """Map a Note G index to the modern index.

    Lovelace numbers only the non-zero Bernoulli numbers, so her ``B1`` is the
    modern ``B2``.  The mapping is ``modern = ada + 1`` for odd Ada indices.

    >>> ada_to_modern_index(1)
    2
    >>> ada_to_modern_index(7)
    8
    """
    if ada_index < 1 or ada_index % 2 == 0:
        raise ValueError(
            f"Note G numbers only odd indices B1, B3, B5, ...; got {ada_index}"
        )
    return ada_index + 1


def modern_to_ada_index(modern_index: int) -> int:
    """Map a modern even index to Lovelace's Note G index.

    >>> modern_to_ada_index(8)
    7
    """
    if modern_index < 2 or modern_index % 2 == 1:
        raise ValueError(
            f"only even modern indices >= 2 have a Note G label; got {modern_index}"
        )
    return modern_index - 1


def ada_bernoulli(ada_index: int) -> Fraction:
    """Return the value Lovelace would have called ``B<ada_index>``.

    >>> ada_bernoulli(1)
    Fraction(1, 6)
    >>> ada_bernoulli(7)
    Fraction(-1, 30)
    """
    return bernoulli(ada_to_modern_index(ada_index))


#: Exact values checked against NIST DLMF Tables 24.2.1 and 24.2.3, whose source
#: note identifies Abramowitz and Stegun (1964), pp. 809-810. Entries omitted
#: here are odd indices, which DLMF equation 24.2.2 states are zero for n > 1.
#: Full bibliographic and comparison details are in course/references.md [M2].
PUBLISHED_TABLE: dict[int, Fraction] = {
    0: Fraction(1, 1),
    1: Fraction(-1, 2),
    2: Fraction(1, 6),
    3: Fraction(0, 1),
    4: Fraction(-1, 30),
    5: Fraction(0, 1),
    6: Fraction(1, 42),
    7: Fraction(0, 1),
    8: Fraction(-1, 30),
    9: Fraction(0, 1),
    10: Fraction(5, 66),
    11: Fraction(0, 1),
    12: Fraction(-691, 2730),
    14: Fraction(7, 6),
    16: Fraction(-3617, 510),
    18: Fraction(43867, 798),
    20: Fraction(-174611, 330),
    22: Fraction(854513, 138),
    24: Fraction(-236364091, 2730),
    26: Fraction(8553103, 6),
    28: Fraction(-23749461029, 870),
    # Included deliberately: a large index where a subtly wrong recurrence will
    # diverge even though B0..B12 looked right.
    30: Fraction(8615841276005, 14322),
}

#: Lovelace's Note G labels against modern indices and values. Note G's worked
#: example computes her ``B7``, i.e. the modern ``B8``.
ADA_NOTE_G_INDEX_MAP: dict[int, int] = {1: 2, 3: 4, 5: 6, 7: 8}


def _recurrence_only(stop: int) -> list[Fraction]:
    """Return ``[B_0, ..., B_(stop-1)]`` from the recurrence alone.

    :func:`bernoulli` returns zero for odd ``n > 1`` without computing it. That
    shortcut states the property, but a check that reads the shortcut back can
    never fail. This computes every index, odd ones included, through the
    recurrence, so the self-check below tests the property instead of
    repeating it.
    """
    values: list[Fraction] = []
    for n in range(stop):
        if n == 0:
            values.append(Fraction(1))
            continue
        total = sum((comb(n + 1, j) * values[j] for j in range(n)), Fraction(0))
        values.append(-total / (n + 1))
    return values


def verify_against_published_table() -> list[str]:
    """Check the recurrence against :data:`PUBLISHED_TABLE`.

    Every index is computed through the recurrence with no shortcut for the
    odd terms, and that is what is compared: with every published value, with
    zero at every odd index from 3 to 31, and with what :func:`bernoulli`
    returns at every index from 0 to 31.

    Returns a list of human-readable problems; an empty list means the oracle
    agrees with every published value.
    """
    problems: list[str] = []
    computed = _recurrence_only(32)
    for n, expected in sorted(PUBLISHED_TABLE.items()):
        got = computed[n]
        if got != expected:
            problems.append(f"B{n}: oracle gave {got}, published table gives {expected}")
        if not isinstance(bernoulli(n), Fraction):
            problems.append(
                f"B{n}: oracle returned {type(bernoulli(n)).__name__}, not Fraction"
            )

    # Property check: every odd index above 1 must come out of the recurrence
    # as exactly zero.
    for n in range(3, 32, 2):
        if computed[n] != 0:
            problems.append(f"B{n}: odd index should be 0, recurrence gave {computed[n]}")

    # The values the oracle serves must be the values the recurrence computes,
    # odd-index shortcut included.
    for n, value in enumerate(computed):
        if bernoulli(n) != value:
            problems.append(
                f"B{n}: oracle returns {bernoulli(n)}, recurrence gives {value}"
            )

    return problems


def _main() -> int:
    problems = verify_against_published_table()

    print(f"Bernoulli oracle - convention: {CONVENTION}")
    print()
    print("  n   B_n")
    print("  --  ---------------------------")
    for n, value in bernoulli_sequence(0, 13):
        print(f"  {n:>2}  {value}")
    print()
    print("Larger index, beyond what small-index checks reach:")
    print(f"  B30 = {bernoulli(30)}")
    print()
    print("Note G indexing (Lovelace's label -> modern label -> value):")
    for ada_index, modern_index in sorted(ADA_NOTE_G_INDEX_MAP.items()):
        print(
            f"  Ada B{ada_index}  ->  modern B{modern_index}  =  {bernoulli(modern_index)}"
        )
    print()

    if problems:
        print("FAIL: oracle disagrees with the published table:")
        for problem in problems:
            print(f"  - {problem}")
        return 1

    checked = len(PUBLISHED_TABLE)
    print(f"OK: oracle agrees with all {checked} published values, odd terms are zero.")
    return 0


if __name__ == "__main__":
    raise SystemExit(_main())