Problem 1007: Alternating Difference

Project Euler Problem 1007 Solution

EulerSolve provides an optimized solution for Project Euler Problem 1007, Alternating Difference, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Place $F_0,F_1,\ldots,F_n$ in order with a minus sign between adjacent terms. A valid expression is a full parenthesization: every one of the $n$ pairs of parentheses encloses exactly one top-level minus sign, although either side may contain further parenthesized subexpressions. The task is to add the values of all syntactically different expressions. There are $C_n$ such expressions, where $C_n$ is the $n$-th Catalan number. Direct construction is therefore impossible at $n=10^7$. We need $$A(10^7)\pmod{p},\qquad p=10^9+9,$$ with the published checks $A(3)=-6$, $A(10)=-177666$, and $A(100)\equiv71792794\pmod p$. Mathematical Approach Parenthesizations are ordered full binary trees A parenthesized subtraction expression with $m=n+1$ operands is an ordered full binary tree with $m$ leaves. Its root chooses a unique split after $k$ leaves: the left subtree uses the first $k$ operands and the right subtree uses the remaining $m-k$. A block of $j+1$ leaves has $$C_j=\frac{1}{j+1}\binom{2j}{j}$$ possible trees. The root split is the key because the value of the whole tree is “left value minus right value,” and every valid tree occurs at exactly one split. A geometric probe records every leaf coefficient at once Temporarily replace the Fibonacci operands by $1,t,t^2,\ldots,t^{m-1}$....

Detailed mathematical approach

Problem Summary

Place $F_0,F_1,\ldots,F_n$ in order with a minus sign between adjacent terms. A valid expression is a full parenthesization: every one of the $n$ pairs of parentheses encloses exactly one top-level minus sign, although either side may contain further parenthesized subexpressions. The task is to add the values of all syntactically different expressions.

There are $C_n$ such expressions, where $C_n$ is the $n$-th Catalan number. Direct construction is therefore impossible at $n=10^7$. We need

$$A(10^7)\pmod{p},\qquad p=10^9+9,$$

with the published checks $A(3)=-6$, $A(10)=-177666$, and $A(100)\equiv71792794\pmod p$.

Mathematical Approach

Parenthesizations are ordered full binary trees

A parenthesized subtraction expression with $m=n+1$ operands is an ordered full binary tree with $m$ leaves. Its root chooses a unique split after $k$ leaves: the left subtree uses the first $k$ operands and the right subtree uses the remaining $m-k$. A block of $j+1$ leaves has

$$C_j=\frac{1}{j+1}\binom{2j}{j}$$

possible trees. The root split is the key because the value of the whole tree is “left value minus right value,” and every valid tree occurs at exactly one split.

A geometric probe records every leaf coefficient at once

Temporarily replace the Fibonacci operands by $1,t,t^2,\ldots,t^{m-1}$. Let $P_m(t)$ be the sum of the values of all parenthesizations of these $m$ formal operands. In particular, $P_1(t)=1$.

For a split after $k$ leaves, the left sum $P_k(t)$ is repeated once for every right tree, while the right sum is repeated once for every left tree and is shifted by $t^k$. Therefore

$$P_m(t)=\sum_{k=1}^{m-1} \left(C_{m-k-1}P_k(t)-C_{k-1}t^kP_{m-k}(t)\right).$$

This recurrence does not merely count trees. It preserves the complete signed contribution of every leaf position. For example,

$$P_1(t)=1,\qquad P_2(t)=1-t,\qquad P_3(t)=2-2t.$$

Worked sign-polynomial example: $n=3$

For four leaves, adding the five tree polynomials gives

$$P_4(t)=5-5t+t^2-t^3.$$

In the Fibonacci ring, $r^2=1+r$ and $r^3=1+2r$. Hence the root coefficient is

$$[r]P_4(r)=-5+1-2=-6,$$

which is the sum $-4-2+0+2-2=-6$ from the five displayed expressions in the problem. This example simultaneously checks the Catalan multiplicities, the right-subtree sign, the exponent shift, and the Fibonacci extraction rule.

