#!/usr/bin/env python3
"""Project Euler Problem 1007 - Alternating Difference."""

from __future__ import annotations

MOD = 1_000_000_009
TARGET = 10_000_000


def normalize(value: int) -> int:
    return value % MOD


def alternating_sum(n: int) -> int:
    """Return A(n) modulo MOD using the quotient ring r^2 = r + 1."""

    inverse_four = pow(4, MOD - 2, MOD)

    # Every ring value is stored as two scalars: constant + root * r.
    a_scaled = 1
    b_constant, b_root = 1, 0
    c_older_constant, c_older_root = 0, 0
    c_previous_constant, c_previous_root = 1, 0
    result_constant, result_root = 0, 0
    factorial = 1

    for m in range(1, n + 2):
        linear_factor = (4 * m - 6) % MOD
        m_mod = m % MOD

        a_next = linear_factor * a_scaled % MOD

        # r(a + br) = b + (a + b)r.
        b_next_constant = linear_factor * b_root % MOD
        b_next_root = linear_factor * (b_constant + b_root) % MOD

        # (r + 1)(a + br) = (a + b) + (a + 2b)r.
        c_next_constant = linear_factor * (
            c_previous_constant + c_previous_root
        ) % MOD
        c_next_root = linear_factor * (
            c_previous_constant + 2 * c_previous_root
        ) % MOD

        quadratic_factor = (
            16 * ((m - 3) % MOD) % MOD * ((m - 1) % MOD) % MOD
        )
        # Subtract quadratic_factor * r * c_older.
        c_next_constant = (
            c_next_constant - quadratic_factor * c_older_root
        ) % MOD
        c_next_root = (
            c_next_root
            - quadratic_factor * (c_older_constant + c_older_root)
        ) % MOD

        # Build n_m from q(z), q(rz), and q(z)q(rz).  The identities
        # (r - 1)(a + br) = (b - a) + ar are expanded componentwise.
        numerator_constant = (
            -a_next
            - m_mod * a_scaled
            + b_next_constant
            + m_mod * (b_root - b_constant)
            - c_next_constant
        ) % MOD
        numerator_root = (
            m_mod * a_scaled
            + b_next_root
            + m_mod * b_constant
            - c_next_root
        ) % MOD
        if m == 1:
            numerator_constant = (numerator_constant + 2) % MOD
            numerator_root = (numerator_root + 2) % MOD

        # Subtract 2m(2-r)u_{m-1}, where
        # (2-r)(a + br) = (2a-b) + (b-a)r.
        twice_m = 2 * m_mod % MOD
        adjusted_constant = (
            numerator_constant
            - twice_m * (2 * result_constant - result_root)
        ) % MOD
        adjusted_root = (
            numerator_root
            - twice_m * (result_root - result_constant)
        ) % MOD

        # u_m = (2-r) * adjusted / 4.
        result_next_constant = (
            (2 * adjusted_constant - adjusted_root) * inverse_four % MOD
        )
        result_next_root = (
            (adjusted_root - adjusted_constant) * inverse_four % MOD
        )

        a_scaled = a_next
        b_constant, b_root = b_next_constant, b_next_root
        c_older_constant, c_older_root = (
            c_previous_constant,
            c_previous_root,
        )
        c_previous_constant, c_previous_root = c_next_constant, c_next_root
        result_constant, result_root = result_next_constant, result_next_root
        factorial = factorial * m_mod % MOD

    return result_root * pow(factorial, MOD - 2, MOD) % MOD


def brute_values(fibonacci: list[int], begin: int, end: int) -> list[int]:
    if end - begin == 1:
        return [fibonacci[begin]]

    values: list[int] = []
    for split in range(begin + 1, end):
        left = brute_values(fibonacci, begin, split)
        right = brute_values(fibonacci, split, end)
        for lhs in left:
            for rhs in right:
                values.append(lhs - rhs)
    return values


def brute_alternating_sum(n: int) -> int:
    fibonacci = [0] * (n + 1)
    if n >= 1:
        fibonacci[1] = 1
    for index in range(2, n + 1):
        fibonacci[index] = fibonacci[index - 1] + fibonacci[index - 2]

    return sum(brute_values(fibonacci, 0, n + 1)) % MOD


def run_checkpoints() -> None:
    for n in range(11):
        assert alternating_sum(n) == brute_alternating_sum(n), (
            f"brute force comparison for n={n}"
        )

    assert alternating_sum(3) == normalize(-6), "published A(3)"
    assert alternating_sum(10) == normalize(-177_666), "published A(10)"
    assert alternating_sum(100) == 71_792_794, "published A(100)"


def main() -> None:
    run_checkpoints()
    print(alternating_sum(TARGET))


if __name__ == "__main__":
    main()