Divstep: Division-Free Algorithm for Reversible Modular Inversion

• Divsteps is an efficient algorithm providing modular inversion via parity operations, additions, subtractions, and value exchanges without traditional division, making it suitable for elliptic-curve cryptography, maintaining determinant invariant `fg_v - gf_v = p` where `f_v` and `g_v` track contributions from modulo operations.
• Implements arithmetic-based branch history recovery for undoing decisions of modular inversion execution, history stack manages ambiguous decisions, with parity counters ensuring accurate measurement of branch diverging.
• Architecture allows constructing a modular inversion circuit using 4.0x10^6 Toffoli gates and 935 qubits, maintaining local compressed execution history via reversible computation logic, increasing parallel performance and minimizing computational footprint.

Divsteps is a division-free variant of the extended Euclidean algorithm that computes modular inverses through parity tests, additions or subtractions, division by two, and occasional exchanges of remainder-like values. In reversible computation, divsteps provides the basis for a prime-field modular-inversion architecture discovered through evaluator-guided large-language-model search. The architecture uses arithmetic-state decoding, bounded branch history, measurement-based uncomputation, packed registers, freeze-and-reverse termination handling, and reverse-walk modular halving. For the 256-bit secp256k1 prime, it yields a complete reversible elliptic-curve point adder using $4.0\times10^6$ Toffolis and 935 qubits, while accepting a bounded failure probability (Gu et al., 15 Oct 2026).

1. Classical algorithm and modular inversion

Given an odd modulus $p$ and a nonzero integer $w<p$, divsteps computes

$w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$

The algorithm is a division-free form of the extended Euclidean algorithm. The ordinary extended Euclidean algorithm maintains integers and coefficients that eventually satisfy

$\alpha p+\beta w=1.$

Reduction modulo $p$ then gives $\beta w\equiv1\pmod p$, so $\beta$ is a modular inverse of $w$.

Divsteps replaces multi-bit quotient calculations and repeated magnitude comparisons with elementary operations. An odd remainder-like value $f$ is retained while another value $p$0 is reduced. If $p$1 is even, it is halved. If $p$2 is odd, $p$3 is added to or subtracted from $p$4 before halving, making the result even. A signed counter $p$5, also called the half-delta, determines when the roles of $p$6 and $p$7 are exchanged without explicitly comparing their magnitudes.

The classical divstep state is

$p$8

initialized as

$p$9

The auxiliary values $w<p$0 and $w<p$1 track the contribution of $w<p$2 in a scale-dependent form. After $w<p$3 active steps,

$w<p$4

The three update cases are

$w<p$5

These cases can be represented by two decision bits: $w<p$6, indicating whether the step swaps the roles of $w<p$7 and $w<p$8, and $w<p$9, indicating whether the incoming $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$0 is odd. Their values are respectively $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$1, $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$2, and $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$3 for the three cases.

The walk terminates when $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$4. Assuming $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$5, the resulting $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$6 is either $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$7 or $w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$8. Divsteps preserves the determinant invariant

$w^{-1}\pmod p, \qquad ww^{-1}\equiv 1\pmod p.$9

At termination this gives $\alpha p+\beta w=1.$0. Together with the scaling invariant, it follows that

$\alpha p+\beta w=1.$1

Thus divsteps computes the inverse up to a power-of-two scale. For the benchmark modulus, safeGCD analysis guarantees termination within 590 active divsteps for every $\alpha p+\beta w=1.$2.

2. Elliptic-curve application and reversible-computation obstacles

The benchmark uses the 256-bit secp256k1 prime

$\alpha p+\beta w=1.$3

The modular inversion occurs in affine elliptic-curve point addition. For a point $\alpha p+\beta w=1.$4 and a classically supplied point $\alpha p+\beta w=1.$5, the slope is

$\alpha p+\beta w=1.$6

followed by

$\alpha p+\beta w=1.$7

The denominator inversion dominates the point-adder cost.

