Explain the bound 3q/4 which applies to all montmul-by-const

To ease the explanation, we cite and add
[Survey_Hwang23](https://eprint.iacr.org/2023/1962) to the bibliography.

In the process, we noticed that the original MLD_INTT_BOUND defined as
(MLDSA_Q * 3 / 4) = floor(3q/4) was too strong to justify concisely, so
we bump it to (MLDSA_Q * 3 / 4 + 1) or ceil(3q/4).

Signed-off-by: jammychiou1 <[email protected]>

Author: jammychiou1

BIBLIOGRAPHY.md

@@ -247,6 +247,15 @@ source code and documentation.
   - [mldsa/src/native/aarch64/src/intt.S](mldsa/src/native/aarch64/src/intt.S)
   - [mldsa/src/native/aarch64/src/ntt.S](mldsa/src/native/aarch64/src/ntt.S)
 
+### `Survey_Hwang23`
+
+* A Survey of Polynomial Multiplications for Lattice-Based Cryptosystems
+* Author(s):
+  - Vincent Hwang
+* URL: https://eprint.iacr.org/2023/1962
+* Referenced from:
+  - [dev/x86_64/src/intt.S](dev/x86_64/src/intt.S)
+
 ### `libmceliece`
 
 * libmceliece implementation of Classic McEliece

BIBLIOGRAPHY.yml

@@ -103,6 +103,14 @@
   - Kannwischer, Matthias J.
   - Yang, Bo-Yin
   - Yang, Shang-Yi
+
+- id: Survey_Hwang23
+  name: "A Survey of Polynomial Multiplications for Lattice-Based Cryptosystems"
+  year: 2023
+  url: https://eprint.iacr.org/2023/1962
+  author:
+  - Hwang, Vincent
+
 - id: mupq
   name: Common files for pqm4, pqm3, pqriscv
   author: 

dev/x86_64/src/intt.S

@@ -15,6 +15,11 @@
  *   CRYSTALS-Dilithium optimized AVX2 implementation
  *   Bai, Ducas, Kiltz, Lepoint, Lyubashevsky, Schwabe, Seiler, Stehlé
  *   https://github.com/pq-crystals/dilithium/tree/master/avx2
+ *
+ * - [Survey_Hwang23]
+ *   A Survey of Polynomial Multiplications for Lattice-Based Cryptosystems
+ *   Vincent Hwang
+ *   https://eprint.iacr.org/2023/1962
  */
 
 #include "../../../common.h"
@@ -47,7 +52,8 @@ vpblendd	$0xAA,%ymm\r1,%ymm\r0,%ymm\r3
  * Compute l + h, montmul(h - l, zh) then store the results back to l, h
  * respectively.
  *
- * The general abs bound of Montgomery multiplication is 3q/4.
+ * The abs bound of "Montgomery multiplication with signed canonical constant"
+ * is ceil(3q/4) (see the end of this file).
  */
 .macro butterfly l,h,zl0=1,zl1=1,zh0=2,zh1=2
 vpsubd		%ymm\l,%ymm\h,%ymm12
@@ -107,7 +113,7 @@ vmovshdup	%ymm3,%ymm1
 vmovshdup	%ymm15,%ymm2
 butterfly	10,11,1,3,2,15
 
-/* 4, 6, 8, 10: abs bound < 2q; 5, 7, 9, 11: abs bound < 3q/4 */
+/* 4, 6, 8, 10: abs bound < 2q; 5, 7, 9, 11: abs bound < ceil(3q/4) */
 /*
  * Note that since 2^31 / q > 256, the sum of all 256 coefficients does not
  * overflow. This allows us to greatly simplify the range analysis by relaxing
@@ -236,7 +242,7 @@ butterfly	5,9
 butterfly	6,10
 butterfly	7,11
 
-/* 4, 5, 6, 7: abs bound < 256q; 8, 9, 10, 11: abs bound < 3q/4 */
+/* 4, 5, 6, 7: abs bound < 256q; 8, 9, 10, 11: abs bound < ceil(3q/4) */
 
 vmovdqa         %ymm8,512+32*\off(%rdi)
 vmovdqa         %ymm9,640+32*\off(%rdi)
@@ -251,8 +257,7 @@ vmovdqa         %ymm11,896+32*\off(%rdi)
  * For ymm{8,9,10,11}, the scaling has been merged into the last butterfly, so
  * only ymm{4,5,6,7} need to be scaled explicitly.
  *
- * The scaling is achieved by computing montmul(-, MLD_AVX2_DIV), so the output
- * will have an abs bound of 3q/4.
+ * The scaling is achieved by computing montmul(-, MLD_AVX2_DIV).
  *
  * 4, 5, 6, 7: abs bound < 256q
  */
@@ -305,7 +310,22 @@ vmovshdup	%ymm7,%ymm7
 vpblendd	$0xAA,%ymm8,%ymm6,%ymm6
 vpblendd	$0xAA,%ymm9,%ymm7,%ymm7
 
-/* 4, 5, 6, 7: abs bound < 3q/4 */
+/*
+ * The bound ceil(3q/4) for this scaling, as well as any other "Montgomery
+ * multiplication with signed canonical constant", is justified as follows.
+ *
+ * In @[Survey_Hwang23, Section 2.2] they showed a bound that works for any
+ * variable input a, as long as the constant b is signed canonical:
+ *
+ *   |montmul(a, b)| <= (|a| (q/2) + (R/2) q) / R = (q/2) (1 + |a|/R).
+ *
+ * Therefore, even if we know nothing about a except that it fits inside
+ * int32_t (thus |a| <= R/2), we still have |montmul(a, b)| <= 3q/4. This can be
+ * strengthened to |montmul_pos(a, b)| <= floor(3q/4) < ceil(3q/4) since LHS is
+ * an integer and 3q/4 isn't.
+ */
+
+/* 4, 5, 6, 7: abs bound < ceil(3q/4) */
 
 vmovdqa         %ymm4,  0+32*\off(%rdi)
 vmovdqa         %ymm5,128+32*\off(%rdi)

dev/x86_64/src/ntt.S

@@ -48,8 +48,9 @@ vpblendd	$0xAA,%ymm\r1,%ymm\r0,%ymm\r3
  * Compute l + montmul(h, zh), l - montmul(h, zh) then store the results back to
  * l, h respectively.
  *
- * Although the general abs bound of Montgomery multiplication is 3q/4, we use
- * the more convenient bound q here.
+ * Although the abs bound of "Montgomery multiplication with signed canonical
+ * constant" is ceil(3q/4) (see the end of dev/x86_64/src/intt.S), we use the
+ * more convenient bound q here.
  *
  * In conclusion, the magnitudes of all coefficients grow by at most q after
  * each layer.

mldsa/src/native/x86_64/src/intt.S

@@ -15,6 +15,11 @@
  *   CRYSTALS-Dilithium optimized AVX2 implementation
  *   Bai, Ducas, Kiltz, Lepoint, Lyubashevsky, Schwabe, Seiler, Stehlé
  *   https://github.com/pq-crystals/dilithium/tree/master/avx2
+ *
+ * - [Survey_Hwang23]
+ *   A Survey of Polynomial Multiplications for Lattice-Based Cryptosystems
+ *   Vincent Hwang
+ *   https://eprint.iacr.org/2023/1962
  */
 
 #include "../../../common.h"

mldsa/src/ntt.h

@@ -21,8 +21,8 @@
 
 /* Absolute exclusive upper bound for the output of the forward NTT */
 #define MLD_NTT_BOUND (9 * MLDSA_Q)
-/* Absolute exclusive upper bound for the output of the inverse NTT*/
-#define MLD_INTT_BOUND (MLDSA_Q * 3 / 4)
+/* Absolute exclusive upper bound for the output of the inverse NTT */
+#define MLD_INTT_BOUND (MLDSA_Q * 3 / 4 + 1) /* ceil(3 * MLDSA_Q / 4) */
 
 #define mld_ntt MLD_NAMESPACE(ntt)
 /*************************************************