Collapse the Catalan convolution with generating functions

Introduce

$$T(z)=\sum_{j\ge0}C_jz^j=\frac{1-\sqrt{1-4z}}{2z}, \qquad P(z,t)=\sum_{m\ge1}P_m(t)z^m.$$

Thus $T(z)$ is simply the ordinary Catalan generating function (also often denoted $C(z)$), while $P(z,t)$, equivalently $P(z,x)$ when the marker is named $x$, is the bivariate sign-polynomial series.

Summing the root-split recurrence over $m$ turns the two convolutions into products:

$$P(z,t)=z+zT(z)P(z,t)-tzT(tz)P(z,t).$$

Writing $q(z)=\sqrt{1-4z}$ and $s(z)=\sqrt{1-4tz}$, we obtain the compact algebraic function

$$\boxed{P(z,t)=\frac{2z}{2+q(z)-s(z)}}.$$

Thus the exponentially large family of trees has been compressed into coefficient extraction from one algebraic power series.

Fibonacci numbers are one coefficient in a quadratic ring

Work modulo $p$ in the quotient ring

$$R=\mathbb{F}_p[r]/(r^2-r-1).$$

Every element has the form $a+br$. The defining relation gives, for $i\ge1$,

$$r^i=F_{i-1}+F_ir,$$

and $r^0=1$ has $r$-coefficient $0=F_0$. Consequently the coefficient of $r$ in $P_m(r)$ is exactly the sum obtained when leaf $i$ receives $F_i$:

$$A(n)=[r]\,P_{n+1}(r).$$

No square root of $5$ and no extension-field inversion is needed. Addition, subtraction, and multiplication by $r$, $r+1$, $r-1$, or $2-r$ reduce to operations on the pair $(a,b)$:

$$\begin{aligned} r(a+br)&=b+(a+b)r,\\ (r+1)(a+br)&=(a+b)+(a+2b)r,\\ (r-1)(a+br)&=(b-a)+ar,\\ (2-r)(a+br)&=(2a-b)+(b-a)r. \end{aligned}$$

Rationalization gives a constant-memory coefficient recurrence

Now set $t=r$, $q(z)=\sqrt{1-4z}$, and $s(z)=\sqrt{1-4rz}$. Define

$$N(z)=1-q(z)+s(z)-q(z)s(z) +(r-1)z\bigl(q(z)+s(z)\bigr)+2(r+1)z.$$

Using $q(z)^2=1-4z$, $s(z)^2=1-4rz$, and $r^2=r+1$, rationalizing the denominator of $P$ yields

$$\bigl(4+2(2-r)^2z\bigr)P(z,r)=(2-r)N(z).$$

This identity is chosen because its denominator is only linear in $z$. It lets the implementation advance one coefficient at a time instead of storing a power-series table.

Factorial scaling removes divisions inside the loop

Let the factorial-scaled coefficients be

$$a_m=m![z^m]q(z),\quad b_m=m![z^m]s(z),\quad c_m=m![z^m]q(z)s(z),\quad u_m=m![z^m]P(z,r).$$

With $a_0=b_0=c_0=1$, $c_{-1}=0$, and $u_0=0$, coefficient comparison gives

$$\begin{aligned} a_m&=(4m-6)a_{m-1},\\ b_m&=(4m-6)r\,b_{m-1},\\ c_m&=2(2m-3)(r+1)c_{m-1} -16(m-3)(m-1)r\,c_{m-2}. \end{aligned}$$

The scaled coefficient of $N(z)$ is

$$d_m=-a_m+b_m-c_m +m(r-1)\bigl(a_{m-1}+b_{m-1}\bigr) +2(r+1)\delta_{m,1},$$

and the rationalized identity becomes the one-step update

$$\boxed{u_m=\frac{2-r}{4} \left(d_m-2m(2-r)u_{m-1}\right)}.$$