A classical divstep overwrites its state, whereas a reversible circuit must preserve enough information to reconstruct every overwritten value. A direct reversible implementation would need to retain branch decisions, overwritten arithmetic values, termination information, variable-width and sign information, temporary carries, and comparison results. Four independent 256-bit registers for $\alpha p+\beta w=1.$8, $\alpha p+\beta w=1.$9, $p$0, and $p$1 would already require approximately 1,024 qubits before counters, ancillas, and history.

Reversible implementation introduces several further complications. Branch predicates must be computed, used, and erased while all three divstep cases remain represented coherently. Swaps cannot be treated as ordinary classical variable exchanges whose identities are forgotten. Signed values require two’s-complement handling, sign extension, arithmetic right shifts, and width management. Input-dependent stopping times require either a worst-case schedule or a mechanism that freezes terminated walks. Finally, the scale factor $p$2 depends on the input-dependent number of active steps.

The architecture addresses these issues without changing the divstep recurrence itself. Its modifications concern representation, control synthesis, history management, uncomputation, stopping-time handling, bounded-resource scheduling, and integration with point-addition arithmetic.

3. State-based branch recovery and excursion decoding

The central reversible principle is that the current arithmetic state contains information about many earlier branch decisions. Throughout the recurrence,

$p$3

After a swap,

$p$4

so $p$5 has exactly one trailing zero. After a nonswap,

$p$6

and because $p$7 is already even, $p$8 is divisible by four. Consequently,

$p$9

Every swap therefore becomes recoverable from the updated coefficient register, with no stored swap history.

Parity decisions can also be recovered locally in many cases. For a nonswap with $\beta w\equiv1\pmod p$0, the incoming $\beta w\equiv1\pmod p$1 must have been even, because an odd $\beta w\equiv1\pmod p$2 would have triggered a swap. If a nonswap first reaches $\beta w\equiv1\pmod p$3, the incoming $\beta w\equiv1\pmod p$4 must have been odd, since halving a nonzero even integer cannot produce zero. These tests recover roughly two thirds of the step decisions. The remaining ambiguous decisions are nonterminating nonswaps entered with $\beta w\equiv1\pmod p$5.

An ambiguous excursion begins with a swap followed by a sequence of ambiguous nonswaps. If

$\beta w\equiv1\pmod p$6

then the associated coefficient update has the form

$\beta w\equiv1\pmod p$7

Multiplication by $\beta w\equiv1\pmod p$8 modulo $\beta w\equiv1\pmod p$9 gives

$\beta$0

The relevant invariant is

$\beta$1

Therefore the top bit of the residue identifies the current branch bit:

$\beta$2

Using the determinant invariant, an equivalent decoder is

$\beta$3

The decoder depends only on the updated state and low-order bits of $\beta$4 and $\beta$5. The underlying window invariant states that if $\beta$6, then $\beta$7 divides $\beta$8; if $\beta$9 and $w$0 with $w$1 odd, then there exists $w$2 satisfying

$w$3

This prevents earlier decisions from contaminating the bit used to decode the current one. Immediately after the initiating swap, the window contains a leading 1, while subsequent ambiguous decisions contribute higher binary digits.

For small depths, the modular decoder is compiled into Boolean truth tables. The numbers of product terms at depths 2, 3, 4, and 5 are respectively 3, 3, 11, and 26; at depths 6, 7, and 8 they grow to 79, 268, and 1,092. This growth motivates a tunable boundary between decoding and explicit history storage.

4. Architecture, history, and uncomputation

The discovered architecture combines state-based branch recovery, a bounded history stack, measurement-based uncomputation, step-dependent register packing, freeze-and-reverse stopping-time handling, reverse-walk scale removal, and prime-specific modular arithmetic.

Let $w$4 be the storage threshold. Decisions with $w$5 are decoded from the arithmetic state; decisions with $w$6 are stored as history bits. The principal configuration uses

$w$7

