NitroSAT: Solving MaxSAT Through Number Theory, Statistical Mechanics, and Persistent Homology
Sethurathienam Iyer
https://orcid.org/0009-0008-5446-2856
https://zenodo.org/records/18214172
2026-02-17
1 Introduction
The Boolean Satisfiability problem (SAT) and its optimization variant MaxSAT occupy a central position in computational complexity theory. As the canonical NP-complete problem, SAT underpins verification, planning, scheduling, and combinatorial optimisation across both theoretical computer science and industrial applications. Complete MaxSAT solvers — those that guarantee optimality — typically rely on branch-and-bound with clause-learning (e.g., RC2, Ignatiev, Morgado, and Marques-Silva 2019; Morgado et al. 2013), achieving strong performance on moderate-sized instances but encountering prohibitive runtimes on large structured problems where symmetry and local minima dominate the search landscape.
This article presents NitroSAT, a continuous relaxation solver for MaxSAT that draws its mathematical machinery from four disparate fields: number theory [Riemann zeta function, prime-weighted clause structures; Titchmarsh (1986); Montgomery and Vaughan (2007)], statistical mechanics [partition function estimation, fracture detection via heat capacity; Mézard and Montanari (2009); Kirkpatrick, Gelatt Jr., and Vecchi (1983)], spectral geometry [heat kernel smoothing on the clause-variable hypergraph; Chung (1997); Grigor’yan (2009)], and algebraic topology [persistent homology for cycle detection in the unsatisfied constraint graph; Edelsbrunner and Harer (2008); Carlsson (2009)]. The solver is implemented in approximately 2,500 lines of standalone LuaJIT, achieving 99–100% clause satisfaction on hard structured instances — including a 650,000-clause graph coloring problem solved in 4.6 seconds and a 354,890-clause clique coloring instance solved perfectly in 3.5 seconds.
NitroSAT is not a complete solver; it cannot certify unsatisfiability. Rather, it operates as a high-quality MaxSAT approximation engine with structural UNSAT awareness: the ability to detect, via thermodynamic phase transitions, when an instance is likely unsatisfiable, and to terminate early rather than cycling indefinitely. This property distinguishes it from standard stochastic local search methods (WalkSAT, Selman, Kautz, and Cohen 1996; probSAT, Hoos 2002) that lack termination criteria on unsatisfiable instances.
The remainder of this article is organised as follows. Section 2 develops the mathematical foundations. Section 3 describes the solver architecture and its multi-phase solution protocol. Section 4 provides an annotated code walkthrough of the key algorithmic components. Section 5 presents benchmark results on structured satisfiable instances, scheduling problems, CNFgen benchmarks, and XOR-SAT stress tests. Section 6 discusses limitations and future directions.
2 Mathematical Foundations
2.1 The Partition Function Factorisation
NitroSAT’s energy landscape is derived from a partition function that factorises into an arithmetic component (encoding clause structure via number theory) and a spectral component (encoding variable connectivity via the constraint hypergraph Laplacian):
\[ Z(\beta) = \zeta(\beta) \cdot \operatorname{Tr}\!\left(e^{-\beta L}\right) \]
where \(\beta\) is the inverse temperature parameter, \(\zeta(\beta)\) denotes the Riemann zeta function evaluated along the real line, and \(L\) is the graph Laplacian of the clause-variable interaction graph.
The logarithmic derivative of this partition function yields the gradient used for continuous optimisation:
\[ \nabla \log Z(\beta) = \nabla \log \zeta(\beta) + \nabla \log \operatorname{Tr}\!\left(e^{-\beta L}\right) \]
The first term introduces arithmetic structure through the von Mangoldt function \(\Lambda(n)\) (Titchmarsh 1986), while the second encodes spectral geometry through the heat kernel \(e^{-\beta L}\) (Grigor’yan 2009; Kondor and Lafferty 2002). The key insight is that at low temperatures (\(\beta \to \infty\)), the zeta component drives symmetry-breaking among degenerate solutions, while the spectral component smooths the energy landscape to facilitate gradient descent.
2.2 Prime-Weighted Clause Energies
Each clause \(c\) in the MaxSAT instance is assigned a weight derived from the \(c\)-th prime number \(p_c\):
\[ w_c = \left(1 + \log p_c\right)^{-\alpha} \]
where \(\alpha \geq 0\) is a tuneable parameter (default \(\alpha = 1\)). This weighting scheme has a number-theoretic motivation: the primes provide an incommensurable weighting basis that prevents resonance between clause weights, thereby reducing the likelihood of the optimiser becoming trapped in symmetric local minima. The logarithmic scaling ensures that the weight distribution has bounded variance even as the clause count grows, since by the Prime Number Theorem \(p_c \sim c \log c\).
# Generate primes via sieve
sieve_primes <- function(n) {
limit <- ceiling(n * log(n) * 1.3) + 100
is_prime <- rep(TRUE, limit)
is_prime[1] <- FALSE
for (i in 2:floor(sqrt(limit))) {
if (is_prime[i]) is_prime[seq(i*i, limit, i)] <- FALSE
}
which(is_prime)[1:n]
}
primes <- sieve_primes(500)
alpha <- 1.0
weights <- (1 + log(primes))^(-alpha)
df_weights <- data.frame(
clause_index = 1:500,
prime = primes,
weight = weights
)
ggplot(df_weights, aes(x = clause_index, y = weight)) +
geom_line(color = "#2E75B6", linewidth = 0.6) +
geom_point(data = df_weights[seq(1, 500, 25), ],
aes(x = clause_index, y = weight),
color = "#C0392B", size = 1.5) +
labs(
title = expression("Prime-Weighted Clause Energies: " * w[c] == (1 + log~p[c])^{-alpha}),
subtitle = expression("Incommensurable weights from the prime number theorem (" * alpha * " = 1.0)"),
x = "Clause Index (c)",
y = expression("Weight " * w[c])
) +
theme_nitrosat()
Prime-weighted clause energies for the first 500 clauses. The logarithmic prime weighting provides an incommensurable basis that prevents symmetric trapping.
The number of prime harmonics used for the zeta perturbation is itself adaptive, scaling as:
\[ K = \left\lfloor 3 \cdot \ln(n) \cdot \ln(\alpha_{\text{ratio}} + 1) \right\rfloor, \quad \alpha_{\text{ratio}} = \frac{m}{n} \]
where \(n\) and \(m\) are the variable and clause counts respectively. This formula couples symmetry resolution (via \(\ln n\)) with landscape ruggedness (via the clause-to-variable ratio).
2.3 Heat Kernel Smoothing on the Constraint Hypergraph
The clause-variable interaction structure defines a hypergraph whose Laplacian \(L\) governs diffusion dynamics. NitroSAT applies heat kernel smoothing to the gradient field: each variable \(i\) receives a multiplier
\[ h_i = 1 + \lambda \cdot \exp\!\left(-\beta \cdot d_i\right) \]
where \(d_i\) is the degree of variable \(i\) in the clause-variable graph (the number of literal occurrences across all clauses), \(\beta\) is the heat kernel bandwidth, and \(\lambda\) is the smoothing intensity. This multiplier amplifies gradients for low-degree (peripheral) variables and dampens gradients for high-degree (hub) variables, implementing an implicit diffusion that smooths local minima in the energy landscape.
degrees <- seq(0, 100, 1)
betas <- c(0.1, 0.3, 0.5, 1.0)
lambda <- 0.1
df_heat <- expand.grid(degree = degrees, beta = betas) %>%
mutate(
multiplier = 1 + lambda * exp(-beta * degree),
beta_label = paste0("β = ", beta)
)
ggplot(df_heat, aes(x = degree, y = multiplier, color = beta_label)) +
geom_line(linewidth = 0.8) +
scale_color_brewer(palette = "Set1", name = "Bandwidth") +
labs(
title = "Heat Kernel Gradient Multiplier",
subtitle = expression("h"[i] == 1 + lambda %.% exp(-beta %.% d[i]) * ", " * lambda * " = 0.1"),
x = expression("Variable Degree " * d[i]),
y = expression("Multiplier " * h[i])
) +
theme_nitrosat()
Heat kernel multiplier as a function of variable degree for various bandwidth parameters. Lower bandwidth (larger β) concentrates smoothing on peripheral variables.
The spectral initialisation step goes further: before the main optimisation loop, NitroSAT computes the leading eigenvector of the XOR-Laplacian via power iteration, then maps this spectral embedding to initial variable assignments. This pre-conditioning aligns the starting point with the principal structure of the constraint hypergraph.
2.4 Persistent Homology for Topological Awareness
A distinctive feature of NitroSAT is its use of persistent homology (Edelsbrunner, Letscher, and Zomorodian 2002; Zomorodian and Carlsson 2005) to detect topological structure in the unsatisfied constraint graph. At periodic intervals during optimisation, the solver constructs a similarity graph over variables: two variables are connected if they co-occur in unsatisfied clauses, with edge weight proportional to co-occurrence frequency normalised by degree. A filtration is applied by thresholding on edge similarity, and the Betti numbers \(\beta_0\) (connected components) and \(\beta_1\) (independent cycles) are computed.
The zeroth Betti number \(\beta_0\) counts the number of disconnected clusters of “stuck” variables. The first Betti number \(\beta_1\) counts independent cycles in the unsatisfied constraint graph — these correspond to circular dependency structures where no single variable flip can resolve the conflict. NitroSAT tracks the persistence of these topological features across optimisation steps, identifying long-lived cycles as structural bottlenecks requiring coordinated multi-variable repair.
The complexity score \(\kappa = \beta_1 / \max(\beta_0, 1)\) provides a scalar measure of constraint entanglement. When \(\kappa\) exceeds a threshold, NitroSAT activates its topological repair phase (Section 3.3), which identifies and targets cycle-breaking variables for coordinated perturbation.
2.5 Fracture Detection via Statistical Mechanics
NitroSAT models the satisfiability landscape as a statistical mechanical system at inverse temperature \(\beta\) (Mézard, Parisi, and Zecchina 2002; Monasson et al. 1999; Krzakała et al. 2007). The partition function is estimated by importance sampling:
\[ \hat{Z}(\beta) = \frac{1}{N_s} \sum_{i=1}^{N_s} e^{-\beta E(\mathbf{x}_i)} \]
where \(E(\mathbf{x})\) counts the number of unsatisfied clauses for assignment \(\mathbf{x}\), and the samples \(\mathbf{x}_i\) are drawn uniformly at random (respecting decimated variables). From \(\hat{Z}(\beta)\), the solver estimates the mean energy \(\langle E \rangle\), the energy variance \(\operatorname{Var}(E)\) (proportional to heat capacity \(C_v = \beta^2 \operatorname{Var}(E)\)), and the third central moment (skewness, indicating asymmetry near branch splitting).
A fracture is declared when two independent thermodynamic signals coincide:
- Slope shock: The rate of change \(|\partial \log Z / \partial \beta|\) exceeds the running mean by more than \(\sigma_f\) standard deviations (default \(\sigma_f = 1.5\)).
- Variance spike: The energy variance exceeds a \(z\)-score threshold relative to its running statistics.
This dual-signal requirement guards against false positives from noise in the sampled partition function. Upon detecting a fracture, NitroSAT uses the Lambert \(W\) function to enumerate alternative solution branches predicted by the free energy landscape, scores them by sampling, and executes a coordinated jump to the most promising branch.
set.seed(42)
n_steps <- 200
beta_seq <- seq(0.5, 10, length.out = n_steps)
# Simulate log Z with a phase transition around beta = 5
logZ <- -0.5 * beta_seq + 2 * exp(-0.5 * (beta_seq - 5)^2) + rnorm(n_steps, 0, 0.15)
slope <- c(0, diff(logZ) / diff(beta_seq))
# Simulate variance with spike at transition
variance <- 0.5 + 3 * exp(-0.8 * (beta_seq - 5.2)^2) + abs(rnorm(n_steps, 0, 0.1))
df_fracture <- data.frame(
beta = beta_seq,
slope = abs(slope),
variance = variance
)
# Detect fracture region
df_fracture <- df_fracture %>%
mutate(
slope_z = (slope - mean(slope)) / sd(slope),
var_z = (variance - mean(variance)) / sd(variance),
fracture = slope_z > 1.5 & var_z > 0.8
)
p1 <- ggplot(df_fracture, aes(x = beta)) +
geom_line(aes(y = slope_z, color = "Slope Shock (z-score)"), linewidth = 0.7) +
geom_line(aes(y = var_z, color = "Variance Spike (z-score)"), linewidth = 0.7) +
geom_hline(yintercept = 1.5, linetype = "dashed", color = "grey50", linewidth = 0.4) +
geom_rect(
data = df_fracture %>% filter(fracture) %>%
summarise(xmin = min(beta), xmax = max(beta)),
aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
fill = "red", alpha = 0.1, inherit.aes = FALSE
) +
annotate("text", x = 5.5, y = 4, label = "FRACTURE\nDETECTED",
color = "red", fontface = "bold", size = 4) +
scale_color_manual(values = c("Slope Shock (z-score)" = "#2E75B6",
"Variance Spike (z-score)" = "#C0392B"),
name = "") +
labs(
title = "Dual-Signal Fracture Detection",
subtitle = "Phase transition identified by simultaneous slope shock and variance spike",
x = expression("Inverse Temperature " * beta),
y = "Standardised Signal (z-score)"
) +
theme_nitrosat()
print(p1)
Schematic of the dual-signal fracture detection. A fracture is declared only when both the slope shock (rate of change of log Z) and the variance spike (heat capacity proxy) exceed their respective thresholds simultaneously.
2.6 Branch Jumping via the Lambert W Function
When a fracture is detected at critical inverse temperature \(\beta_c\), NitroSAT uses the Lambert \(W\) function (Corless et al. 1996) to enumerate alternative solution branches. The key equation arises from the free energy landscape: defining \(u = \beta - \beta_c\), the branch structure satisfies \(u \cdot e^u = \xi\) for some parameter \(\xi\) determined by the current state. The Lambert \(W\) function inverts \(w \cdot e^w = z\), yielding:
\[ \beta_0 = \beta_c + W_0(\xi), \quad \beta_{-1} = \beta_c + W_{-1}(\xi) \]
The principal branch \(W_0\) returns to a nearby basin, while the \(W_{-1}\) branch (defined for \(\xi \in [-1/e, 0)\)) identifies a distant alternative basin in the energy landscape. Each candidate branch is scored by sampling and evaluated via Metropolis–Hastings acceptance, with the optimizer state (momentum, velocity) reset upon a successful jump to avoid carrying stale gradient information into the new basin.
3 Solver Architecture
3.1 Overview
NitroSAT’s solution protocol operates in multiple phases, each activating when the previous phase plateaus:
NitroSAT’s multi-phase solution protocol. Each phase activates when the previous phase’s satisfaction rate crosses the indicated threshold.
Phase 1 (Gradient Descent + BAHA). Continuous relaxation with Nadam optimisation (Kingma and Ba 2015; Dozat 2016), prime-weighted log-barrier energy, heat kernel smoothing, zeta zero perturbation, and BAHA fracture detection with branch jumping (Iyer 2026a; Wales and Doye 1997). This phase handles the bulk of the optimisation, typically reaching 95–99% satisfaction.
Zeta Sweep Finisher. A stochastic repair pass using prime-weighted “Adelic forces” to flip variables in unsatisfied clauses with probability proportional to a sigmoid of the accumulated force intensity.
Phase 2 (Topological Repair). Activated at \(\geq 95\%\) satisfaction. Computes persistent homology to identify cycle-breaking variables, then applies focused perturbation and restricted WalkSAT on the topological bottleneck region.
Phase 3 (Adelic Saturation). Activated at \(\geq 98\%\) satisfaction. Binary search on inverse temperature \(\beta\) with aggressive zeta-weighted forces combined with unit propagation on decimated variables.
Core Decomposition. Identifies the unsatisfied core (remaining unsatisfied clauses), amplifies their weights by \(100\times\), and runs restricted WalkSAT on the core variables.
BAHA-WalkSAT. A phase-transition-aware discrete local search that combines Metropolis–Hastings acceptance with BAHA fracture detection. Upon detecting a phase transition in the discrete search, it executes coordinated multi-variable flips on the most unstable variables, with jump size scaled by the variance \(z\)-score.
3.2 The Continuous Relaxation
Each Boolean variable \(x_i \in \{0, 1\}\) is relaxed to a continuous variable \(x_i \in [0, 1]\). The energy function for a clause \(c\) with literals \(\ell_1, \ldots, \ell_k\) is the product of literal violations:
\[ V_c(\mathbf{x}) = \prod_{j=1}^{k} \bar{\ell}_j(\mathbf{x}), \quad \text{where } \bar{\ell}_j = \begin{cases} 1 - x_i & \text{if } \ell_j = x_i \\ x_i & \text{if } \ell_j = \neg x_i \end{cases} \]
The log-barrier objective is then:
\[ \mathcal{L}(\mathbf{x}) = -\sum_{c=1}^{m} w_c \cdot \log\!\left(\varepsilon + 1 - V_c(\mathbf{x})\right) + \eta \cdot H(\mathbf{x}) \]
where \(H(\mathbf{x}) = -\sum_i \left[x_i \log x_i + (1-x_i)\log(1-x_i)\right]\) is the binary entropy regulariser (encouraging exploration away from boundary values), and \(\eta\) is the entropy weight. The gradient with respect to variable \(x_i\) accumulates contributions from all clauses containing \(x_i\), weighted by the prime-weighted clause energy and multiplied by the heat kernel smoothing factor.
3.3 Oscillating Annealing Schedule
The learning rate follows a multi-frequency damped oscillation inspired by physical phase transitions across different energy scales:
\[ \text{lr}(t) = A_0 \cdot \left[\sin\!\left(\frac{\pi t}{20}\right)e^{-\pi t/20} + \sin\!\left(\frac{\pi t}{10}\right)e^{-\pi t/10} + \sin\!\left(\frac{\pi t}{9}\right)e^{-\pi t/9} + \sin\!\left(\frac{t}{\varphi^2}\right)e^{-t/\varphi^2}\right] \]
where \(\varphi = (1+\sqrt{5})/2\) is the golden ratio. Each term represents a “phase of matter” (solid, liquid, gas, golden), with the irrational frequency \(1/\varphi^2\) ensuring quasi-periodic coverage of the parameter space without exact repetition.
phi <- (1 + sqrt(5)) / 2
t_seq <- seq(0, 30, length.out = 1000)
A0 <- 0.1
solid <- sin(pi/20 * t_seq) * exp(-pi/20 * t_seq)
liquid <- sin(pi/10 * t_seq) * exp(-pi/10 * t_seq)
gas <- sin(pi/9 * t_seq) * exp(-pi/9 * t_seq)
golden <- sin(1/(phi^2) * t_seq) * exp(-1/(phi^2) * t_seq)
total <- A0 * (solid + liquid + gas + golden)
df_anneal <- data.frame(
t = rep(t_seq, 5),
lr = c(solid, liquid, gas, golden, total),
component = rep(c("Solid (π/20)", "Liquid (π/10)", "Gas (π/9)",
paste0("Golden (1/φ²)"), "Total lr(t)"), each = length(t_seq))
)
ggplot(df_anneal, aes(x = t, y = lr, color = component)) +
geom_line(aes(linewidth = component)) +
scale_linewidth_manual(values = c(0.4, 0.4, 0.4, 0.4, 1.0), guide = "none") +
scale_color_manual(values = c("#95A5A6", "#95A5A6", "#95A5A6", "#95A5A6", "#C0392B"),
name = "") +
labs(
title = "Oscillating Annealing Schedule",
subtitle = "Multi-frequency damped oscillation with golden ratio quasi-periodicity",
x = "Rescaled Time (t)",
y = "Learning Rate Intensity"
) +
theme_nitrosat()
The oscillating annealing schedule. Four damped sinusoidal components with incommensurable frequencies produce a quasi-periodic learning rate that avoids exact repetition.
4 Annotated Code Walkthrough
This section provides an annotated walkthrough of the key algorithmic
components in the NitroSAT source (nitrosat.lua). The full
implementation is approximately 2,500 lines of LuaJIT.
4.1 Data Structures and Memory Layout
NitroSAT uses LuaJIT’s Foreign Function Interface (FFI) for all performance-critical arrays, avoiding Lua’s garbage-collected heap and enabling the solver to handle instances with millions of clauses:
-- FFI Helpers for massive scale
local function ffi_double(n) return ffi.new("double[?]", n + 1) end
local function ffi_int(n) return ffi.new("int32_t[?]", n + 1) end
local function ffi_uint8(n) return ffi.new("uint8_t[?]", n + 1) endClauses are stored in a flat Compressed Sparse Row (CSR) format. The
array clauses_flat holds all literals contiguously, and
clauses_offsets[c] gives the start index of clause \(c\).
This layout achieves cache-friendly
sequential access during gradient computation and satisfaction
checking:
self.clauses_flat = ffi_int(total_literals)
self.clauses_offsets = ffi_int(self.num_clauses + 1)
local curr_idx = 0
for c = 1, self.num_clauses do
self.clauses_offsets[c-1] = curr_idx
local clause = instance.clauses[c]
for i = 1, #clause do
self.clauses_flat[curr_idx] = clause[i]
curr_idx = curr_idx + 1
end
end
self.clauses_offsets[self.num_clauses] = curr_idxA variable-to-clause index (v2c_ptr,
v2c_data) is also stored in CSR format, enabling O(1)
lookup of all clauses containing a given variable — critical for
incremental WalkSAT operations.
4.2 The Log-Barrier Gradient
The gradient computation is the inner loop of Phase 1. For each clause, the product of literal violations \(V_c\) is computed, and the log-barrier gradient is accumulated:
function NitroSat:compute_gradients()
local grads = ffi_double(self.num_vars)
for i = 1, self.num_vars do grads[i] = 0.0 end
local unsat_count = 0
for c = 0, self.num_clauses - 1 do
local violation = 1.0
local is_discrete_sat = false
for i = start, stop - 1 do
local lit = self.clauses_flat[i]
local var = math.abs(lit)
local val = self.x[var]
-- Check discrete satisfaction
if (lit > 0 and val > 0.5) or (lit < 0 and val <= 0.5) then
is_discrete_sat = true
end
-- Accumulate continuous violation product
local lit_val = lit > 0 and (1.0 - val) or val
violation = violation * lit_val
end
-- Log-barrier gradient: -w_c * ∇log(ε + 1 - V_c)
local slack = 1.0 - violation
local barrier = 1.0 / (1e-6 + slack)
local coef = self.clause_weights[c]
if self.is_hard_constraint[c] ~= 0 then coef = coef * 50.0 end
for i = start, stop - 1 do
local lit = self.clauses_flat[i]
local var = math.abs(lit)
local lit_val = ...
local term_grad = (violation / lit_val) * (lit > 0 and -1 or 1)
grads[var] = grads[var] + coef * barrier * term_grad
end
end
-- Heat kernel smoothing + entropy regularisation
local multipliers = self:compute_heat_multipliers()
for i = 1, self.num_vars do
grads[i] = grads[i] * multipliers[i]
local val = clamp(self.x[i], 1e-9, 1.0 - 1e-9)
grads[i] = grads[i] + self.entropy_weight * math.log((1.0 - val) / val)
end
return grads, unsat_count
endThe entropy gradient \(\partial H / \partial x_i = \log\!\left((1-x_i)/x_i\right)\) acts as a regulariser that pushes variables away from 0 and 1 during early optimisation, encouraging exploration before the annealing schedule drives convergence.
4.3 Zeta Zero Perturbation
The zeta perturbation applies number-theoretic corrections to the variable assignments. Each variable \(i\) receives a perturbation derived from prime harmonics along the Riemann critical line \(\operatorname{Re}(s) = 1/2\):
local function zeta_zero_perturbation(var_idx, step, max_steps, primes, prime_count)
local progress = step / max_steps
local critical_line = 0.5 -- Re(s) = 1/2
local zero_guidance = 0
for i = 0, prime_count - 1 do
local p = primes[i]
if p then
local zeta_map = math.sin(2 * math.pi * (p / var_idx)) * progress
zero_guidance = zero_guidance + zeta_map * (1 / math.sqrt(i + 1))
end
end
local critical_osc = math.sin(2 * math.pi * critical_line * progress) * 0.1
return zero_guidance + critical_osc
endThe perturbation amplitude scales with optimisation progress (stronger near convergence) and decays as \(1/\sqrt{k+1}\) across the prime harmonic series, ensuring convergence of the total perturbation while maintaining sufficient symmetry-breaking power.
4.4 BAHA-WalkSAT: Phase-Transition-Aware Local Search
The BAHA-WalkSAT finisher combines Metropolis–Hastings acceptance with BAHA fracture detection in the discrete search space. The key innovation is the coordinated multi-flip upon fracture detection:
-- Upon fracture detection:
-- 1. Rank variables by "instability" (count of unsatisfied clause appearances)
-- 2. Scale jump size by variance z-score
-- 3. Enumerate Lambert W branches around critical beta
-- 4. Execute coordinated multi-flip with M-H acceptance/revert
local base_K = floor(sqrt(#unsat_list) + 0.5)
local jump_scale = 1.0 + self.bw_var_jump_gain * max(0.0, var_z)
local K = floor(base_K * jump_scale + 0.5)
-- Flip top-K most unstable variables simultaneously
for i = 1, K do
self:execute_flip(ranked[i].var, clause_sat_counts, unsat_list, unsat_pos)
end
-- Accept or revert via Metropolis-Hastings at branch beta
if energy_after > energy_before then
if rng() >= exp(-branch_beta * (energy_after - energy_before)) then
-- Revert: restore saved state
...
end
endThe jump size \(K = \lceil \sqrt{|\text{unsat}|} \cdot (1 + \gamma \cdot z_{\text{var}}) \rceil\) adapts to both the current unsatisfied count and the strength of the detected phase transition, enabling aggressive exploration near critical points while remaining conservative in smooth regions.
5 Benchmark Results
5.1 Structured Satisfiable Instances
benchmarks <- data.frame(
Instance = c("Graph Coloring", "Clique Coloring (80,10,10)",
"Ramsey R(5,5,5)", "Scheduling (500 jobs)",
"Parity (CNFgen)", "Counting (CNFgen)",
"Matching (100 nodes)", "Ramsey (17,4)"),
Clauses = c(650000, 354890, 402752, 64570,
485200, 78760, 400, 4760),
Result = c("100%", "100% (5/5 seeds)", "100%", "100%",
"100%", "100%", "100%", "100%"),
Time = c("4.6s", "3.5s", "7.5s", "172ms",
"2.6s", "0.5s", "21ms", "<1ms"),
Hardware = c("Ryzen 5 5600H", "Laptop", "Laptop", "Laptop",
"Laptop", "Laptop", "Laptop", "Laptop")
)
kable(benchmarks,
caption = "NitroSAT performance on structured satisfiable instances. All results represent 100% clause satisfaction.",
align = c("l", "r", "c", "r", "l")) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
row_spec(1:3, bold = TRUE)| Instance | Clauses | Result | Time | Hardware |
|---|---|---|---|---|
| Graph Coloring | 650000 | 100% | 4.6s | Ryzen 5 5600H |
| Clique Coloring (80,10,10) | 354890 | 100% (5/5 seeds) | 3.5s | Laptop |
| Ramsey R(5,5,5) | 402752 | 100% | 7.5s | Laptop |
| Scheduling (500 jobs) | 64570 | 100% | 172ms | Laptop |
| Parity (CNFgen) | 485200 | 100% | 2.6s | Laptop |
| Counting (CNFgen) | 78760 | 100% | 0.5s | Laptop |
| Matching (100 nodes) | 400 | 100% | 21ms | Laptop |
| Ramsey (17,4) | 4760 | 100% | <1ms | Laptop |
scaling <- data.frame(
instance = c("Matching", "Ramsey 17,4", "Scheduling 500", "Counting",
"Clique Color", "Ramsey R(5,5,5)", "Parity", "Graph Color"),
clauses = c(400, 4760, 64570, 78760, 354890, 402752, 485200, 650000),
time_s = c(0.021, 0.001, 0.172, 0.5, 3.5, 7.5, 2.6, 4.6)
)
ggplot(scaling, aes(x = clauses, y = time_s)) +
geom_point(size = 3, color = "#C0392B") +
geom_text(aes(label = instance), hjust = -0.1, vjust = -0.5, size = 3.2) +
scale_x_log10(labels = comma, breaks = c(100, 1000, 10000, 100000, 1000000)) +
scale_y_log10(labels = function(x) ifelse(x < 1, paste0(x*1000, "ms"), paste0(x, "s"))) +
labs(
title = "NitroSAT Scaling: Clause Count vs. Solve Time",
subtitle = "All instances solved to 100% satisfaction (log-log scale)",
x = "Number of Clauses (log scale)",
y = "Solve Time (log scale)"
) +
theme_nitrosat() +
expand_limits(x = 2000000)
Clause count vs. solve time for perfectly solved instances. NitroSAT demonstrates near-linear scaling on structured satisfiable problems up to 650K clauses.
5.2 Scheduling Benchmarks
The scheduling instances encode real-world job scheduling with resource constraints: each job has a limited set of available time slots, and conflicting jobs cannot share a slot. These are formulated as list coloring problems in CNF.
sched <- data.frame(
Jobs = c(30, 50, 100, 200, 500, 1000),
Slots = c(5, 6, 8, 10, 10, 10),
Density = c(0.3, 0.2, 0.1, 0.05, 0.03, 0.02),
Clauses = c(1095, 2522, 7516, 21150, 64570, 154760),
Satisfaction = c(100, 100, 100, 100, 100, 99.99),
Time_ms = c(28, 22, 26, 85, 172, 225000)
)
p_sched <- ggplot(sched, aes(x = Jobs, y = Time_ms)) +
geom_col(fill = "#2E75B6", width = 0.6) +
geom_text(aes(label = paste0(Satisfaction, "%")), vjust = -0.5, fontface = "bold", size = 3.5) +
scale_y_log10(labels = function(x) {
ifelse(x < 1000, paste0(x, "ms"), paste0(x/1000, "s"))
}) +
labs(
title = "Scheduling (List Coloring) Benchmarks",
subtitle = "100% satisfaction through 500 jobs; 1000-job instance detected as UNSAT",
x = "Number of Jobs",
y = "Solve Time (log scale)"
) +
theme_nitrosat()
print(p_sched)
Scheduling benchmark results. NitroSAT achieves 100% satisfaction up to 500 jobs, with the 1000-job instance correctly identified as structurally unsatisfiable.
The 1000-job instance achieves 99.99% satisfaction and is correctly flagged as structurally unsatisfiable — NitroSAT’s thermodynamic phase transition detection identifies the impossibility rather than continuing to search indefinitely.
5.3 XOR-SAT Stress Test
XOR-SAT represents a critical stress test for continuous relaxation solvers. The flat energy landscape of XOR constraints produces zero gradient almost everywhere, defeating standard continuous methods. NitroSAT’s persistent homology module detects the cycle structure inherent in XOR constraints:
xor_results <- data.frame(
Instance = c("XOR (SAT, 200 clauses)", "XOR (planted, 8K clauses)", "XOR (UNSAT, 2K clauses)"),
Clauses = c(200, 8000, 2000),
Result = c("100%", "100%", "95.15%"),
Time = c("3.9ms", "9.5ms", "802ms"),
Beta1_Cycles = c("98 detected", "trivial", "2-26 cycles"),
Complexity = c("1.0 (maximum)", "low", "high")
)
kable(xor_results,
caption = "XOR-SAT stress test results. β₁ detects cycle structure; complexity score reaches maximum on the satisfiable XOR instance.",
align = c("l", "r", "c", "r", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Instance | Clauses | Result | Time | Beta1_Cycles | Complexity |
|---|---|---|---|---|---|
| XOR (SAT, 200 clauses) | 200 | 100% | 3.9ms | 98 detected | 1.0 (maximum) |
| XOR (planted, 8K clauses) | 8000 | 100% | 9.5ms | trivial | low |
| XOR (UNSAT, 2K clauses) | 2000 | 95.15% | 802ms | 2-26 cycles | high |
The detection of \(\beta_1 = 98\) independent cycles on the 200-clause XOR instance, combined with solution in 52 gradient steps, demonstrates that the persistent homology module provides genuine structural insight rather than serving as a passive diagnostic. The heat kernel propagates information along the parity chain structure identified by the topological analysis, enabling the continuous solver to navigate a landscape that would be impenetrable to gradient-only methods.
5.4 UNSAT Awareness
mirage <- data.frame(
Instance = c("Mirage 200", "Mirage 300"),
Peak_Satisfaction = c(85.2, 95.9),
Detection = c("Trap Detected", "Structural Awareness")
)
ggplot(mirage, aes(x = Instance, y = Peak_Satisfaction, fill = Detection)) +
geom_col(width = 0.5) +
geom_hline(yintercept = 100, linetype = "dashed", color = "red", linewidth = 0.5) +
annotate("text", x = 2.4, y = 100, label = "100% (impossible)",
color = "red", vjust = -0.5, size = 3.5) +
geom_text(aes(label = paste0(Peak_Satisfaction, "%")), vjust = -0.5, fontface = "bold") +
scale_fill_manual(values = c("#E67E22", "#C0392B")) +
coord_cartesian(ylim = c(0, 110)) +
labs(
title = "UNSAT Awareness: The Mirage Trap",
subtitle = "NitroSAT detects structural impossibility via thermodynamic phase transitions",
x = "",
y = "Peak Clause Satisfaction (%)"
) +
theme_nitrosat()
NitroSAT’s UNSAT awareness on ‘mirage’ instances designed to have high near-satisfaction but no complete solution. The solver detects structural impossibility and terminates early.
6 Discussion and Limitations
NitroSAT demonstrates that physics-informed mathematical design choices — number-theoretic clause weighting (Iyer 2026b), spectral pre-conditioning (Iyer 2026c, 2025b), topological cycle detection, and statistical mechanical phase transition awareness (Iyer 2025a) — can produce a MaxSAT solver that achieves state-of-the-art approximation quality on hard structured instances. The solver’s principal strengths are its performance on large structured satisfiable problems (graph coloring, Ramsey theory, scheduling) and its ability to detect unsatisfiable instances via thermodynamic signatures.
Several limitations warrant discussion. First, NitroSAT is not a complete solver: it cannot certify unsatisfiability, and its UNSAT detection is heuristic rather than provably sound. Second, NitroSAT does not consistently reach exact 100% on random 3-SAT near the phase transition [\(\alpha \approx 4.27\); Mitchell, Selman, and Levesque (1992); Achlioptas and Coja-Oghlan (2008); Mézard, Mora, and Zecchina (2005)], but post-release stress tests show stable high-quality approximation in this regime: 80 instances across \(n=300,500,1000\) all remained above 99.37% with scale-stable means (99.65%, 99.64%, 99.65%). Third, the current implementation is single-threaded LuaJIT; a C++ or Rust implementation with SIMD vectorisation and multi-threading would substantially improve wall-clock performance on industrial-scale instances.
Future directions include formal analysis of the partition function factorisation’s approximation quality, integration with clause-learning for hybrid complete/incomplete solving, and extension to weighted partial MaxSAT with explicit hard/soft clause separation.
7 Post-Release Empirical Validation
After the initial release, NitroSAT was subjected to extensive adversarial stress testing across multiple dimensions: scale, randomness, permutation invariance, and novel problem classes. All tests used identical solver parameters — no re-tuning was performed.
7.1 Full Benchmark Suite (358 Instances)
NitroSAT was evaluated across 358 benchmark instances spanning 19 problem types, including 5 novel problem types (N-Queens, Hamiltonian Cycle, Mutilated Chessboard, Planted 3-SAT, Exact Cover) that the solver had never been tuned for.
| Metric | Value |
|---|---|
| Instances tested | 358 |
| Solved at 100% | 115 (32.1%) |
| At 99%+ | 340 (95.0%) |
| At 96%+ | 353 (98.6%) |
| Average satisfaction | 99.58% |
| Largest perfect solve | 354,890 clauses (cliquecol) |
| Fastest 10K+ clause perfect solve | 22,521 clauses in 33ms |
Full results: full_benchmark_results.csv.
7.2 Phase-Transition Robustness (Random 3-SAT)
Random 3-SAT at the critical ratio \(\alpha \approx 4.26\) is the standard worst-case regime — no spectral coherence, no community structure, no symmetry. 80 instances were generated across three scales:
| Scale | Seeds | Avg Satisfaction | Min | Max | Std Dev |
|---|---|---|---|---|---|
| \(n=300\) (1,278 clauses) | 50 | 99.65% | 99.37% | 99.84% | 0.11% |
| \(n=500\) (2,130 clauses) | 20 | 99.64% | 99.44% | 99.81% | 0.10% |
| \(n=1000\) (4,260 clauses) | 10 | 99.65% | 99.53% | 99.72% | 0.06% |
80/80 instances above 99%. The mean does not degrade from \(n=300\) to \(n=1000\). Variance decreases with scale — this is scale-invariant approximate MaxSAT behavior on the hardest random regime.
7.3 Permutation Invariance
Variable renumbering and clause shuffling across 20 random permutations of the same instance:
| Metric | Value |
|---|---|
| Permutations tested | 20 |
| Perfect (100%) | 20/20 |
| Standard deviation | 0.0000% |
Zero variance. The solver operates on spectral properties of the clause hypergraph, not syntactic encoding order. This is a direct architectural consequence of diffusion- and topology-based dynamics.
7.4 Quasigroup / Latin Square Completion
A new algebraic CNF class — Latin-square completion (quasigroup) — was added with paired SAT/UNSAT instances at orders \(n=7,9,11,13\). This class features heavy algebraic symmetry and hard combinatorial structure fundamentally different from graph-coloring or Ramsey families.
Single-seed benchmark (8 instances):
| Metric | Value |
|---|---|
| Instances | 8 (4 SAT + 4 UNSAT) |
| Average satisfaction | 99.9763% |
| At 99.9%+ | 8/8 |
| SAT average | 99.9850% |
| UNSAT average | 99.9675% (4/4 plateau below 100%) |
5-seed distribution sweep (40 runs):
| Metric | Value |
|---|---|
| Average satisfaction | 99.9599% |
| Standard deviation | 0.0360% |
| At 99.9%+ | 38/40 |
| Min / Max | 99.8154% / 100.0000% |
| SAT runs at 99.9%+ | 20/20 |
| UNSAT runs below 100% | 20/20 |
Results: quasigroup_results.csv, quasigroup_seed_sweep.csv.
7.5 Zero-Tuning Generality on Novel Problem Types
Same code, same parameters, 100% on never-seen problem types:
| Problem | Clauses | Result | Time |
|---|---|---|---|
| N-Queens 25 | 24,825 | 100% | 3.7s |
| Exact Cover | 1,461 | 100% | 4ms |
| Planted 3-SAT (ratio 4.26) | 1,065 | 100% | 5ms |
| Hamiltonian Cycle | 13,240 | 99.99% | — |
| Mutilated Chessboard (UNSAT) | — | 99.42% (correct plateau) | — |
Planted List Coloring — 4/4 perfect across all scales:
| Vertices | Colors | Clauses | Result | Time |
|---|---|---|---|---|
| 50 | 5 | 699 | 100% | 105ms |
| 200 | 6 | 10,285 | 100% | 1.14s |
| 500 | 8 | 35,817 | 100% | 61ms |
| 1,000 | 10 | 122,090 | 100% | 119ms |
122,090 clauses, 10,000 variables, solved at step 66 in 119ms.
8 Reproducibility
NitroSAT is open-source under the Apache 2.0 license. The complete source code, benchmark scripts, and test instances are available at:
Repository: https://github.com/sethuiyer/NitroSAT