The loop runs through $m=n+1$. Finally it undoes the factorial scaling and extracts the root component:

$$A(n)=[r]\left(u_{n+1}\,((n+1)!)^{-1}\right)\pmod p.$$

Because $p$ is prime and $n+1=10000001\lt p$, both $4$ and $(n+1)!$ are nonzero modulo $p$. Their inverses are computed with Fermat's little theorem, $x^{-1}\equiv x^{p-2}\pmod p$.

Why the small brute-force check is valuable

The optimized recurrence is compact but algebraically dense. For $0\le n\le10$, the implementations independently enumerate every root split, construct every parenthesized value, and compare its direct sum with the ring recurrence. They also check the three published values. This catches sign errors at a right subtree, incorrect index shifts, and mistakes in either ring component before the large instance is evaluated.

Correctness Argument

Lemma 1. The recurrence for $P_m(t)$ sums every valid parenthesization exactly once. Every full binary tree has one root split $k$; the Catalan factors count the choices for the opposite subtree, and the right block acquires exactly the shift $t^k$ and the root minus sign.

Lemma 2. The generating function $P(z,t)$ has coefficient $P_m(t)$. Multiplying by $zC(z)$ and $tzC(tz)$ reproduces the two split convolutions of Lemma 1, so solving the resulting linear equation gives $2z/(2+q-s)$.

Lemma 3. In $R$, the $r$-coefficient of $r^i$ is $F_i$. Linearity therefore changes the formal leaf weights $t^i$ in $P_m(t)$ into the required Fibonacci weights when $t=r$.

Lemma 4. The updates for $a_m,b_m,c_m,d_m,u_m$ are coefficient identities of the square-root series and the rationalized equation. Factorial scaling is reversible because the relevant factorial is invertible modulo $p$.

By Lemmas 1 and 2, $P_{n+1}$ represents the sum over all valid expressions; by Lemma 3 its root component is $A(n)$; and by Lemma 4 the loop computes precisely that component. Hence the printed residue is $A(10^7)\bmod(10^9+9)$.

How the Code Works

alternating_sum keeps only the current $a_m$, $b_m$, the last two $c_m$ values, the previous $u_m$, and the running factorial. A ring element is stored as its constant and root components. The four short multiplication identities above are expanded directly in Python and Java to avoid allocating millions of temporary objects; they are the same operations performed by the C++ RingElement helpers.

mod_pow performs binary exponentiation for the two modular inverses. brute_values is deliberately separate from the optimized derivation: it recursively realizes all root splits for small inputs. run_checkpoints compares both methods for $n=0,\ldots,10$, then verifies the published residues.

Complexity Analysis

The main recurrence performs a constant number of modular ring operations for each $m=1,\ldots,n+1$, so it takes $O(n+\log p)$ time. The exponentiation term is negligible beside $n=10^7$. Only a fixed number of scalar residues is retained, giving $O(1)$ auxiliary memory.

Direct enumeration would require $C_n\sim4^n/(n^{3/2}\sqrt{\pi})$ expressions and is exponentially impossible. Even storing all coefficients through degree $n$ would use $O(n)$ memory; the linear-denominator recurrence avoids that storage.

Footnotes and References

Problem page: Project Euler 1007 - Alternating Difference
Catalan numbers: Wikipedia - Catalan number
Generating functions: Wikipedia - Generating function
Fibonacci sequence: Wikipedia - Fibonacci sequence
Quotient rings: Wikipedia - Quotient ring
Fermat's little theorem: Wikipedia - Fermat's little theorem

Problem 1007 source code

C++

#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;

constexpr i64 MOD = 1'000'000'009LL;
constexpr int TARGET = 10'000'000;

i64 normalize(i64 value) {
    value %= MOD;
    return value < 0 ? value + MOD : value;
}

i64 multiply_mod(const i64 lhs, const i64 rhs) {
    return lhs * rhs % MOD;
}