In sampled walks, 50.1% of ambiguous decisions occur at depth 2, 25.0% at depth 3, 12.5% at depth 4, and 12.4% at depth 5 or greater. Consequently, only 12.4% of ambiguous decisions are stored in the main configuration. The average required history is approximately 21 bits, although the unrestricted worst case can require 263 stored decisions for the 533-step schedule. A stack of 60 decision bits is allocated, with overflow treated as a rare failure event having probability $w$8 per inversion.

During forward execution, the parity bit $w$9 controls the arithmetic update. Stored decisions are pushed onto the history stack. Decoded decisions are subjected to measurement-based uncomputation in the $f$0 basis. The measurement removes the qubit but leaves a phase dependent on $f$1 and the measurement result. During reverse execution, the decoder reconstructs $f$2, applies the corresponding phase correction, and uses the restored $f$3 parity to clear the residual branch information. The reverse step either pops a stored decision or reconstructs one from the updated arithmetic state before applying the inverse divstep transformation.

The arithmetic transformations, determinant identity, and branch-decoding rules are exact. Failures arise from deliberately bounded resource limits rather than approximation of the divstep algebra.

The architecture also packs data into two 256-qubit registers rather than maintaining four permanent 256-bit registers. One register shares storage between $f$4 and $f$5; the other shares storage between $f$6 and $f$7. Remainder-like values shrink while coefficient-like values grow, allowing growing coefficients to occupy positions vacated by shrinking remainders. The widths and locations are fixed when the circuit is generated; inputs exceeding those assigned widths count as failures.

5. Termination, scale removal, and prime-specific arithmetic

The circuit schedules 533 divstep positions, although the mathematical worst-case bound is 590 active steps. It freezes a walk when $f$8 using a frozen flag and a step-index register. Because a frozen state is identical at every later scheduled step, the state alone does not reveal when freezing began. The expensive zero test is enabled only after step 498, and the circuit accepts failures for walks that terminate before the zero-test window, fail to terminate by the end of the schedule, or overflow counters and integer widths.

Sampled stopping times are concentrated roughly between steps 505 and 530. The reported early-stopping probability is approximately $f$9 per inversion, while the $p$00-register overflow probability is approximately $p$01 per inversion.

The input-dependent factor $p$02 is removed during reverse execution rather than by retaining the active-step count. If the desired result is $p$03, termination produces

$p$04

Each active reverse step performs modular halving:

$p$05

The frozen flag identifies active reverse steps. After exactly $p$06 halvings, the result is $p$07, while the walk registers are restored and no explicit stopping-time register remains.

The sparse form of the secp256k1 prime enables specialized reduction:

$p$08

High product bits can therefore be folded into low bits using the small constant $p$09. The point adder fuses coordinate arithmetic into destination registers, avoids materializing a complete $p$10 register, and splits

$p$11

so that

$p$12

The three products are accumulated and cleared one at a time.

Carry propagation is truncated within the permitted failure budget. Modular-reduction carries propagate through only the low 58 bits, or 60 bits in the fused square, and comparison tests inspect only the top 27 bits. These truncations are resource-limited approximations in the complete circuit and are distinct from the exact divstep identities.

6. Resource frontier, comparisons, and limitations

The main architecture is embedded in a complete reversible affine point adder implementing

$p$13

with a classically supplied addend $p$14, cleaned ancillas, and the correct phase except on a specified bounded fraction of inputs. Its resource counts are

$p$15

with

$p$16

The four divstep walks contribute approximately $p$17 Toffolis, or 63% of the total. Surrounding field multiplication contributes $p$18 Toffolis, or 35%, while fused square and coordinate additions contribute approximately $p$19.

The history-decoding threshold produces an explicit Toffoli–qubit trade-off.

Configuration	Toffolis	Qubits
Larger-history / lower-decoding-cost	$p$20	1,031
Main $p$21 configuration	$p$22	935
Low-qubit $p$23 configuration	$p$24	893
Higher-qubit configuration	$p$25	1,191

At $p$26, every ambiguous decision is stored, requiring approximately 198 history qubits and producing a total of 1,065 qubits. At $p$27, decoders handle depths through 8, leaving 24 history qubits and producing the lowest-qubit point reported.

