feat(bench): add Hen 2026 planted-solution factorization Ising/QUBO benchmark

Implements the integer-factorization construction of

    Hen, "Planted-solution SAT and Ising benchmarks from integer
    factorization", arXiv:2604.09837 (2026)

so QQA4CO can be exercised against problems whose ground state is known
*by construction*. The semi-prime N = p * q is encoded by long
multiplication into a circuit of AND/XOR clauses, each compiled to the
three- or four-spin energy gadgets of Eqs. 10-11 (verified directly by
the new spectrum tests). The bits of N and the leading bits of p, q are
pinned and folded into a constant offset, so the resulting QUBO is
strictly free-variable and the planted assignment evaluates to
``loss = 0``.

* `src/qqa/problems/factorization.py`
  - `IntegerFactorizationIsing(p, q)` — full QUBOProblem with
    `loss_fn`, `score_summary`, decoded-factor reporting (`p_hat`,
    `q_hat`, `N_hat`, `is_ground_state`, `decoded_N_correct`).
  - Half-adder reduction fuses the carry and XOR auxiliary spins,
    cutting variable count by ~30% relative to the naive encoding.
  - Ising -> QUBO conversion is documented in-source with the algebra
    that produced the `4 * J_sym` off-diagonal coefficient.
  - Helpers: `_is_probable_prime` (Miller-Rabin), `random_prime`,
    `random_semiprime`, `random_factorization_problems`.
* `tests/test_factorization.py` (39 tests):
  - Number-theory helpers, planted-solution invariants
    (`loss == 0`, decoder round-trip, gap >= 2),
  - **`test_and_gadget_full_spectrum`** and
    **`test_xor_gadget_full_spectrum`** check the per-clause Ising
    energies against Eqs. 10/11 byte-for-byte,
  - End-to-end `qqa.anneal` smoke test that recovers (p, q) for small N.
* `scripts/bench_factorization.py` reports both `ground_state_count`
  (E = 0) and `decoded_n_count` (factors multiply to N) so partial
  successes are visible.
* `data/factorization/README.md` documents the construction, the
  spin-count scaling table, and the BibTeX citation.
* README.md / docs/problems.md add a top-level section + the
  `@article{hen2026planted, ...}` entry under the `## Cite` block.

Author: Yuma-Ichikawa

README.md

@@ -435,6 +435,31 @@ See [`data/discs/README.md`](data/discs/README.md) for the full layout,
 optional Hugging Face source, the `qqa.datasets.discs_*` Python loaders,
 and the batched-instance (`parallel=True`) API.
 
