NIKX-Tech · erfanjazebnikoo · Jul 12, 2026
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
@@ -48,6 +48,9 @@ and in `tests/test_runner.c` and `benchmarks/bench_runner.c`.
  * Mathematical definition:
  *   result = integral from a to b of f(x) dx
  *
+ * Reference: Author, Title, Venue, Year. Full citation:
+ * docs/algorithms/[module].md#references.
+ *
  * @param[in]  param1  Description, valid range.
  * @param[out] result  Pointer to store result. Must not be NULL.
  *
@@ -68,11 +71,39 @@ numx_status_t numx_function(numx_real_t param1, numx_real_t *result);
 
 ---
 
+## Citing sources
+
+Every function implementing a published algorithm cites that source inline in its
+docstring (short form, per the template above); the full bibliographic entry lives in
+the module's `docs/algorithms/[name].md#references`, not duplicated in the code. This
+applies equally to a function you wrote from a textbook and to a fix contributed by
+someone outside the core team (e.g. a community-submitted correctness or performance
+fix): the docstring should say what it's based on either way. `docs/algorithms/*.md`
+formula sections carry the same `[n]` marker inline, next to the formula, not only in
+a bottom References list, so a reader doesn't have to hunt for which citation backs
+which equation.
+
+Modules with an available external reference implementation or published dataset
+should include a `validation/reference/<module>/` directory: a small generator script
+(any language/tool) that produces a frozen test-vector fixture from that external
+source, checked into `tests/vectors/`, with a README documenting provenance (source,
+pinned commit/version, regeneration steps). This is validation tooling only; it is
+never a build or runtime dependency of the library itself. See
+`validation/reference/ntt/` for the reference instance of this pattern.
+
+---
+
 ## Test file structure
 
 Tests are organized into four levels:
 
-- **L1** — Known-answer tests: verified against analytical truth or a reference implementation.
+- **L1** — Known-answer tests: verified against analytical truth, or a reference
+  implementation. For crypto/security-critical modules, "reference implementation"
+  means an external, independently-citable one (published test vectors, or an
+  established third-party implementation), not code written for this project. An
+  internally-written reference (e.g. a naive O(n^2) implementation checked against an
+  optimized one) is still valid for general numerical modules, but is self-consistency,
+  not independent proof, and is not sufficient on its own for crypto modules.
 - **L2** — Property tests: mathematical invariants that must always hold.
 - **L3** — Edge cases: zero, negative, extremes, near-singular, empty input.
 - **L4** — Error handling: null pointers and invalid args return the correct error code.

diff --git a/ROADMAP.md b/ROADMAP.md
@@ -34,21 +34,36 @@ fixes: Cholesky decomposition (contributed by Erfan Esmaeili), a branchless
 Barrett reduction closing the NTT module's one known constant-time gap
 (contributed by u/robchroma on r/C_Programming, verified exhaustively
 against 22,164,483 inputs before merging), function/data-section garbage
-collection, and Kahan compensated summation across 8 hot loops. Remaining:
+collection, Kahan compensated summation across 8 hot loops, and cross-
+validating the `ntt` module bit-for-bit against PQClean's official reference
+implementation (closing a self-consistency-only test gap). Remaining:
 multi-platform Cholesky validation, a post-Kahan benchmark refresh, and the
 release itself.
 
 ---
 
 ## Next
 
+### [Academic paper](https://github.com/NIKX-Tech/numx/milestone/9)
+
+Promoted from backlog: an outside academic reviewer confirmed the algorithms
+and cross-platform validation methodology are publication-grade, provided
+every implementation is traceable to a citable scientific source (in the
+docstring, not just the docs) and independently validated, not just
+self-consistency-tested, the way the NTT/PQClean cross-validation above now
+does. Covers outline, target venue, a citation audit of the crypto modules,
+and an AI-assisted-development disclosure section.
+
 ### [v1.2.0 - Post-quantum crypto and rng](https://github.com/NIKX-Tech/numx/milestone/4)
 
 CRYSTALS-Kyber KEM and CRYSTALS-Dilithium signatures on top of `numx_ntt_*`,
 making numx a verified, allocation-free, pure C99 post-quantum stack for
 embedded systems. Includes a new `rng` module (portable PRNG plus a
 caller-supplied entropy hook) that also closes a known ESP32 test