The divstep construction is not the lowest-Toffoli point-adder construction in the reported comparison. A binary-gcd construction with a dialog/history trick uses approximately $p$28 Toffolis and 1,446 qubits. The discovered divstep circuit uses more Toffolis but fewer qubits. Packed Euclidean approaches use 838 to 843 qubits but require approximately $p$29 to $p$30 Toffolis.

A second discovered construction specializes a packed Euclidean walk to the secp256k1 modulus. It uses $p$31 Toffolis and 843 qubits, compared with a reconstructed baseline using $p$32 Toffolis and 838 qubits. The specialization reduces Toffoli count by a factor of approximately 3.4 at the cost of five qubits. The divstep architecture is therefore the lower-Toffoli alternative among the two discovered constructions, whereas the specialized Euclidean architecture occupies the extreme low-qubit end.

The permitted failure probability is approximately $p$33 per point addition. The reported union-bound contributions are:

Failure source	Probability per point addition
Multiplier carry and comparison truncation	$p$34
Stopping-time and integer-size limits	$p$35
Early stopping and $p$36-register limits	$p$37
Branch-history overflow	$p$38
Fused-square carry limit	$p$39
Union bound	$p$40

In 100,000 random tests of the full point-addition circuit, four failures were observed. The point-addition input distribution samples $p$41 and $p$42 independently and nearly uniformly from curve points, excluding affine-formula exceptional cases. For a 16-bit-window algorithm invoking the point adder 28 times, a union bound gives less than 0.6% retry overhead.

The primary score is $p$43, with Toffoli and CCZ gates counted equally. This metric is a convenient uniform objective rather than a universal hardware criterion. Toffoli depth, active volume, magic-state throughput, or a more heavily qubit-weighted objective could be more appropriate under other fault-tolerant architectures. No complete modular-inversion circuit-depth or Toffoli-depth estimate is reported.

7. LLM-guided discovery and broader significance

The architecture was discovered through evolutionary search over Python programs that construct circuits rather than through ordinary gate-level peephole optimization. Candidate programs contain named quantum registers, reusable arithmetic functions, loops, control logic, ancilla allocation and freeing, reversible gates, measurement-based uncomputation, and the overall algorithmic organization.

The search procedure provided a hand-written task specification, requested several algorithm families from a LLM, evolved each family independently, periodically shared strong candidates across families, tested candidates with an independent evaluator, and retained valid circuits according to their resource scores. The search used OpenEvolve with Claude Mythos 5 and MAP-Elites. MAP-Elites preserved different points in the $p$44 resource plane rather than only the single minimum-$p$45 candidate.

The evaluator checked basis-state input/output relations, preservation of input registers, phase correctness, cleanup of all ancillas, measurement-record independence of the logical permutation, Toffoli count, and peak qubit count. The circuit model permits reversible classical gates, diagonal sign gates, and destructive $p$46-basis measurements with phase correction. A candidate is valid only if every ancilla returns to $p$47 and the measurement record factors from the logical output. The simulator follows computational-basis trajectories rather than full amplitudes because the allowed gates preserve computational-basis states up to phase.

The significance of the divstep result lies in its reversible architecture and representation rather than in a local gate rewrite. The current arithmetic state acts as a compressed, partially decodable record of the execution history. The determinant and scaling invariants provide global constraints, while the $p$48-adic structure of the coefficient registers supplies local branch information. Explicit history is retained only where decoding becomes more expensive than storage.

The construction demonstrates that divsteps, previously impractical for reversible integer arithmetic under a straightforward implementation, can be incorporated into a reversible prime-field inversion circuit for the secp256k1 benchmark. Its scope is nevertheless constrained by the prime-specific sparse reduction, the assumed input distribution, bounded failure probability, fitted stopping-time and carry limits, selected resource objective, and extensive state-decoding logic. It establishes a viable benchmark architecture rather than a universal result that divsteps dominates binary-gcd or Euclidean inversion for all moduli, operand sizes, or fault-tolerant hardware models.