Header-only C++20 library for arithmetic in ℤ₂[[t]]: carry chains, formal derivatives, forward differences, Witt vectors, basis conversions, power series composition and reversion. Umbrella header dyadic.h includes all sub-headers from include/dyadic/.
#include<dyadic.h>
#include<cstdio>usingnamespacedyadic;intmain() {
Polynomial<4, uint64_t> p{{1, 2, 3, 4}};
// Formal derivative: D(P)(t) = dP/dtauto dp = formal_derivative(p); // {2, 6, 12, 0}for (auto c : dp) std::printf("%lu ", c);
// Forward difference: Δ(P)(t) = P(t+1) − P(t)auto fp = forward_difference(p); // {9, 18, 12, 0}// Basis conversions (exact roundtrip)auto ff = change_basis<FallingFactorialBasis>(p);
auto back = change_basis<MonomialBasis>(ff); // == p// Evaluate P(5) = 1 + 2·5 + 3·25 + 4·125 = 586auto y = p.eval(5);
// Witt vector ring operations
WittVector<3, uint64_t> a{{1, 2, 3}}, b{{4, 5, 6}};
auto sum = a + b; // Witt addition (via ghost map)auto prod = a * b; // Witt multiplication// Power series reversion (Lagrange inversion)
Polynomial<6, uint64_t> P{{0, 1, 1}}; // t + t²auto Q = reversion(P);
// Q = t − t² + 2t³ − 5t⁴ + 14t⁵ − 42t⁶ + …// (stored as unsigned ℤ₂ values: −1 ≡ 2⁶⁴−1, etc.)
}- C++20 compiler (GCC 12+, Clang 17+, MSVC 2022+)
- No external dependencies
detail::uint128_t(software 128-bit pair) providesuint64_tword support withoutunsigned __int128- Clang needs
-fconstexpr-steps=50000000for compile-time proofs; GCC needs-fconstexpr-ops-limit=200000000
| Area | What |
|---|---|
| 2-Adic Primitives | v2 (valuation), modinv_odd, div_2k_adic, artin_schreier (℘(x)=x²−x) |
| Polynomials | Polynomial<N,W,Basis> with eval, +, −, *, basis conversion (Monomial / FallingFactorial / Taylor), GCD, resultant, discriminant, pseudo-remainder, divmod, square-free check |
| Calculus | formal_derivative (D), forward_difference (Δ), taylor_shift, indefinite_sum (Σ = Δ⁻¹) |
| Witt Vectors | WittVector<N,W> with ghost map, Frobenius, Verschiebung, +, *, exp, log, inverse, adams_operation, teichmueller_lift |
| Carry Chain | Full-width carry propagation C = (I−N)⁻¹ — converges in one pass |
| Compose / Reversion | Power series composition P(Q(t)) and Lagrange inversion |
| Compile-Time Proofs | 24 named static_assert proofs (~55 total assertions) verifying ring axioms, basis roundtrips, D∘Δ=Δ∘D, ghost homomorphism, carry idempotence, Witt exp/log roundtrip (see dyadic/verify.h) |
| Combinatorial | binom, stirling_2, stirling_1, stirling_1_unsigned — all constexpr, cached at compile time |
WittVector<3, uint64_t> a{{1, 2, 3}}, b{{4, 5, 6}};
auto sum = a + b;
auto prod = a * b;
// Frobenius and Verschiebungauto fa = a.frobenius();
auto vb = b.verschiebung();
// F(V(a)) == V(F(a)) — verified at compile time// Adams operation ψⁿauto psi = adams_operation(a, 3);
// Teichmüller lift: τ(x) = (x, 0, 0, …) (N is required)auto tau = teichmueller_lift<4>(uint64_t{7});
// ghost_j(τ(ab)) = ghost_j(τ(a)) · ghost_j(τ(b)) — verified at compile time// NOTE: operator* uses carry-chain (2-adic ring).// poly_divmod_cw / poly_gcd_cw use coefficient-wise (standard ring).// Use poly_mul_cw() and verify_divmod_cw() for verification.
Polynomial<3, uint64_t> A{{1, 2, 1}}; // (x+1)²
Polynomial<2, uint64_t> B{{1, 1}}; // (x+1)auto [Q, R] = poly_divmod_cw(A, B); // Q=(x+1), R=0bool ok = verify_divmod_cw(A, Q, B, R); // true: A=Q·B+R ✓auto G = poly_gcd_cw(A, B); // G=(x+1) — divides bothPolynomial<6, uint64_t> P{{0, 1, 1}}; // P(t) = t + t²
Polynomial<6, uint64_t> Q{{0, 1, 2}}; // Q(t) = t + 2t²autoPQ = compose(P, Q); // P(Q(t)) — degree (N-1)*(M-1)+1 = 26auto R = reversion(P); // Lagrange inverse: P(R(t)) = t// R(t) = t − t² + 2t³ − 5t⁴ + 14t⁵ − 42t⁶ + …// (stored in ℤ₂: negative coefficients wrap modulo 2^W)#include<dyadic/verify.h>// triggers all static_asserts at compile time// If it compiles, all 22 named static_assert proofs passed.Day-to-day compiles include a sampled subset. For the full exhaustive suite (256² cases, 8K+ multiplication cases), define DYADIC_HEAVY_PROOFS.
Polynomial<6, uint8_t> p{{0, 0, 0, 0, 0, 255}};
if (!check_taylor_roundtrip_precision(p))
// T₅ = 5! · 255 = 30600 > 256 — high bits lost
WittVector<4, uint32_t> w{{1, 2, 3, 4}};
if (!check_witt_recovery_precision(w))
// Some rⱼ ≥ 2³²⁻ʲ — ghost recovery inexact| File | What |
|---|---|
docs/precision.md | Five precision windows unified under a single principle |
docs/theory.md | Four-domain unified theory: operator calculus ↔ Witt ghosts ↔ Mersenne ghosts ↔ thermodynamic classification |
As a header collection — copy dyadic.h and the include/dyadic/ directory into your project's include path, then #include <dyadic.h> (umbrella) or #include <dyadic/core.h> (individual header).
With CMake:
cmake -B build -DDYADIC_BUILD_TESTS=ON
cmake --build build
ctest --test-dir buildOptional: -DDYADIC_HEAVY_PROOFS=ON for exhaustive compile-time proofs.
Via CI script:
./ci_compile_check.sh| File | What |
|---|---|
test_core.cpp | Core axiom unit tests: six axioms, carry chain, arithmetic, evaluation |
test_verify.cpp | 24 named compile-time proofs (~55 assertions) + runtime verification across 9 (N,W) combos |
test_property.cpp | Randomized property-based tests: 18 invariants × 10 (N,W) combos |
test_full.cpp | 17 functional test groups covering the entire API surface |
test_negatives.cpp | Fork-based assertion verification (6 precondition checks) |
benchmark.cpp | Runtime benchmarks for key operations (build manually: g++ -O2 -std=c++20 -I. -Iinclude test/benchmark.cpp) |
All tests pass under GCC 14+ and Clang 17+ with ASan+UBSan. CI covers GCC (light + heavy proofs), Clang, and MSVC on Windows.
- Two ring semantics:
operator*uses carry-chainpoly_mul(correct ℤ₂ arithmetic with overflow propagation). Functions with_cwsuffix (poly_divmod_cw,poly_gcd_cw,poly_mul_cw) operate coefficient-wise (standard ring). These are different rings — usepoly_mul_cw()andverify_divmod_cw()to verify divmod/gcd results. quad_width<W>alias: Shorthand forwiden_t<dword_t<W>>(up to 4× word width), used by ghost-map accumulation and unsaturated polynomial products.assertpreconditions: All silent-failure points (poly_divmod_cwwith zero/even divisor,witt_expwith v₂(a₀) < 2,witt_log/witt_inversewith even a₀,reversionwith even P[1]) now useassert(). Invalid inputs triggerSIGABRTin debug builds. A few functions return 0 for invalid inputs by design rather than asserting:poly_discriminant_cwreturns 0 when degree < 1 or leading coefficient is even;det_laplacereturns 0 for dim > 6 (use Bareiss viadeterminant_cwinstead).- Compile-time cached Stirling/Pascal tables:
detail::STIRLING_CACHE<N,W>anddetail::PASCAL_CACHE<N,W>areinline constexprglobals — computed once per (N,W) at compile time, shared across all basis conversion and difference calls. - Auto-vectorization hints:
poly_mul_unsaturatedandpoly_muluseDYADIC_RESTRICT(__restrict__) and#pragma GCC ivdepon their inner multiply-accumulate loops, enabling the compiler to generate SIMD code for runtime invocations without breakingconstexpr.
See docs/complexity.md for a complete function-by-function breakdown.
| Operation | Time | Space | Notes |
|---|---|---|---|
v2, modinv_odd, exact_divide, div_2k_adic | O(1) | O(1) | Single-instruction or fixed-iteration 2-adic primitives |
eval (Monomial) | O(N) | O(1) | Horner's method |
eval (FF / Taylor) | O(N) | O(1) | Incremental falling product / binomial coefficient |
formal_derivative (Monomial) | O(N) | O(N) | Direct coefficient scaling |
formal_derivative (FF / Taylor) | O(N²) | O(N) | Convert → monomial D → convert back |
forward_difference (FF / Taylor) | O(N) | O(N) | FF diagonalizes Δ: scale+shift; Taylor: pure shift |
forward_difference (Monomial) | O(N²) | O(N) | Pascal-cache binomial transform |
taylor_shift | O(N²) | O(N) | Pascal-cache / falling-factorial binomial transform |
indefinite_sum (FF) | O(N) | O(N) | Σ(FFₙ₋₁) = FFₙ / n |
indefinite_sum (Monomial / Taylor) | O(N²) | O(N) | Converts to FF, sums, converts back |
poly_exp / poly_log | O(N²) | O(N) | Double-loop recurrence with dword intermediates |
poly_mul / operator* | O(N·M) | O(1) chunked | Naive convolution + carry chain; chunked 256-bit buffer; unsigned __int128 for uint64_t (~2.7× speedup deg 63) |
mul_unsigned | O(N·M) | O(1) chunked | Big-integer unsigned multiply (N+M limbs); same chunked unsigned __int128 pattern as poly_mul |
poly_mul_cw | O(N·M) | O(N+M) | Standard coefficient-wise convolution |
compose | O(N²·M²) | O(N·M) | Naive power-series accumulation; result ≤ 4095 coeffs |
reversion | O(N³) | O(N²) | Incremental Lagrange inversion (not Newton) |
change_basis | O(N²) | O(N) | Stirling-cache triangular matrix multiply; 1–2 conversions per call |
witt_add / witt_mul | O(N²) | O(N²) | Ghost-map + Newton recovery; quad_width accumulators |
witt_log / witt_exp | O(N² + N·T) | O(N²) | T ≈ 2·bit_width (ghost p-adic series terms) |
witt_inverse | O(N²) | O(N²) | Ghost inverse via modinv_odd per component |
adams_operation | O(N² + N·n) | O(N²) | Ghost ^ n powering |
poly_divmod_cw | O(N·(N+M)) | O(max(N,M)) | Long division (requires odd lc) |
pseudo_remainder_cw | O(N·(N+M)) | O(N) | Subresultant PRS (no lc requirement) |
poly_gcd_cw | O(K³) worst | O(K) | Euclidean PRS; K = max(N,M) |
polynomial_resultant_cw | O(1) bounded | O(1) | Sylvester matrix (dim ≤ 6) + Laplace det (≤ 6! = 720) |
poly_discriminant_cw | O(N) + O(1) | O(N) | Derivative + resultant |
poly_is_square_free_cw | O(N³) | O(N) | gcd(P, P′) |
Matrix::operator* (M×N × N×P) | O(M·N·P) | O(M·P) | i-k-j loop with zero-skip |
Matrix::determinant | O(M³) | O(M·N) | Bareiss fraction-free elimination |
Matrix::inverse / solve / rref | O(M³) | O(M²) | Gauss–Jordan with odd-pivot modinv_odd |
pade_approximant [M/N] | O(N³ + M·N) | O(N² + M+N) | Gaussian elimination on N×N Toeplitz system |
cf_convergent | O(n²) | O(n) | Classical recurrence; degree grows linearly |
div_unsigned | O(NL·NR) | O(NL+NR) | Knuth Algorithm D (NR ≥ 2) or divmod_single (NR=1) |
reciprocal_newton | O(NR²·log NR) | O(NR) | Newton iteration; ceil(log₂((NR+1)·W)) iterations |
div_newton | O(NR²·log NR + NL·NR) | O(NL+NR) | Newton reciprocal division with correction step |
binom / stirling_2 / stirling_1 | O(k·log n) / O(n²) | O(1) / O(n) | GCD-multiplicative / DP recurrence |
StirlingCache / PascalCache ctor | O(N²) compile-time | O(N²) | inline constexpr cached once per (N,W) |
All operations are constexpr; runtime performance matches compile-time complexity bounds. The carry-chain poly_mul and mul_unsigned are auto-vectorized with SIMD hints (#pragma GCC ivdep, DYADIC_RESTRICT) and use hardware unsigned __int128 accumulation for uint64_t where available.
Polynomial::degree()renamed tomax_degree(): The old name was misleading (it returnedN−1, the maximum possible degree, not the actual degree). Addedactual_degree()member function.- Taylor basis roundtrip:
T_k = k! · FF_kwraps whenFF_k ≥ 2^W / k!. Use small coefficients for exact roundtrips. FallingFactorial basis has no such limitation. - Witt precision window: Recovery
r_j = (G_j − S_j) / 2^jrequiresr_j < 2^{W−j}. - Witt exp/log term truncation: Uses valuation-aware dynamic term counting with 2× bit-width budget. Requires
v₂(x) ≥ 2for exp convergence (mathematical limit —v₂(x) = 1stalls at≤ log₂(n)+1). Log converges forv₂(y) ≥ 1(~135 terms at 128-bit). detail::uint128_tis a software 128-bit pair — nounsigned __int128required.__int128is used as a hot-path optimization inbinom()(dyadic/combinatorial.h),poly_mul()andmul_unsigned()(dyadic/core.h/arith.h, ~2.7× speedup at deg=63 for uint64_t),div_unsigned_knuth()(dyadic/arith.h, hardware trial-division / multiply-subtract for ~3-4× speedup), and the p-adic series inwitt_log/witt_exp(dyadic/witt.h, hardware multiply-accumulate). All paths guarded by__SIZEOF_INT128__.
This project is made available under the terms of the MIT License.