-portability issue. Gated on an [open invitation for independent security
+portability issue. Architecture design and `rng` can proceed now; the Kyber
+and Dilithium implementation issues themselves are gated on the citation
+audit and reference-validation groundwork above landing first. Also gated on
+an [open invitation for independent security
 review](https://github.com/NIKX-Tech/numx/issues/68) of the `ntt` module.
 
 ### [v1.3.0 - Estimation and control](https://github.com/NIKX-Tech/numx/milestone/5)
@@ -69,7 +84,6 @@ require hand-derived Jacobians.
 
 - [Fixed-point DSP (q15/q31)](https://github.com/NIKX-Tech/numx/milestone/7) - vector ops for FPU-less targets (Cortex-M0, AVR)
 - [RISC-V hardware validation](https://github.com/NIKX-Tech/numx/milestone/8) - closing the gap between the platform list and what's actually been run on real silicon
-- [Academic paper](https://github.com/NIKX-Tech/numx/milestone/9) - documenting the algorithms and cross-platform validation methodology
 - [v2.0.0 - API evolution](https://github.com/NIKX-Tech/numx/milestone/10) - breaking changes, applied together with a migration guide (first candidate: a context parameter on `numx_func1d_t` callbacks)
 
 ---

diff --git a/docs/algorithms/ntt.md b/docs/algorithms/ntt.md
@@ -74,28 +74,30 @@ $$x^{256} + 1 = \prod_{k=0}^{127} (x^2 - \omega^{2k+1})$$
 
 Each factor $\mathbb{Z}_q[x]/(x^2 - c_k)$ is a degree-2 ring with modulus $c_k = \omega^{2k+1} \bmod q$.
 
-### Forward NTT (Cooley-Tukey)
+### Forward NTT (Cooley-Tukey) [2, 5]
 
 The forward NTT maps $f \in R_q$ to its 128 CRT residues in seven Cooley-Tukey
-butterfly stages. Each butterfly with twiddle factor $\zeta$:
+butterfly stages [5], specialized for negacyclic lattice-based cryptography
+following Longa & Naehrig [2]. Each butterfly with twiddle factor $\zeta$:
 
 $$f[j]       \leftarrow f[j] + \zeta \cdot f[j + \ell] \pmod{q}$$
 $$f[j + \ell] \leftarrow f[j] - \zeta \cdot f[j + \ell] \pmod{q}$$
 
 where $\ell$ is the current half-group length (128, 64, 32, 16, 8, 4, 2 across stages).
 The output pairs $(f[2i], f[2i+1])$ represent $f \bmod (x^2 - c_i)$ for $i = 0, \ldots, 127$.
 
-### Inverse NTT (Gentleman-Sande)
+### Inverse NTT (Gentleman-Sande) [1, 3]
 
-The inverse NTT uses Gentleman-Sande (decimation-in-frequency) butterflies, processing
+The inverse NTT uses Gentleman-Sande (decimation-in-frequency) butterflies [3], processing
 stages in reverse order ($\ell = 2, 4, \ldots, 128$). Each inverse butterfly:
 
 $$t           \leftarrow f[j]$$
 $$f[j]        \leftarrow t + f[j + \ell] \pmod{q}$$
 $$f[j + \ell] \leftarrow \zeta^{-1} \cdot (t - f[j + \ell]) \pmod{q}$$
 
 After all stages, every coefficient is multiplied by the normalization constant
-$n^{-1} = 128^{-1} \equiv 3303 \pmod{3329}$, recovering the original polynomial.
+$n^{-1} = 128^{-1} \equiv 3303 \pmod{3329}$, recovering the original polynomial,
+following the Kyber specification's inverse-NTT algorithm [1].
 
 **Key invariant:** within each inverse stage, butterfly groups must use the same twiddle
 indices as the corresponding forward stage (FORWARD order within each stage, not
@@ -114,16 +116,17 @@ Because the NTT is a linear map over $\mathbb{Z}_q$, these operations are valid
 
 Both functions apply Barrett reduction to keep outputs in $[0, q-1]$.
 
-### Pointwise multiplication (basemul)
+### Pointwise multiplication (basemul) [1]
 
 Once both polynomials are in NTT domain, their product is computed in O(n) time via
-128 independent degree-2 ring multiplications. For each pair $i$ with twiddle $c_i$:
+128 independent degree-2 ring multiplications, following the Kyber specification's
+basemul algorithm [1]. For each pair $i$ with twiddle $c_i$:
 
 $$(a_0 + a_1 x)(b_0 + b_1 x) \bmod (x^2 - c_i) = (a_0 b_0 + c_i a_1 b_1) + (a_0 b_1 + a_1 b_0)\, x$$
 
-### Barrett reduction
+### Barrett reduction [4]
 
-All modular multiplications use Barrett reduction to avoid the hardware division
+All modular multiplications use Barrett reduction [4] to avoid the hardware division
 instruction, which is slow or absent on embedded processors:
 
 $$\text{barrett}(a) = a - \left\lfloor \frac{(a \ggg 10) \cdot v}{2^{15}} \right\rfloor \cdot q, \qquad v = 10079$$
@@ -149,6 +152,18 @@ known data-dependent branches. Contributed by u/robchroma (r/C_Programming, 2026
 and exhaustively verified against `a % q` for all 22,164,483 valid inputs before
 adoption.
 
+### Validation provenance
+
+The known-answer tests in `tests/test_ntt.c` originally checked the fast NTT
+only against a naive O(n^2) reference multiplication written for this
+project: self-consistency, not proof against an independent ground truth.
+`numx_ntt_forward` and `numx_ntt_polymul` are now additionally cross-validated
+bit-for-bit against PQClean's `ml-kem-512` reference implementation (the
+FIPS 203-standardized successor to CRYSTALS-Kyber-512), an external,
+independently-citable implementation not written by this project. See
+`validation/reference/ntt/README.md` for the pinned commit, the Montgomery-
+domain reasoning, and regeneration steps.
+
 ---
 
 ## Complexity
@@ -214,6 +229,8 @@ $f \bmod (x^2 - c_i)$ for $i = 0, \ldots, 127$. Coefficients in $[0, q-1]$.
    *AFIPS Fall Joint Computing Conf.*, Vol. 29, pp. 563-578, 1966.
 4. Barrett, P. — "Implementing the Rivest Shamir and Adleman Public Key Encryption
    Algorithm on a Standard Digital Signal Processor", *CRYPTO 1986*, LNCS 263, 1987.
+5. Cooley, J. W. & Tukey, J. W. — "An Algorithm for the Machine Calculation of
+   Complex Fourier Series", *Mathematics of Computation*, 19(90), pp. 297-301, 1965.
 
 ---
 

diff --git a/include/numx/ntt.h b/include/numx/ntt.h
@@ -36,6 +36,9 @@
  * The output is 128 pairs (f[2i], f[2i+1]) representing f mod (x^2 - c_i).
  * All coefficients are reduced to [0, q-1] after the transform.
  *
+ * References: Cooley & Tukey, 1965; Longa & Naehrig, CANS 2016 (NTT-specific
+ * speedup). Full citations: docs/algorithms/ntt.md#references.
+ *
  * @param[in,out] f  Array of NUMX_NTT_N coefficients. Must not be NULL.
  *                   Each value must be in [0, NUMX_NTT_Q - 1] on entry.
  *
@@ -52,6 +55,10 @@ numx_status_t numx_ntt_forward(numx_q15_t *f);
  * normalization constant 128^{-1} = 3303 (mod 3329). The output coefficients
  * are reduced to [0, q-1].
  *
+ * References: Gentleman & Sande, AFIPS 1966; Avanzi et al. (Kyber spec),
+ * NIST PQC Round 3, 2021, for the normalization constant. Full citations:
+ * docs/algorithms/ntt.md#references.
+ *
  * @param[in,out] f  Array of NUMX_NTT_N NTT-domain coefficients. Must not be NULL.
  *
  * @return NUMX_OK on success.
@@ -72,6 +79,9 @@ numx_status_t numx_ntt_inverse(numx_q15_t *f);
  * Combined with numx_ntt_inverse, this implements polynomial multiplication
  * in Z_3329[x]/(x^256 + 1) in O(n) time after the O(n log n) NTT steps.
  *
+ * Reference: Avanzi et al. (Kyber spec), NIST PQC Round 3, 2021, degree-2
+ * ring basemul formula. Full citation: docs/algorithms/ntt.md#references.
+ *
  * @param[in]  a    First NTT-domain array of NUMX_NTT_N coefficients. Must not be NULL.
  * @param[in]  b    Second NTT-domain array of NUMX_NTT_N coefficients. Must not be NULL.
  * @param[out] out  Output array of NUMX_NTT_N coefficients. Must not be NULL.
@@ -98,6 +108,10 @@ numx_status_t numx_ntt_pointwise_mul(
  * All temporaries are stack-allocated (512 bytes). The caller must supply
  * distinct output buffer out; a and b may be the same pointer.
  *
+ * Cross-validated bit-for-bit against PQClean's ml-kem-512 reference
+ * implementation (FIPS 203), not just numx's own naive reference multiply.
+ * See validation/reference/ntt/README.md for provenance.
+ *
  * @param[in]  a    First polynomial: NUMX_NTT_N coefficients in [0, q-1]. Must not be NULL.
  * @param[in]  b    Second polynomial: NUMX_NTT_N coefficients in [0, q-1]. Must not be NULL.
  * @param[out] out  Product polynomial: NUMX_NTT_N coefficients in [0, q-1]. Must not be NULL.
@@ -117,6 +131,9 @@ numx_status_t numx_ntt_polymul(
  * values in [0, 2*q-1] to the canonical range [0, q-1]. Useful after
  * external coefficient arithmetic before passing to numx_ntt_forward.
  *
+ * Reference: Barrett, CRYPTO 1986. Full citation:
+ * docs/algorithms/ntt.md#references.
+ *
  * @param[in,out] f  Array of NUMX_NTT_N coefficients. Must not be NULL.
  *
  * @return NUMX_OK on success.

diff --git a/src/ntt.c b/src/ntt.c
@@ -428,6 +428,7 @@ static const int16_t priv_basemul[128] = {
 
 /*
  * Barrett reduction mod 3329 for inputs in [0, 2*q^2].
+ * Reference: Barrett, CRYPTO 1986 (docs/algorithms/ntt.md#references [4]).
  *
  * Branchless canonicalization: uses only unsigned arithmetic and a boolean
  * comparison (well-defined in C99), not the reference Kyber technique's

diff --git a/tests/test_ntt.c b/tests/test_ntt.c
@@ -12,6 +12,7 @@
 
 #include "unity.h"
 #include "numx/ntt.h"
+#include "vectors/ntt_pqclean_kat.h"
 
 /* q = 3329 */
 #define NTT_Q 3329
@@ -88,6 +89,33 @@ void test_ntt_forward_zero(void)
         TEST_ASSERT_EQUAL_INT16(0, f[i]);
 }
 
+/* L1 */
+void test_ntt_forward_vs_pqclean_reference(void)
+{
+    /* Cross-validated against PQClean's ml-kem-512 (FIPS 203) reference
+     * implementation, not numx's own code -- see
+     * validation/reference/ntt/README.md for provenance. */
+    static const numx_q15_t *const inputs[5] = {
+        ntt_pqclean_a_0, ntt_pqclean_a_1, ntt_pqclean_a_2,
+        ntt_pqclean_a_3, ntt_pqclean_a_4
+    };
+    static const numx_q15_t *const expected[5] = {
+        ntt_pqclean_ntt_a_0, ntt_pqclean_ntt_a_1, ntt_pqclean_ntt_a_2,
+        ntt_pqclean_ntt_a_3, ntt_pqclean_ntt_a_4
+    };
+    numx_q15_t f[256];
+    int v, i;
+
+    for (v = 0; v < 5; v++)
+    {
+        for (i = 0; i < 256; i++) f[i] = inputs[v][i];
+        TEST_ASSERT_EQUAL(NUMX_OK, numx_ntt_forward(f));
+        for (i = 0; i < 256; i++)
+            TEST_ASSERT_EQUAL_INT16_MESSAGE(expected[v][i], f[i],
+                                             "forward NTT != PQClean reference");
+    }
+}
+
 /* L4 */
 void test_ntt_forward_null(void)
 {
@@ -335,6 +363,37 @@ void test_ntt_polymul_random_vs_ref(void)
     }
 }
 
+/* L1 */
+void test_ntt_polymul_vs_pqclean_reference(void)
+{
+    /* Full negacyclic polynomial multiply cross-validated against
+     * PQClean's ml-kem-512 (FIPS 203) reference implementation, not
+     * numx's own naive poly_mul_ref() -- see
+     * validation/reference/ntt/README.md for provenance. */
+    static const numx_q15_t *const inputs_a[5] = {
+        ntt_pqclean_a_0, ntt_pqclean_a_1, ntt_pqclean_a_2,
+        ntt_pqclean_a_3, ntt_pqclean_a_4
+    };
+    static const numx_q15_t *const inputs_b[5] = {
+        ntt_pqclean_b_0, ntt_pqclean_b_1, ntt_pqclean_b_2,
+        ntt_pqclean_b_3, ntt_pqclean_b_4
+    };
+    static const numx_q15_t *const expected[5] = {
+        ntt_pqclean_polymul_0, ntt_pqclean_polymul_1, ntt_pqclean_polymul_2,
+        ntt_pqclean_polymul_3, ntt_pqclean_polymul_4
+    };
+    numx_q15_t out[256];
+    int v, i;
+
+    for (v = 0; v < 5; v++)
+    {
+        TEST_ASSERT_EQUAL(NUMX_OK, numx_ntt_polymul(inputs_a[v], inputs_b[v], out));
+        for (i = 0; i < 256; i++)
+            TEST_ASSERT_EQUAL_INT16_MESSAGE(expected[v][i], out[i],
+                                             "NTT polymul != PQClean reference");
+    }
+}
+
 /* L4 */
 void test_ntt_polymul_null(void)
 {
@@ -485,6 +544,7 @@ void numx_test_ntt(void)
     /* Forward NTT */
     RUN_TEST(test_ntt_forward_delta);
     RUN_TEST(test_ntt_forward_zero);
+    RUN_TEST(test_ntt_forward_vs_pqclean_reference);
     RUN_TEST(test_ntt_forward_null);
 
     /* Inverse NTT */
@@ -507,6 +567,7 @@ void numx_test_ntt(void)
     RUN_TEST(test_ntt_polymul_commutativity);
     RUN_TEST(test_ntt_polymul_known_linear);
     RUN_TEST(test_ntt_polymul_random_vs_ref);
+    RUN_TEST(test_ntt_polymul_vs_pqclean_reference);
     RUN_TEST(test_ntt_polymul_null);
 
     /* Reduce */