+### Planted-solution factorization Ising/QUBO benchmark (Hen 2026)
+
+A second, complementary benchmark is the **planted-solution
+factorization Ising** of Hen, *Planted-solution SAT and Ising
+benchmarks from integer factorization* (arXiv:2604.09837, 2026). The
+construction encodes ``N = p · q`` as a circuit of AND/XOR Ising
+gadgets whose unique ground state is the bit-string of ``(p, q)``, so
+**every solver claim is verifiable** against the known optimum.
+
+```python
+import qqa
+prob = qqa.IntegerFactorizationIsing(p=11, q=13)         # N = 143
+result = qqa.anneal(prob, sol_size=500, num_epochs=3000,
+                    learning_rate=0.5,
+                    schedule=qqa.LinearBGSchedule(min_bg=-2, max_bg=0.5))
+print(result.score["extra"]["p_hat"], result.score["extra"]["q_hat"])
+```
+
+```bash
+uv run python scripts/bench_factorization.py --bits 4 --instances 5
+```
+
+See [`data/factorization/README.md`](data/factorization/README.md) for
+the construction details, scaling table, and the citation.
+
 ## Streamlit dashboard
 
 ```bash
@@ -716,3 +741,21 @@ Reference implementation: <https://github.com/google-research/discs>.
 The unified ``data/discs/`` layout (one ``.gpickle`` per instance plus a
 ``manifest.jsonl`` sidecar) is QQA4CO-specific and described in
 [`data/discs/README.md`](data/discs/README.md).
+
+If you use the **planted-solution factorization Ising benchmark**
+(`qqa.IntegerFactorizationIsing`, `scripts/bench_factorization.py`,
+`data/factorization/`), please cite the paper that introduced the
+construction and the gadget formulas this implementation follows:
+
+```bibtex
+@article{hen2026planted,
+  title   = {Planted-solution {SAT} and Ising benchmarks from integer factorization},
+  author  = {Hen, Itay},
+  journal = {arXiv preprint arXiv:2604.09837},
+  year    = {2026}
+}
+```
+
+Reference SAT/Ising compiler: <https://github.com/itay-hen/pq-SAT-benchmark>.
+The QQA4CO implementation in `src/qqa/problems/factorization.py` is
+independent and follows Eqs. 10–12 of the paper directly.

data/factorization/README.md

@@ -0,0 +1,150 @@
+# Planted-solution factorization Ising/QUBO benchmark
+
+Implementation of the **planted-solution Ising / QUBO benchmark from
+integer factorization** introduced by:
+
+> Itay Hen, *Planted-solution SAT and Ising benchmarks from integer
+> factorization*, **arXiv:2604.09837** (2026).
+
+The construction encodes the long-multiplication identity ``N = p · q``
+as a Boolean circuit of AND and XOR gates, then lowers each gate to one
+of the two Ising **energy gadgets** of the paper (Eqs. 10 and 11). The
+resulting QUBO has a **provably optimal, planted ground state** equal to
+the bits of ``(p, q)`` — making it an unusually clean benchmark for any
+discrete optimiser (SAT solvers, simulated / quantum annealers,
+gradient-based samplers, …) because every reported solution can be
+*verified* against the known minimum.
+
+Unlike random ``k``-SAT (where the ground state energy is unknown and
+solver "successes" cannot be cross-checked), every instance here has
+
+* a known optimum ``H(s*) = 0`` (after offset normalisation),
+* a single tunable difficulty knob: the bit-length ``d = max(n_p, n_q)``,
+* a strictly positive spectral gap ``≥ 2`` between the planted optimum
+  and any other configuration.
+
+## Quick start
+
+The data is generated at runtime; no download needed.
+
+```bash
+# Tiny smoke (4-bit primes, ~50 free spins; CPU, < 30 s)
+uv run python scripts/bench_factorization.py --bits 4 --instances 5
+
+# Larger sweep on GPU
+uv run python scripts/bench_factorization.py --bits 6 --instances 10 \
+    --device cuda --sol-size 1000 --num-epochs 5000 \
+    --output bench_factorization_d6.json
+```
+
+Or programmatically:
+
+```python
+import qqa
+
+# Single instance with known factors
+prob = qqa.IntegerFactorizationIsing(p=11, q=13)        # N = 143
+print(prob.num_nodes, prob.num_and, prob.num_xor)        # ~150 spins
+
+# Random benchmark suite
+suite = qqa.random_factorization_problems(
+    bit_length=5, num_instances=10, seed=0,
+)
+
+# Solve with QQA
+result = qqa.anneal(
+    prob,
+    sol_size=500, num_epochs=3000, learning_rate=0.5,
+    schedule=qqa.LinearBGSchedule(min_bg=-2, max_bg=0.5),
+    curve_rate=4, div_param=0.2,
+)
+print(result.score["extra"]["p_hat"], result.score["extra"]["q_hat"])
+```
+
+`prob.score_summary(x)` reports the residual energy above the planted
+optimum, the **decoded** factor pair ``(p̂, q̂)``, ``N̂ = p̂ · q̂``, and
+the bit-Hamming distance to ``s*``. ``score['feasible']`` is `True` iff
+the optimiser reached the planted ground state and ``p̂ · q̂ = N``.
+
+## What the construction looks like inside
+
+| Pipeline step                                   | Where in the code                          |
+|-------------------------------------------------|--------------------------------------------|
+| Random equal-bit-length prime sampling          | `random_prime`, `random_semiprime`         |
+| Long-multiplication → AND/XOR clause graph      | `_compile`                                 |
+| AND gadget (Eq. 10), XOR gadget (Eq. 11)        | `_add_and_gadget`, `_add_xor_gadget`       |
+| Pin substitution + Ising → QUBO                 | `_build_qubo_from_compilation`             |
+| Public ``qqa.IntegerFactorizationIsing``        | top of `src/qqa/problems/factorization.py` |
+| End-to-end smoke + planted-invariant tests      | `tests/test_factorization.py`              |
+
+The implementation is **gadget-faithful** to the paper but does *not*
+ship the preprocessing pipeline of Sec. III. Problem sizes therefore
+scale as ``O(d^4)`` in the symmetric case ``n_p = n_q = d``:
+
+| ``d`` | free spins (typical) | suitable for |
+|-------|----------------------|--------------|
+| 2–3   |  16 – 40             | CPU smoke    |
+| 4–5   | 100 – 250            | CPU bench    |
+| 6–7   | 500 – 1 500          | single-GPU   |
+| ≥ 8   | ≥ 3 000              | multi-GPU / preprocessing |
+
+Concrete sizes (``carry`` and XOR-aux fused per Eq. 12 footnote — see
+``_compile`` in `src/qqa/problems/factorization.py`):
+
+| ``N``  | ``(p, q)`` | free spins | AND clauses | XOR clauses |
+|--------|-----------|------------|-------------|-------------|
+| 15     | (3, 5)    | 16         | 12          | 6           |
+| 35     | (5, 7)    | 37         | 24          | 15          |
+| 143    | (11, 13)  | 102        | 60          | 44          |
+| 221    | (13, 17)  | 148        | 85          | 65          |
+| 667    | (23, 29)  | 213        | 120         | 95          |
+
+For larger ``d`` the upstream paper's preprocessing is essential; an
+optional pin-propagation / AND-XOR simplification pass is a natural
+follow-up that would slot directly into ``_build_qubo_from_compilation``.
+
+## What QQA can do today
+
+The construction is **correct for any ``(p, q)``**: ``planted_x``
+evaluates to exactly ``loss = 0``, any single bit-flip raises the
+loss by at least 2 (the gadget gap; see
+`tests/test_factorization.py::test_planted_is_strict_minimum`), and the
+gadget energy tables match Eqs. (10)–(11) of the paper element-wise
+(`test_and_gadget_full_spectrum`, `test_xor_gadget_full_spectrum`).
+
+The bench script reports **two** orthogonal success metrics:
+
+* `gs` (ground state) — `loss == 0`, i.e. the entire spin string sits
+  on the planted optimum. This is the *strong* contract.
+* `decoded_N` — the bits on the input wires `p_0..p_{n_p-1}` and
+  `q_0..q_{n_q-1}` decode to a pair `(p̂, q̂)` with `p̂ · q̂ = N`,
+  even if internal pp/sum/carry spins are still inconsistent (so
+  `loss > 0`). This is the *weak* contract: the optimiser found the
+  factorisation but did not finish the auxiliary bookkeeping.
+
+Empirically, the default ``qqa.anneal`` recipe consistently reaches
+`gs` for ``d ≤ 3`` (``N ≤ 7·7``) and gets `decoded_N` for ``d = 4``
+(``N ≈ 11·13``) within a few seconds on CPU; ``d ≥ 5`` benefits from
+larger ``sol_size`` and ``num_epochs`` and may still leave a non-zero
+gap, exactly as the paper predicts (``T_med`` grows as ``2^{β d}``
+with ``β ≈ 1``). This is the *whole point* of the benchmark: it is
+genuinely difficult for general-purpose CO solvers even at moderate
+``d``.
+
+## Citation
+
+Please cite the original paper:
+
+```bibtex
+@article{hen2026planted,
+  title   = {Planted-solution {SAT} and Ising benchmarks from integer factorization},
+  author  = {Hen, Itay},
+  journal = {arXiv preprint arXiv:2604.09837},
+  year    = {2026},
+}
+```
+
+The reference SAT/Ising compiler is open source at
+<https://github.com/itay-hen/pq-SAT-benchmark>; the QQA4CO
+implementation here is independent and follows the published gadget
+formulas (Eqs. 10–12) directly.

docs/problems.md

@@ -38,6 +38,28 @@ The batched path is supported by ``qqa.anneal`` only — SA / PA backends
 require a single dense ``Q_mat`` and intentionally reject
 batched-instance problems.
 
+## Planted-solution factorization Ising/QUBO (Hen 2026)
+
+| Class                                   | Inputs            | Loss / property                                  |
+| --------------------------------------- | ----------------- | ------------------------------------------------ |
+| `IntegerFactorizationIsing(p, q)`       | two primes        | $x^T Q x + E_0$, **min = 0** on the planted bits |
+| `random_factorization_problems(d, k)`   | bit-length, count | list of $k$ random instances of width $d$        |
+
+Implements arXiv:2604.09837. The semi-prime $N = p \cdot q$ is encoded
+via long multiplication into a circuit of AND and XOR clauses; each
+clause becomes the three- or four-spin energy gadget of Eqs. 10–11 of
+the paper. Pinned spins (the bits of $N$ and the leading bits of
+$p$, $q$) are folded into the constant offset, so the resulting QUBO is
+strictly free-variable.
+
+`score_summary(x)` returns the residual energy above the planted optimum,
+the **decoded** factor pair `(p̂, q̂)`, the product `N̂`, and the
+bit-Hamming distance to the planted assignment. `feasible == True`
+exactly when the optimiser found the planted solution and `p̂ · q̂ = N`.
+
+See `data/factorization/README.md` for the construction details, problem
+sizes vs. bit-length, and the citation.
+
 ## Classic CO (binary, problem-specific losses)
 
 | Class | Inputs | Loss / feasibility |

scripts/bench_factorization.py

@@ -0,0 +1,137 @@
+"""Benchmark runner for the planted-solution factorization Ising/QUBO suite.
+
+Usage examples
+--------------
+$ uv run python scripts/bench_factorization.py --bits 4 --instances 5
+$ uv run python scripts/bench_factorization.py --bits 5 --instances 3 \
+      --device cuda --sol-size 500 --num-epochs 5000
+
+For each instance the script reports:
+
+* the planted ``(p, q)`` and ``N``;
+* the number of free Ising spins (problem size after pin folding);
+* QQA's residual energy ``H(x̂) − E_0`` (zero if the planted is found);
+* the decoded ``(p̂, q̂)``, ``N̂`` and the bit-Hamming distance to ``x*``;
+* wall-clock time.
+
+The benchmark is intentionally **CPU-friendly** for ``bits ≤ 5`` so that
+contributors can iterate without a GPU; for ``bits ≥ 6`` the recommended
+defaults are ``--device cuda --sol-size 1000 --num-epochs 5000``.
+
+Citation
+--------
+Hen, *Planted-solution SAT and Ising benchmarks from integer
+factorization*, arXiv:2604.09837 (2026).
+"""
+
+from __future__ import annotations
+
+import argparse
+import sys
+import time
+from pathlib import Path
+from typing import Any
+
+import qqa
+
+# `scripts/` is on sys.path when this file is invoked as
+# ``python scripts/bench_factorization.py``; importing as a sibling module
+# works without packaging.
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+import _bench_common as bench  # noqa: E402
+
+
+def main(argv: list[str] | None = None) -> int:
+    p = argparse.ArgumentParser(description=__doc__)
+    p.add_argument("--bits", type=int, default=4,
+                   help="bit-length d of each prime factor (≥ 2). Problem size "
+                        "scales as O(d^4) without preprocessing.")
+    p.add_argument("--instances", type=int, default=5,
+                   help="number of random semiprimes to benchmark.")
+    p.add_argument("--seed", type=int, default=0)
+    p.add_argument("--device", default="auto", help="cpu / cuda / auto")
+    p.add_argument("--sol-size", type=int, default=200)
+    p.add_argument("--num-epochs", type=int, default=2000)
+    bench.add_qqa_hp_args(p)
+    p.add_argument("--output", type=Path, default=None,
+                   help="optional JSON file to dump per-instance results.")
+    p.add_argument("-v", "--verbose", action="store_true")
+    args = p.parse_args(argv)
+
+    log = bench.setup_logging(args.verbose, name="bench_factorization")
+    bench.setup_device(args)
+    log.info("bits=%d instances=%d device=%s", args.bits, args.instances, args.device)
+    hp = bench.qqa_hp_kwargs(args)
+    log.info(
+        "qqa hp: lr=%g temp=%g curve_rate=%d gamma=[%g,%g] div=%g sol_size=%d num_epochs=%d",
+        hp["learning_rate"], hp["temp"], hp["curve_rate"],
+        hp["gamma_min"], hp["gamma_max"], hp["div_param"],
+        args.sol_size, args.num_epochs,
+    )
+
+    suite = qqa.random_factorization_problems(
+        bit_length=args.bits, num_instances=args.instances, seed=args.seed,
+        device=args.device,
+    )
+
+    rows: list[dict[str, Any]] = []
+    n_ground_state = 0  # E_min = 0 AND decoded N̂ = N
+    n_decoded = 0       # decoded N̂ = N (input bits right, internal bits maybe not)
+    for k, prob in enumerate(suite):
+        t0 = time.time()
+        result = bench.run_qqa_anneal(
+            prob, device=args.device, sol_size=args.sol_size,
+            num_epochs=args.num_epochs, **hp,
+        )
+        wall = time.time() - t0
+        s = result.score
+        ex = s["extra"]
+        # Two distinct success criteria:
+        # - `decoded`: the optimiser placed the right bits on the input wires
+        #   p_0..p_{n_p-1} and q_0..q_{n_q-1} (cheap to satisfy locally).
+        # - `gs`: the entire spin string sits at the planted ground state
+        #   (energy 0), which additionally requires every internal pp/sum/
+        #   carry spin to be consistent with the multiplication circuit.
+        decoded_ok = bool(ex["matches_planted"])
+        gs_ok = bool(s["feasible"])
+        n_decoded += int(decoded_ok)
+        n_ground_state += int(gs_ok)
+        log.info(
+            "  [%d/%d] N=%d=%d*%d  free=%d  E_min=%.3f  decoded=%d*%d=%d  "
+            "decoded_N=%s  gs=%s  t=%.2fs",
+            k + 1, args.instances, ex["N"], ex["p"], ex["q"], ex["num_free_spins"],
+            float(s["value"]), ex["p_hat"], ex["q_hat"], ex["N_hat"],
+            "Y" if decoded_ok else "N", "Y" if gs_ok else "N", wall,
+        )
+        rows.append({
+            "k": k, "N": ex["N"], "p": ex["p"], "q": ex["q"],
+            "num_free_spins": ex["num_free_spins"],
+            "energy_above_planted": float(s["value"]),
+            "p_hat": ex["p_hat"], "q_hat": ex["q_hat"], "N_hat": ex["N_hat"],
+            "decoded_N_correct": decoded_ok,
+            "is_ground_state": gs_ok,
+            "hamming_to_planted": ex["hamming_to_planted"],
+            "wall_s": wall,
+        })
+
+    log.info(
+        "== SUMMARY == ground-state %d/%d (%.1f%%)  decoded-N %d/%d (%.1f%%)",
+        n_ground_state, len(suite), 100.0 * n_ground_state / max(1, len(suite)),
+        n_decoded, len(suite), 100.0 * n_decoded / max(1, len(suite)),
+    )
+    if args.output:
+        payload = {
+            "bits": args.bits, "instances": args.instances, "device": args.device,
+            "qqa_hp": {**hp, "sol_size": args.sol_size, "num_epochs": args.num_epochs},
+            "ground_state_count": n_ground_state,
+            "decoded_n_count": n_decoded,
+            "total": len(suite),
+            "results": rows,
+        }
+        bench.dump_results_json(args.output, payload)
+        log.info("wrote %s", args.output)
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())