Integer Sequence Research

Autonomous Integer Sequence Research System

Adversarial + generative pipeline for OEIS-grade sequence discovery. The system generates candidates, cross-validates them, aggressively tries to break conjectured formulas, and scores survivors on a confidence scale.

Reference implementation: seq_research.py
Run with python seq_research.py (full sweep) or
python seq_research.py --single <NAME> (single sequence).

Pipeline Overview

Generator → Evaluator → Conjecturer → Falsifier → Classifier → Report
                ↑_____________________________|
                     (loop until validated or discarded)

Each stage is described below. Read all sections before starting work; the adversarial falsifier (§4) is the most important stage and the one most commonly skipped too early.

§1 · Sequence Generation

Grammar

Build candidates compositionally from:

Bias Rules

Prefer candidates that are likely to be multiplicative — they are richer, easier to characterize at prime powers, and more likely to be novel. Structural signals to pursue:

• φ compositions over divisors
• Products involving unitary divisors
• Dirichlet convolutions (μ ★ f, id ★ φ, etc.)
• Functions of the prime signature only (depend on exponents, not primes)

Avoid: trivial linear combinations of known sequences, constant multiples of φ or σ, anything obviously reducible by inspection.

Example Constructions Worth Exploring

φ(∏_{d|n} φ(d))          — §2 shows this is non-multiplicative, interesting growth
∑_{d|n} d·φ(d)            — multiplicative, near-quadratic growth, strong candidate
∏_{d|n, gcd(d,n/d)=1} φ(d)  — unitary variant; check against dual impl
φ(n) · 2^{ω(n)}           — multiplicative; geometric prime-power profile
(μ ★ φ)(n)                 — Dirichlet convolution; vanishes on many n
∑_{d|n} gcd(d, n/d)        — multiplicative; sub-linear growth
∏_{p^a ∥ n} p^{φ(a)}      — depends only on prime signature

§2 · Multi-Implementation Evaluation

Rule: every sequence must have ≥2 independent implementations before conjecture work begins.

Implementation Types

A. Definition-based — iterate over divisors(n) directly.

B. Factorization-based — exploit the prime factorization n = ∏ pᵢ^aᵢ. For multiplicative functions this gives a cleaner formula: f(n) = ∏ f(pᵢ^aᵢ), computable from factorint(n) alone.

C. Hybrid — cached + partial symbolic where helpful for large n.

Differential Validation

Run both implementations on n = 1…30 plus the weak zones (§4). If they disagree at any n: the sequence definition is ambiguous or one implementation is wrong. Fix before proceeding. Do not conjecture on a sequence with unresolved implementation disagreements.

§3 · Multiplicativity Detection

Test f(ab) = f(a)·f(b) for gcd(a,b) = 1 across:

• Pairs of distinct primes: (2,3), (2,5), (3,7), (5,11), …
• Prime × prime square: (2, 9), (4, 3), (4, 25), …
• Larger coprime products: (4, 9), (8, 25), (16, 27), …

If multiplicative: proceed to prime-power profiling (§3.1) and enforce multiplicative modeling in conjecture work (§5).

If not multiplicative: compute the interaction residual f(ab)/(f(a)·f(b)) for coprime pairs and check whether the residual has structure. If it does, the function may be "almost multiplicative" in a useful sense. If it doesn't, the function is unlikely to have a clean closed form.

3.1 · Prime-Power Profiling

For any multiplicative (or candidate-multiplicative) function, compute:

f(p^k) for p ∈ {2, 3, 5, 7, 11}, k = 1…8

Look for:

Once the prime-power formula is identified, the full multiplicative function follows immediately from f(n) = ∏ f(p^a) over the factorization.

§4 · Adversarial Falsification (Most Important Stage)

The goal is to find the smallest n where the conjecture fails before claiming success. Always run this before reporting any formula.

4.1 · Weak Zones (Test These First)

These inputs break naive formulas most often:

