The backtracking search in the paper computes const(p) exactly up to p = 14. Replacing the enumeration by a dynamic programme over subsets, carried out modulo primes, reached p = 18: const(18) has 463 digits and was computed on a single node of the Hábrók cluster. The walkthrough below builds up both steps for small p, and the code follows it.
The step-by-step walkthrough requires JavaScript. In short: the backtracking search calls itself once for every partial permutation, while the dynamic programme computes one value for each set of values still to be placed. At p = 5 that is 156,297 calls against 282 states. Each state is stored as a bitmask, each value modulo primes, and const(p) is rebuilt by the Chinese remainder theorem.
The walkthrough is limited to p ≤ 5. The solvers in the repository below handle the larger cases.
The solvers are available at
kiancshah/wronskian-constant
under the MIT licence. Two independently written solvers agree residue by residue.
verify.py compares them with each other, with reference.py (a direct
implementation of the definition) and with the stored values. identity_check.py
expands the operator identity symbolically and recovers const(1) = 1 and const(2) = 2 without
using the recursion, which checks the recursion against the original definition.
wronskian-constant/ ├── constp.cpp flat array over every subset of steps ├── constp2.cpp the same recursion, one layer at a time ├── driver.py picks primes, runs a solver, rebuilds by CRT ├── reference.py the definition, transcribed directly ├── verify.py solvers against each other, reference and data ├── identity_check.py symbolic expansion of the identity ├── data/const.txt computed values ├── docs/ algorithm.md, equivalence.md └── cluster/ Slurm scripts, one prime per task
The algorithm of the paper, shown step by step on the research page, reaches p = 14. For larger p the enumeration of permutations is replaced by a dynamic programme over subsets. It stores only the states that are needed, reduces each transition to a table lookup and an addition, and uses a symmetry of the recursion (reversal combined with negation) to stop halfway through the computation. About half of all subset states are reached for every p tested, so the speed-up comes from making each state cheap rather than from sparsity.
The computation is carried out modulo primes below 229, and the integer is reconstructed by the Chinese remainder theorem. An improved proved bound on the size of const(p) reduces the number of primes needed at p = 18 from 86 to 63. Per prime, the solver is 15 to 20 times faster than the previous one and uses about a fifth of the memory; on two cores, const(15) takes 37 seconds instead of 12 minutes. const(16) and const(17) were computed in minutes on two cores, and const(18), with 463 digits, on a single node of the Hábrók cluster.