i64 mod_pow(i64 base, i64 exponent) {
    i64 result = 1;
    while (exponent > 0) {
        if ((exponent & 1LL) != 0) {
            result = multiply_mod(result, base);
        }
        base = multiply_mod(base, base);
        exponent >>= 1LL;
    }
    return result;
}

struct RingElement {
    i64 constant;
    i64 root;
};

RingElement add(const RingElement lhs, const RingElement rhs) {
    return {normalize(lhs.constant + rhs.constant), normalize(lhs.root + rhs.root)};
}

RingElement subtract(const RingElement lhs, const RingElement rhs) {
    return {normalize(lhs.constant - rhs.constant), normalize(lhs.root - rhs.root)};
}

RingElement scale(const RingElement value, const i64 factor) {
    return {multiply_mod(value.constant, factor), multiply_mod(value.root, factor)};
}

// Represents a + b*r modulo r^2-r-1; the r coefficient supplies Fibonacci weights.
RingElement multiply_by_root(const RingElement value) {
    return {value.root, normalize(value.constant + value.root)};
}

RingElement multiply_by_root_plus_one(const RingElement value) {
    return {normalize(value.constant + value.root),
            normalize(value.constant + 2 * value.root)};
}

RingElement multiply_by_root_minus_one(const RingElement value) {
    return {normalize(value.root - value.constant), value.constant};
}

RingElement multiply_by_two_minus_root(const RingElement value) {
    return {normalize(2 * value.constant - value.root),
            normalize(value.root - value.constant)};
}

i64 alternating_sum(const int n) {
    const i64 inverse_four = mod_pow(4, MOD - 2);
    i64 a_scaled = 1;
    RingElement b_scaled{1, 0};
    RingElement c_older{0, 0};
    RingElement c_previous{1, 0};
    RingElement result_scaled{0, 0};
    i64 factorial = 1;

    for (i64 m = 1; m <= static_cast<i64>(n) + 1; ++m) {
        const i64 linear_factor = normalize(4 * m - 6);
        const i64 a_next = multiply_mod(linear_factor, a_scaled);
        const RingElement b_next = scale(multiply_by_root(b_scaled), linear_factor);

        RingElement c_next =
            scale(multiply_by_root_plus_one(c_previous), normalize(2 * (2 * m - 3)));
        const i64 quadratic_factor =
            multiply_mod(normalize(16 * normalize(m - 3)), normalize(m - 1));
        c_next = subtract(c_next, scale(multiply_by_root(c_older), quadratic_factor));

        RingElement numerator{normalize(-a_next), 0};
        numerator = add(numerator,
                        scale(multiply_by_root_minus_one({a_scaled, 0}), m % MOD));
        numerator = add(numerator, b_next);
        numerator = add(numerator,
                        scale(multiply_by_root_minus_one(b_scaled), m % MOD));
        numerator = subtract(numerator, c_next);
        if (m == 1) {
            numerator = add(numerator, {2, 2});
        }

        const RingElement adjusted =
            subtract(numerator,
                     scale(multiply_by_two_minus_root(result_scaled), (2 * m) % MOD));
        const RingElement result_next =
            scale(multiply_by_two_minus_root(adjusted), inverse_four);

        a_scaled = a_next;
        b_scaled = b_next;
        c_older = c_previous;
        c_previous = c_next;
        result_scaled = result_next;
        factorial = multiply_mod(factorial, m % MOD);
    }

    return multiply_mod(result_scaled.root, mod_pow(factorial, MOD - 2));
}

std::vector<i64> brute_values(const std::vector<i64>& fibonacci,
                              const int begin,
                              const int end) {
    if (end - begin == 1) {
        return {fibonacci[static_cast<std::size_t>(begin)]};
    }

    std::vector<i64> values;
    for (int split = begin + 1; split < end; ++split) {
        const std::vector<i64> left = brute_values(fibonacci, begin, split);
        const std::vector<i64> right = brute_values(fibonacci, split, end);
        for (const i64 lhs : left) {
            for (const i64 rhs : right) {
                values.push_back(lhs - rhs);
            }
        }
    }
    return values;
}