WEAK_ZONES = (
    # Prime powers
    [p**k for p in [2,3,5,7,11,13] for k in range(1,6)] +
    # Products of two primes and their powers
    [p**a * q**b for p,q in [(2,3),(2,5),(3,5)] for a,b in [(1,1),(2,1),(1,2)]] +
    # Highly composite: 12, 24, 60, 120, 360, 720, 840, 2520
    [12, 24, 60, 120, 360, 720, 840, 2520] +
    # Powers of 2 and near-powers
    [2**k for k in range(1,20)] + [2**k - 1 for k in range(2,20)] +
    # Three-prime products
    [2*3*5, 2*3*7, 2*5*7, 3*5*7, 2*3*5*7]
)

4.2 · Mutation Attack

Given a failing n₀, propagate:

n₀ → n₀ · p   (extend by small prime)
n₀ → n₀ / p   (reduce if p² | n₀)
n₀ → n₀²      (square)

This often reveals whether the failure is isolated or structural.

4.3 · Delta Debugging

If a counterexample is found at composite n, find the smallest n' where the formula fails. Prime powers and small semiprimes are the most informative.

4.4 · Red Team Checklist

Before accepting any formula, verify:

1. Does it work for n = 1? (edge case: φ(1)=1, σ(1)=1, ω(1)=0, μ(1)=1)
2. Does it work for prime n? For p²? For p³?
3. Does it work for pq (two distinct primes)?
4. Does it work for highly composite n (60, 120, 720)?
5. Is the formula actually a known sequence in disguise? (→ §6 OEIS check)
6. Does the conjectured formula match the factorization-based dual? (→ §2)
7. Is the growth rate consistent across the tested range? (→ §5)

§5 · Conjecture Engine

Only conjecture after multiplicativity testing (§3) and adversarial testing (§4) have passed.

5.1 · Template Hierarchy

Try formulas in this order (simpler first):

1. f(n) = n, f(n) = φ(n), f(n) = σ(n), f(n) = τ(n), f(n) = rad(n)
2. f(n) = φ(n)^a, f(n) = n^a · φ(n)^b
3. f(n) = 2^{g(ω(n))} · h(n)
4. For multiplicative f: derive from the prime-power formula directly

5.2 · Log-Space Regression

For growth estimation and exponent fitting:

alpha = cov(log n, log f(n)) / var(log n)

This gives the asymptotic exponent in f(n) ~ C · n^alpha. Use n > 10 to avoid distortion from small-n irregularities.

5.3 · Multiplicative Conjecture Construction

If f is multiplicative with prime-power formula f(p^k) = g(p,k), then:

f(n) = ∏_{p^a ∥ n}  g(p, a)

This is the canonical form. Test the fully reconstructed multiplicative function against the definition-based implementation on n = 1…100 and all weak zones.

§6 · OEIS Collision Filter

Compute the first 100 terms and check against:

• Standard arithmetic functions: φ, σ, τ, rad, ω, Ω, μ, id, id²
• Their Dirichlet convolutions
• Known transforms: partial sums, Euler transform, Möbius transform

If the first 30 terms match a known sequence exactly, the candidate is not novel — report the collision and move on.

If the first 30 terms are unique: compute 100 terms and search OEIS directly (format: 1, 2, 4, 4, 8, 8, 12, 8, 12, 16 — comma-separated, no spaces around commas).

§7 · Confidence Scoring

Each sequence receives a score 0–100:

Confidence levels:

§8 · Output Standard (OEIS-Grade)

A sequence report must include:

1. Definition — at least two equivalent formulations
2. First 100–500 terms — verified across ≥2 implementations
3. Multiplicativity — proved or disproved with witness if not
4. Prime-power characterization — f(p^k) formula with evidence
5. Closed form — if found, with falsification report (how many n tested, result)
6. Counterexample search — what was tested, what survived
7. Growth exponent — asymptotic α from log-regression
8. OEIS collision check — which known sequences were compared
9. Confidence score — with breakdown
10. Reference implementation — clean Python, importable

§9 · Known Failure Modes and Mitigations

§10 · Meta-Insight Tracking

After running multiple candidates, record which constructions tend to produce useful sequences and which are dead ends. Updated after Batch 1 and Batch 2.

Productive:

• φ(n)·c^{ω(n)} family — multiplicative for any constant c; f(p^k) = c·p^{k-1}(p-1); c=2-6 confirmed OEIS-READY
• n·2^{Ω(n)} — "totally doubled" variant of n; f(p^k) = (2p)^k; exact closed form
• Σ_{d|n} d²·φ(n/d) — Dirichlet id²★φ; f(p^k) = p^{2k} + p^{k-1}(p^k - 1); near-quadratic (α≈2)
• Σ_{d|n} φ(d)² — Dirichlet φ²★1; f(p^k) = 1+(p-1)²(p^{2k}-1)/(p²-1); near-quadratic (α≈2)
• Σ_{d|n} gcd(d, n/d) — multiplicative; piecewise formula for even/odd exponents; slow growth (α≈0.2)
• ∏_{p^a ∥ n} φ(p^a + 1) — product of φ of prime-power+1; α≈0.8; OEIS-READY
• Completely additive g(a) family (Σ_{p^a ∥ n} g(a)); prime-signature only; slow growth (α≈0.1) · g(a) = a(a+1)/2 — triangular of exponent; f(p^k) = k(k+1)/2 (same all p) · g(a) = 2^a — exponential exponent; f(p^k) = 2^k (same all p) · g(a) = a² — squared exponent; f(p^k) = k² (same all p)
• Dirichlet convolutions involving μ and φ (when convergent)
• Products over unitary divisors (when α < 3.0)

Dead ends:

• φ ∘ (∏_{d|n} φ(d)) — non-multiplicative, hard to characterize
• gcd(σ(n), φ(n)) — non-multiplicative, already in OEIS
• lcm(σ(n), φ(n)) — non-multiplicative, complex growth
• σ(rad(n)), rad(σ(n)), Σ rad(d) — sympy's rad is symbolic (radians!), not integer radical; implement manually: rad(n) = prod(p for p in factorint(n))
• φ(n)·σ(n)/n — only integer at specific n; not a well-defined integer sequence
• Σ_{d|n, d sqfree} φ(d) = rad(n) — A007947; see §11
• Σ_{d|n, d sqfree} d = σ(rad(n)) = ∏_{p|n}(1+p) — derivable from known; see §11
• ∏_{d unitary|n} (d+1) — α>3.0 (too fast, out of range)

Watch list (experimental, unresolved):

• ∑_{d|n} (d XOR n/d) — non-multiplicative but interesting XOR structure
• ω(n)^{φ(n)} mod n — slow growth, heavily prime-power-dominated
• φ(n)·ω(n) — not multiplicative, but residual f(ab)/(f(a)f(b)) is governed exactly by the harmonic mean of ω(a) and ω(b); might be worth characterizing as a type
• Σ_{d unitary} φ(d) = ∏(1 + φ(p^a)) — clean but α≈2.9, too fast growth
• (μ ★ φ)(n) = φ(n) — Möbius inversion, known, score 55 (weak zones problematic)
• ∏_{d unitary} φ(d) — non-multiplicative, α≈2.9, experimental

§11 · Confirmed Classical Identities (Do Not Resubmit)

These were rediscovered by the pipeline and confirmed as known. Do not submit.

§12 · Parameterized Families

When a sequence is confirmed as part of a parameterized family, note which parameter values are already in OEIS to guide novelty targeting.

Family: φ(n)·c^{ω(n)} for integer c ≥ 1

Multiplicative for any constant c. Prime-power formula: f(p^k) = c·p^{k-1}(p-1).

Family: Σ_{p^a ∥ n} g(a) for various g (completely additive, prime-signature functions)

Family: Σ_{d|n} φ(d)^k for k≥1

Family: Σ_{d|n} d^k·φ(n/d)

When testing a new member of a known family, reduce the novelty check burden by verifying only the first 30 terms against OEIS before full pipeline.

§13 · Combinatorics Research