i64 brute_alternating_sum(const int n) {
    std::vector<i64> fibonacci(static_cast<std::size_t>(n + 1), 0);
    if (n >= 1) {
        fibonacci[1] = 1;
    }
    for (int index = 2; index <= n; ++index) {
        fibonacci[static_cast<std::size_t>(index)] =
            fibonacci[static_cast<std::size_t>(index - 1)] +
            fibonacci[static_cast<std::size_t>(index - 2)];
    }

    i64 total = 0;
    for (const i64 value : brute_values(fibonacci, 0, n + 1)) {
        total = normalize(total + value);
    }
    return total;
}

void require_checkpoint(const bool condition, const std::string& description) {
    if (!condition) {
        std::cerr << "Checkpoint failed: " << description << '\n';
        std::exit(EXIT_FAILURE);
    }
}

void run_checkpoints() {
    for (int n = 0; n <= 10; ++n) {
        require_checkpoint(alternating_sum(n) == brute_alternating_sum(n),
                           "brute force comparison for n=" + std::to_string(n));
    }

    require_checkpoint(alternating_sum(3) == normalize(-6), "published A(3)");
    require_checkpoint(alternating_sum(10) == normalize(-177'666), "published A(10)");
    require_checkpoint(alternating_sum(100) == 71'792'794, "published A(100)");
}

}  // namespace

int main() {
    run_checkpoints();
    std::cout << alternating_sum(TARGET) << '\n';
    return 0;
}

Python

#!/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()

Java

import java.util.ArrayList;
import java.util.List;

public class Euler1007 {
    private static final long MOD = 1_000_000_009L;
    private static final int TARGET = 10_000_000;

    private static long normalize(long value) {
        value %= MOD;
        return value < 0 ? value + MOD : value;
    }

    private static long multiplyMod(long lhs, long rhs) {
        return lhs * rhs % MOD;
    }

    private static long modPow(long base, long exponent) {
        long result = 1;
        while (exponent > 0) {
            if ((exponent & 1L) != 0) {
                result = multiplyMod(result, base);
            }
            base = multiplyMod(base, base);
            exponent >>= 1;
        }
        return result;
    }

    private static long alternatingSum(int n) {
        long inverseFour = modPow(4, MOD - 2);

        // Every quotient-ring value is represented by constant + root * r,
        // with r^2 = r + 1.  Keeping the components in primitive locals avoids
        // allocating temporary objects in the ten-million-step loop.
        long aScaled = 1;
        long bConstant = 1;
        long bRoot = 0;
        long cOlderConstant = 0;
        long cOlderRoot = 0;
        long cPreviousConstant = 1;
        long cPreviousRoot = 0;
        long resultConstant = 0;
        long resultRoot = 0;
        long factorial = 1;

        for (long m = 1; m <= (long) n + 1; ++m) {
            long linearFactor = normalize(4 * m - 6);
            long mMod = m % MOD;

            long aNext = multiplyMod(linearFactor, aScaled);

            // r(a + br) = b + (a + b)r.
            long bNextConstant = multiplyMod(linearFactor, bRoot);
            long bNextRoot = multiplyMod(linearFactor, normalize(bConstant + bRoot));

            // (r + 1)(a + br) = (a + b) + (a + 2b)r.
            long cNextConstant = multiplyMod(
                    linearFactor, normalize(cPreviousConstant + cPreviousRoot));
            long cNextRoot = multiplyMod(
                    linearFactor, normalize(cPreviousConstant + 2 * cPreviousRoot));

            long quadraticFactor = multiplyMod(
                    normalize(16 * normalize(m - 3)), normalize(m - 1));
            // Subtract quadraticFactor * r * cOlder.
            cNextConstant = normalize(
                    cNextConstant - multiplyMod(quadraticFactor, cOlderRoot));
            cNextRoot = normalize(
                    cNextRoot
                            - multiplyMod(
                                    quadraticFactor,
                                    normalize(cOlderConstant + cOlderRoot)));

            // Build n_m.  Here (r-1)(a+br) = (b-a) + ar.
            long numeratorConstant = normalize(
                    -aNext
                            - multiplyMod(mMod, aScaled)
                            + bNextConstant
                            + multiplyMod(mMod, normalize(bRoot - bConstant))
                            - cNextConstant);
            long numeratorRoot = normalize(
                    multiplyMod(mMod, aScaled)
                            + bNextRoot
                            + multiplyMod(mMod, bConstant)
                            - cNextRoot);
            if (m == 1) {
                numeratorConstant = normalize(numeratorConstant + 2);
                numeratorRoot = normalize(numeratorRoot + 2);
            }

            // Subtract 2m(2-r)u_{m-1};
            // (2-r)(a+br) = (2a-b) + (b-a)r.
            long twiceM = 2 * mMod % MOD;
            long adjustedConstant = normalize(
                    numeratorConstant
                            - multiplyMod(
                                    twiceM,
                                    normalize(2 * resultConstant - resultRoot)));
            long adjustedRoot = normalize(
                    numeratorRoot
                            - multiplyMod(
                                    twiceM,
                                    normalize(resultRoot - resultConstant)));

            // u_m = (2-r) * adjusted / 4.
            long resultNextConstant = multiplyMod(
                    normalize(2 * adjustedConstant - adjustedRoot), inverseFour);
            long resultNextRoot = multiplyMod(
                    normalize(adjustedRoot - adjustedConstant), inverseFour);

            aScaled = aNext;
            bConstant = bNextConstant;
            bRoot = bNextRoot;
            cOlderConstant = cPreviousConstant;
            cOlderRoot = cPreviousRoot;
            cPreviousConstant = cNextConstant;
            cPreviousRoot = cNextRoot;
            resultConstant = resultNextConstant;
            resultRoot = resultNextRoot;
            factorial = multiplyMod(factorial, mMod);
        }

        return multiplyMod(resultRoot, modPow(factorial, MOD - 2));
    }

    private static List<Long> bruteValues(
            List<Long> fibonacci, int begin, int end) {
        if (end - begin == 1) {
            return List.of(fibonacci.get(begin));
        }

        List<Long> values = new ArrayList<>();
        for (int split = begin + 1; split < end; ++split) {
            List<Long> left = bruteValues(fibonacci, begin, split);
            List<Long> right = bruteValues(fibonacci, split, end);
            for (long lhs : left) {
                for (long rhs : right) {
                    values.add(lhs - rhs);
                }
            }
        }
        return values;
    }

    private static long bruteAlternatingSum(int n) {
        List<Long> fibonacci = new ArrayList<>();
        for (int i = 0; i <= n; ++i) {
            fibonacci.add(0L);
        }
        if (n >= 1) {
            fibonacci.set(1, 1L);
        }
        for (int index = 2; index <= n; ++index) {
            fibonacci.set(index, fibonacci.get(index - 1) + fibonacci.get(index - 2));
        }

        long total = 0;
        for (long value : bruteValues(fibonacci, 0, n + 1)) {
            total = normalize(total + value);
        }
        return total;
    }

    private static void requireCheckpoint(boolean condition, String description) {
        if (!condition) {
            throw new AssertionError("checkpoint failed: " + description);
        }
    }

    private static void runCheckpoints() {
        for (int n = 0; n <= 10; ++n) {
            requireCheckpoint(
                    alternatingSum(n) == bruteAlternatingSum(n),
                    "brute force comparison for n=" + n);
        }

        requireCheckpoint(alternatingSum(3) == normalize(-6), "published A(3)");
        requireCheckpoint(alternatingSum(10) == normalize(-177_666), "published A(10)");
        requireCheckpoint(alternatingSum(100) == 71_792_794, "published A(100)");
    }

    public static void main(String[] args) {
        runCheckpoints();
        System.out.println(alternatingSum(TARGET));
    }
}