Same generator → falsifier → classifier loop as §1–§5, applied to counting sequences instead of arithmetic functions. The key difference: ground truth usually comes from direct enumeration on small n, not from a second symbolic formula, so brute-force enumeration is the "definition-based" implementation (§2.A) and the closed form is the "factorization-based" one (§2.B).

13.1 · Grammar

13.2 · Ground-Truth-First Workflow

1. Brute-force enumerate the combinatorial objects for small n (n ≤ 12–15 depending on growth rate) — generate every object explicitly, don't just trust a recursive count.
2. Derive or guess a closed form / generating function from the counted values.
3. Cross-check the closed form against the brute-force counts on the same range before extending further (mirrors §2's differential validation).
4. Only then extend via the closed form or recurrence to get 100+ terms for the OEIS collision check (§6).

13.3 · Adversarial Checks Specific to Combinatorics

• Off-by-one on n=0: does the sequence start counting at the empty structure? Half of all combinatorics bugs live here.
• Symmetry checks: many combinatorial arrays are symmetric or palindromic by row (Narayana numbers, Eulerian numbers) — a broken symmetry is a fast falsifier.
• Row-sum checks: if the object is a triangle (Stirling-like, Narayana), verify row sums against a known sequence (e.g. Narayana row sums = Catalan, Eulerian row sums = n!).
• Generating function sanity: if a GF is conjectured, expand it symbolically (sympy.series) to as many terms as were enumerated and diff against the brute-force list.

13.4 · Known Identities (Do Not Resubmit)

§14 · Prime Numbers and Prime Gaps Research

14.1 · Core Objects

14.2 · Generation and Ground Truth

Use sympy.nextprime / sympy.prevprime / sympy.primerange for ground truth — never hand-roll a primality sieve for validation purposes, only for performance-critical generation that is then cross-checked against sympy on an overlapping range (mirrors §2's dual-implementation rule).

from sympy import primerange, isprime, nextprime

def prime_gaps(limit):
    primes = list(primerange(2, limit))
    return [primes[i+1] - primes[i] for i in range(len(primes) - 1)]

14.3 · Falsification Specific to Prime Gaps

• Record claims are the highest-risk output. A "new maximal gap" claim must be verified by regenerating the prime list independently (e.g. via Miller–Rabin at high witness count, or a second sieve implementation) — do not report a gap record from a single code path.
• Compare against the Cramér model (expected max gap near x scales as (ln x)^2) to sanity-check whether an observed gap is plausible or a bug.
• Off-by-one at the boundary: verify gap counting includes/excludes the correct endpoints; check the first few known values (2,1,2,2,4,2,4,2,4,6...) against A001223 before trusting a bulk computation.
• Twin-prime-like patterns: verify the constellation pattern is admissible (no residue class mod every small prime is fully blocked) before searching for it — inadmissible patterns provably have finitely many (usually zero) instances and will silently return empty or garbage.

14.4 · Verified Reference Records (Dogfood-Confirmed)

Independently regenerated via sympy.primerange up to 300,000 and cross-checked against the pipeline's own gap-tracking logic; use this as a fast local sanity table before trusting a bulk maximal-gap computation elsewhere in this range:

(p, gap, merit=gap/ln(p)):
(2,1,1.443)  (3,2,1.821)   (7,4,2.056)   (23,6,1.914)   (89,8,1.782)
(113,14,2.962) (523,18,2.876) (887,20,2.946) (1129,22,3.130) (1327,34,4.728)
(9551,36,3.928) (15683,44,4.555) (19609,52,5.261) (31397,72,6.954) (155921,86,7.192)

All 15 match A005250/A005669's published prefix. If a new run disagrees with this table anywhere in p ≤ 155921, the bug is in the new run, not here.

14.5 · Worked Constellation Example (Dogfood-Confirmed)

Pattern (p, p+2, p+6, p+8) (prime quadruplets) tested for admissibility (§14.3) and then actually searched, rather than left as a described-but- unexercised procedure:

• Admissibility: no prime q ≤ 50 blocks every residue class → pattern is admissible.
• Real search to 2,000,000: 295 quadruplets found; first few start at p = 5, 11, 101, 191, 821, 1481, 1871, 2081, 3251, 3461.
• Hardy–Littlewood cross-check: numerically integrating the singular- series estimate (C₄ ≈ 4.1512) over the same range gives ≈365, so the actual/estimate ratio is ≈0.81 — the right order of magnitude for a heuristic asymptotic at this range (convergence is slow; treat as a plausibility check per §14.4's Statistical/Heuristic caveat, not a precision match).

14.6 · Statistical / Heuristic Checks (Not Proof)

Use these to flag "interesting" candidates for further exploration, never as a substitute for exact verification:

• Gap distribution vs. Cramér–Granville heuristic model
• Hardy–Littlewood constant estimates for k-tuple density
• Bertrand's postulate as a trivial sanity floor (p(n+1) < 2 p(n) must always hold — if a generated table violates this, the table is wrong)

§15 · PSLQ / Integer Relation Search

Used to discover (or rule out) linear integer relations among a vector of real constants — e.g. testing whether some combination of ζ(3), π², ln 2, Catalan's constant G, and γ satisfies Σ c_i x_i = 0 for integers c_i. This is the constants-analogue of §1–§6: generate candidate constant vectors, run the relation search, then adversarially confirm or falsify.

15.1 · Tooling

from mpmath import mp, mpf, pslq, zeta, pi, euler, catalan, log

mp.dps = 100  # decimal digits of precision; see §15.3

vec = [zeta(3), pi**2, log(2), catalan, euler, mpf(1)]
relation = pslq(vec, maxsteps=10**6)

pslq returns None if no relation is found within the given precision and bound, or a list of integer coefficients c_i such that Σ c_i · vec[i] ≈ 0 to the working precision.

15.2 · Candidate Constant Construction

• Mix constants from different "families" (zeta values, logs of small integers, Catalan-type constants, algebraic numbers, powers of π) — same-family combinations mostly rediscover trivial/known identities. Search first for a term like daedalus/skills history of resubmitted relations before reporting.
• Include a constant 1 in the vector so PSLQ can find relations with a free rational/integer term, not only homogeneous ones.
• Prefer small vectors (4–8 constants) — PSLQ's cost and false-positive risk both grow with vector length.

15.3 · Precision Is the Entire Ballgame

The single most important rule for this section: a relation found at precision P digits is only credible if a coefficient vector with max|c_i| bits requires roughly n · max|c_i-bits| digits of precision to distinguish from a numerical coincidence (n = vector length). Rules of thumb:

1. Always run PSLQ at two different precisions (e.g. 50 and 100 digits). If the same relation (same integer vector, possibly up to sign/scale) appears at both, it survives the first falsification pass.
2. After finding a candidate relation, substitute it back and evaluate the residual at 2–3× the discovery precision. The residual must shrink roughly in proportion to the added precision (i.e., actually go to zero), not merely stay "small."
3. Reject any relation whose smallest nonzero coefficient is implausibly large relative to the precision used — that is the classic PSLQ false-positive signature (an artifact of insufficient precision, not a real identity).
4. Treat rediscovery of a known identity (e.g. Euler's reflection formula, standard zeta value relations from §20) as a successful validation of the pipeline, not a novel result — check §20.4 and standard constant tables before claiming novelty.

15.4 · Adversarial Falsification

• Run PSLQ on a vector of algebraically independent constants (e.g. π, e, ln 2, ln 3) as a negative control — it must return None (or only the trivial/garbage relation at insufficient precision). If it "finds" a relation here, precision is too low or the implementation is buggy.
• Randomize the order of the input vector; a real relation's coefficients permute consistently, a precision artifact often doesn't reproduce.

15.5 · Confirmed Negative Results (Dogfood Reference)

Broader constant vectors actually run (not just described) at 50 and 100 digits, both returning no relation at either precision — recorded so a future pass doesn't re-spend a search budget rediscovering the same absence of a low-degree relation:

A clean None at two precisions for a plausible-looking vector is itself a useful (negative) result: it's evidence against a simple low-degree relation existing among those specific constants at that vector length, not just an untested gap. Widening the vector (more constants) or its precision is the natural next step if this combination is revisited.

§16 · Stirling Numbers

16.1 · Definitions

16.2 · Recurrences (Ground Truth A)

c(n,k) = c(n-1,k-1) + (n-1)·c(n-1,k)     # unsigned first kind
s(n,k) = s(n-1,k-1) - (n-1)·s(n-1,k)      # signed first kind
S(n,k) = S(n-1,k-1) + k·S(n-1,k)          # second kind

Base cases: c(0,0)=s(0,0)=S(0,0)=1; c(n,0)=s(n,0)=S(n,0)=0 for n>0; c(n,k)=s(n,k)=S(n,k)=0 for k>n.

16.3 · Closed Form (Ground Truth B, Second Kind Only)

S(n,k) = (1/k!) · Σ_{i=0}^{k} (-1)^i · C(k,i) · (k-i)^n

Use this as the independent cross-check against the recurrence (§2-style dual implementation) — sympy exposes both via sympy.functions.combinatorial.numbers.stirling.

16.4 · Known Row/Column Identities (Do Not Resubmit)

16.5 · Falsification Checklist

1. Verify recurrence and closed form (second kind) agree for n=1…20, all k.
2. Verify row sums against factorial / Bell numbers.
3. Check sign convention explicitly before submitting — mixing up s(n,k) (signed) and c(n,k) (unsigned) is the single most common submission error for this family.

§17 · Bernoulli Numbers

17.1 · Definition

Generating function:

t / (e^t - 1) = Σ_{n=0}^∞ B_n · t^n / n!

Sign convention warning: this convention gives B_1 = -1/2. The "other" convention (t/(1-e^{-t})) gives B_1 = +1/2. State the convention explicitly in any report — this is the #1 source of disagreement when cross-checking against a second source.

Dogfood-confirmed: implementing the §17.2 recurrence directly and diffing against sympy.bernoulli(n) for n=0..20 produces exactly one disagreement — at n=1 (-1/2 from the recurrence vs sympy's +1/2) — with n=0 and all n=2..20 agreeing exactly (including B₁₂). This is the textbook signature of a sign-convention mismatch rather than an implementation bug: a real bug would also perturb the even-indexed terms or fail von Staudt–Clausen (§17.4), which it does not. Treat a lone n=1 disagreement as "check your convention," not "something is broken."

17.2 · Recurrence (Ground Truth A)

Σ_{k=0}^{n} C(n+1, k) · B_k = 0   for n ≥ 1,   B_0 = 1

Solve for B_n by isolating the k=n term.

17.3 · Closed Form via Zeta (Ground Truth B)

B_{2n} = (-1)^{n+1} · 2 · (2n)! · ζ(2n) / (2π)^{2n}     for n ≥ 1

Cross-check numerator/denominator against sympy.bernoulli(n) (exact rational) — this ties directly into §20 (Riemann zeta).

17.4 · Structural Facts (Use as Falsifiers, Not Just Trivia)

• Odd-index vanishing: B_n = 0 for all odd n ≥ 3. Any nonzero "novel" odd-index Bernoulli-like sequence should be treated as a bug or a different (non-Bernoulli) object.
• von Staudt–Clausen theorem: the denominator of B_{2n} is Π_{(p-1) | 2n} p (product over primes p such that p-1 divides 2n). Use this to verify a computed denominator without needing the full numerator.
• Irregular primes: a prime p is irregular if it divides the numerator of some B_{2k} for 2k ≤ p-3. Relevant if the research touches Bernoulli numerators as a sequence in their own right (A092132 numerators, A002445 denominators) — check against known irregular prime tables (37, 59, 67, 101, 103, ... ) before claiming a new pattern.

17.5 · Weak Zones

Test n = 0,1,2,3,4,6,8,10,12 explicitly — B_{12} = -691/2730 is the classical stress test because 691 is the first "irregular-adjacent" numerator prime that trips up naive implementations using floating point instead of exact rationals.

§18 · Factorial Decomposition

18.1 · Legendre's Formula (Ground Truth)

The exponent of prime p in the factorization of n!:

e_p(n!) = Σ_{i=1}^{∞} floor(n / p^i)          (finite sum, stops when p^i > n)
        = (n - s_p(n)) / (p - 1)               (s_p(n) = digit sum of n in base p)

Both forms must agree — use the digit-sum form as the independent cross-check (§2-style dual implementation) against the floor-sum form.

def legendre_exponent(n, p):
    e, pk = 0, p
    while pk <= n:
        e += n // pk
        pk *= p
    return e

def digit_sum_base(n, p):
    s = 0
    while n:
        s += n % p
        n //= p
    return s

def legendre_via_digitsum(n, p):
    return (n - digit_sum_base(n, p)) // (p - 1)

18.2 · Derived Sequences

18.3 · Falsification Checklist

1. Verify e_2(n!) ≥ e_p(n!) for all p > 2 at the same n (2 is always the most abundant prime factor of a factorial) — a violation means a bug.
2. Cross-check trailing-zero counts against direct string manipulation (str(math.factorial(n)).rstrip('0')) for n up to a few hundred before trusting the Legendre-based formula for larger n.
3. Weak zones: n = powers of p and p-1 below/above a power (e.g. for p=5: n = 24, 25, 26, 124, 125, 126) — these are where floor(n/p^i) terms change discontinuously and where off-by-one bugs surface.

§19 · Harmonic Numbers

19.1 · Definitions

H_n       = Σ_{k=1}^{n} 1/k                  (ordinary harmonic number)
H_n^{(m)} = Σ_{k=1}^{n} 1/k^m                 (generalized / order-m harmonic number)

H_n^{(1)} = H_n. As m → ∞ the generalized harmonic numbers relate directly to ζ(m) (§20) since H_n^{(m)} → ζ(m) as n → ∞ for m > 1.

19.2 · Exact Arithmetic Is Mandatory

Always compute with fractions.Fraction (or sympy.Rational), never floats — harmonic number numerators/denominators are OEIS objects in their own right (A001008 numerators, A002805 denominators), and float accumulation silently corrupts the denominator's prime structure.

from fractions import Fraction

def harmonic(n, m=1):
    return sum(Fraction(1, k**m) for k in range(1, n + 1))

19.3 · Structural Falsifiers

• Wolstenholme's theorem: for prime p ≥ 5, the numerator of H_{p-1} is divisible by p², and the numerator of the second-order sum Σ_{k=1}^{p-1} 1/k² is divisible by p. Use this as a targeted correctness check at p = 5, 7, 11, 13, ... — a computed numerator that fails this divisibility signals a bug, not a mathematical exception.
• Denominator growth: denominator(H_n) | lcm(1, 2, ..., n), with equality failing only at specific n — verify this containment rather than assuming exact equality.
• Digamma connection: H_n = γ + ψ(n+1) where ψ is the digamma function — useful as a high-precision numeric cross-check via mpmath.digamma independent of the exact-fraction computation.

19.4 · Known Identities (Do Not Resubmit)

§20 · Riemann Zeta Function

20.1 · Scope and Caution

This section covers numerical and special-value work with ζ(s) — computing special values, verifying identities, and locating zeros to high precision. It does not cover, and this skill should not be used to claim, progress on the Riemann Hypothesis itself. Zero-location work is purely numerical verification (checking known zeros / hunting for counterexamples to RH within a numerically verified range), never a symbolic proof attempt.

20.2 · Special Values (Ground Truth via Known Closed Forms)

ζ(2)  = π²/6
ζ(4)  = π⁴/90
ζ(2k) = (-1)^{k+1} · B_{2k} · (2π)^{2k} / (2 · (2k)!)     (ties to §17.3)
ζ(-1) = -1/12
ζ(-2k) = 0   for k ≥ 1  (trivial zeros)
ζ(0)  = -1/2

Cross-check any computed ζ(2k) against the Bernoulli-number closed form (§17.3) — this is the primary dual-implementation check for even integer arguments.

20.3 · Numerical Tooling

from mpmath import mp, zeta, mpc

mp.dps = 50  # working precision; raise for zero-hunting (§20.5)

# Euler product sanity check (finite truncation) for Re(s) > 1
def euler_product_approx(s, num_primes=2000):
    from sympy import primerange
    from mpmath import mpf
    result = mpf(1)
    for p in primerange(2, num_primes):
        result *= 1 / (1 - mpf(p) ** (-s))
    return result

The Euler product only converges for Re(s) > 1; use it strictly as a sanity check against mpmath.zeta(s) there, never as a definition for Re(s) ≤ 1 (that requires analytic continuation, which mpmath.zeta already handles correctly).

20.4 · Known Relations (Do Not Resubmit)

20.5 · Falsification / Verification for Zero Work

1. Always verify against a published table of the first N nontrivial zeros' imaginary parts (widely available to high precision) before reporting a "new" zero location — the first few are approximately 14.134725, 21.022040, 25.010858...; a computed value that doesn't match these to the working precision indicates an implementation bug, not a discovery. Dogfood-confirmed through zero 10 (imaginary parts up to ≈49.7738), all matching the known table and all on Re(s)=1/2 — see §22 Run 3.
2. Zero-finding must report the precision used and the verification method (e.g. sign change of Z(t) the Riemann–Siegel function, or root isolation via mpmath.findroot seeded near a known approximate location) — a zero claimed without a stated precision is not verifiable and should be discarded per §7's scoring philosophy.
3. Treat any claimed zero off the critical line Re(s) = 1/2 within verified numerical range as overwhelmingly likely to be a bug (RH has been numerically verified for a very large number of zeros) — apply the same "guilty until proven robust" prior from the skill's core philosophy, and re-derive with an independent high-precision implementation before reporting anything.

§21 · Integer Partitions

21.1 · Core Objects

21.2 · Generating Function (Conceptual Ground Truth)

Σ p(n) x^n = ∏_{k=1}^∞ 1/(1 - x^k)

Restricted variants swap the product range/step: distinct parts → ∏(1+x^k); odd parts only → ∏_{k odd} 1/(1-x^k); parts ≤ m → truncate the product at k=m.

21.3 · Two Independent Implementations (§2-Style Dual Check)

A. Pentagonal number recurrence (fast, exact, no partition objects generated — analogous to the "factorization-based" implementation in §2.B):

def pentagonal_recurrence(N):
    p = [0]*(N+1)
    p[0] = 1
    for n in range(1, N+1):
        total, k = 0, 1
        while True:
            g1 = k*(3*k-1)//2
            g2 = k*(3*k+1)//2
            if g1 > n and g2 > n:
                break
            sign = 1 if k % 2 == 1 else -1
            if g1 <= n: total += sign * p[n-g1]
            if g2 <= n: total += sign * p[n-g2]
            k += 1
        p[n] = total
    return p

Uses Euler's pentagonal number theorem: p(n) = Σ_{k≠0} (-1)^{k+1} p(n - k(3k-1)/2), generalized pentagonal numbers k(3k-1)/2 for k = 1,-1,2,-2,3,-3,....

B. Direct enumeration / DP over parts (definition-based, §2.A analogue) — restricted-partition DP, e.g. for distinct-parts or bounded-part variants:

def count_partitions_odd_parts(n):
    dp = [0]*(n+1); dp[0] = 1
    for o in range(1, n+1, 2):
        for s in range(o, n+1):
            dp[s] += dp[s-o]
    return dp[n]

Cross-check A against sympy.functions.combinatorial.numbers.partition (the modern location; sympy.npartitions is deprecated) before trusting a hand-rolled recurrence at scale.

21.4 · Falsification Specific to Partitions

• p(0) = 1, not 0 — the empty partition is a valid partition of 0; this is the single most common off-by-one in partition code (mirrors the n=0 edge case rule of §4.4).
• Weak zones: pentagonal numbers themselves (1, 2, 5, 7, 12, 15, 22, 26, ...) and their n

…(truncated)