ring.go raw

   1  // Package ring implements polynomial arithmetic over R_q = Z_q[x]/(x^n + 1).
   2  //
   3  // This is the algebraic engine for the hamadryad cryptosystem. All lattice
   4  // operations — hash (SIS), signatures (GPV), encryption (LWE), homomorphic
   5  // evaluation — reduce to polynomial arithmetic in this ring.
   6  //
   7  // The ring R_q = Z_q[x]/(x^n + 1) is a quotient of the polynomial ring Z_q[x]
   8  // by the cyclotomic polynomial Phi_{2n}(x) = x^n + 1. When n is a power of 2
   9  // and q ≡ 1 (mod 2n), this ring supports a negacyclic NTT, enabling O(n log n)
  10  // polynomial multiplication.
  11  //
  12  // Parameters are selected for cryptographic use:
  13  //
  14  //	Falcon-512:    n=512,  q=12289       (signature scheme, ~128-bit security)
  15  //	Kyber-512:     n=256,  q=3329        (KEM, ~128-bit security)
  16  //	Dilithium-II:  n=256,  q=8380417     (signature, ~128-bit security)
  17  //
  18  // The default parameter set targets Falcon-512 (n=512, q=12289) because:
  19  //   - Hash-and-sign naturally extends hamadryad (both are SIS-based)
  20  //   - n=512 gives 128-bit post-quantum security
  21  //   - q=12289 is NTT-friendly: 12289 ≡ 1 (mod 1024), primitive 1024th root exists
  22  //   - All coefficients fit in uint16 (q < 2^14), enabling fast arithmetic
  23  //
  24  // Security: all constructions reduce to SVP on ideal lattices in R = Z[x]/(x^n+1).
  25  package ring
  26  
  27  // Params defines the ring R_q = Z_q[x]/(x^n + 1).
  28  type Params struct {
  29  	// N is the ring dimension (degree of cyclotomic polynomial).
  30  	// Must be a power of 2.
  31  	N int
  32  
  33  	// Q is the coefficient modulus. Must be prime with Q ≡ 1 (mod 2N)
  34  	// to guarantee existence of primitive 2N-th roots of unity for NTT.
  35  	Q uint32
  36  
  37  	// RootOfUnity is a primitive 2N-th root of unity mod Q.
  38  	// Used for the negacyclic NTT: psi such that psi^(2N) ≡ 1 (mod Q)
  39  	// and psi^N ≡ -1 (mod Q).
  40  	RootOfUnity uint32
  41  
  42  	// MontR is the Montgomery parameter R = 2^k where k is chosen
  43  	// so that R > Q. For Q < 2^16, k=16. For Q < 2^32, k=32.
  44  	MontR uint64
  45  
  46  	// QInv is -Q^{-1} mod R (Montgomery reduction constant).
  47  	QInv uint64
  48  }
  49  
  50  // Falcon512 returns parameters for the Falcon-512 ring.
  51  //
  52  //	R_q = Z_12289[x]/(x^512 + 1)
  53  //	12289 = 12·1024 + 1 = 3·2^12 + 1
  54  //	Primitive 1024th root of unity: 49 (i.e., 49^1024 ≡ 1, 49^512 ≡ -1 mod 12289)
  55  func Falcon512() Params {
  56  	return Params{
  57  		N:           512,
  58  		Q:           12289,
  59  		RootOfUnity: 49,
  60  		MontR:       1 << 16,
  61  		QInv:        qinv(12289, 16),
  62  	}
  63  }
  64  
  65  // NewHope256 returns parameters for a NewHope-style ring with n=256.
  66  //
  67  //	R_q = Z_7681[x]/(x^256 + 1)
  68  //	7681 = 15·512 + 1, so 512 | 7680 = q-1
  69  //	Primitive 512th root of unity: 4055
  70  //
  71  // Suitable for KEM (Ring-LWE encryption) at ~128-bit post-quantum security.
  72  // Unlike Kyber's q=3329 which requires an incomplete NTT (512 ∤ 3328),
  73  // q=7681 supports the full negacyclic NTT.
  74  func NewHope256() Params {
  75  	return Params{
  76  		N:           256,
  77  		Q:           7681,
  78  		RootOfUnity: 4055,
  79  		MontR:       1 << 16,
  80  		QInv:        qinv(7681, 16),
  81  	}
  82  }
  83  
  84  // HE64 returns parameters for homomorphic evaluation with n=64.
  85  //
  86  //	R_q = Z_10000769[x]/(x^64 + 1)
  87  //	q = 10000769 (24-bit prime, q ≡ 1 mod 128)
  88  //	Primitive 128th root of unity: 6028202
  89  //
  90  // n=64 matches the hamadryad ring dimension. The 24-bit q gives enough
  91  // noise budget for depth-1 multiplication:
  92  //   Fresh noise ≈ 2 * n * eta = 256 (eta=2)
  93  //   Tensor product noise ≈ 4 * n * 256^2 ≈ 16.8M < q/2 ≈ 5M
  94  //
  95  // Wait — that's still too tight. Let's check:
  96  //   With eta=1: noise ≈ 2 * 64 * 1 = 128
  97  //   Product: 4 * 64 * 128^2 = 4,194,304 < 5,000,384 (q/2)  ✓
  98  //
  99  // So eta=1 gives exactly one level of multiplication.
 100  func HE64() Params {
 101  	return Params{
 102  		N:           64,
 103  		Q:           10000769,
 104  		RootOfUnity: 6028202,
 105  		MontR:       1 << 32,
 106  		QInv:        qinv(10000769, 32),
 107  	}
 108  }
 109  
 110  // Falcon1024 returns parameters for the Falcon-1024 ring (~256-bit security).
 111  //
 112  //	R_q = Z_12289[x]/(x^1024 + 1)
 113  //	12289-1 = 12288 = 6·2048, so 2048 | 12288
 114  //	Primitive 2048th root of unity: 1945
 115  func Falcon1024() Params {
 116  	return Params{
 117  		N:           1024,
 118  		Q:           12289,
 119  		RootOfUnity: 1945,
 120  		MontR:       1 << 16,
 121  		QInv:        qinv(12289, 16),
 122  	}
 123  }
 124  
 125  // qinv computes -Q^{-1} mod 2^k using Hensel lifting.
 126  // Starts from Q^{-1} mod 2 = 1 (Q is odd) and doubles precision each step.
 127  func qinv(q uint32, k int) uint64 {
 128  	// Q is odd, so Q^{-1} mod 2 = 1.
 129  	inv := uint64(1)
 130  	for i := 1; i < k; i++ {
 131  		// Hensel lift: if inv·Q ≡ 1 (mod 2^i), then
 132  		// inv' = inv·(2 - inv·Q) satisfies inv'·Q ≡ 1 (mod 2^{2i}).
 133  		inv = inv * (2 - inv*uint64(q))
 134  	}
 135  	// We want -Q^{-1} mod 2^k.
 136  	mask := (uint64(1) << k) - 1
 137  	return (-inv) & mask
 138  }
 139  
 140  // Poly is a polynomial in R_q, stored as N coefficients in [0, Q).
 141  // Can be in either coefficient form or NTT (evaluation) form.
 142  // Operations document which form they expect.
 143  type Poly struct {
 144  	Coeffs []uint32
 145  	params Params
 146  	isNTT  bool
 147  }
 148  
 149  // New creates a zero polynomial in coefficient form.
 150  func New(p Params) *Poly {
 151  	return &Poly{
 152  		Coeffs: make([]uint32, p.N),
 153  		params: p,
 154  	}
 155  }
 156  
 157  // NewFromCoeffs creates a polynomial from a coefficient slice.
 158  // Coefficients are reduced mod Q.
 159  func NewFromCoeffs(p Params, coeffs []uint32) *Poly {
 160  	poly := New(p)
 161  	n := min(len(coeffs), p.N)
 162  	for i := 0; i < n; i++ {
 163  		poly.Coeffs[i] = coeffs[i] % p.Q
 164  	}
 165  	return poly
 166  }
 167  
 168  // Clone returns a deep copy.
 169  func (a *Poly) Clone() *Poly {
 170  	b := &Poly{
 171  		Coeffs: make([]uint32, len(a.Coeffs)),
 172  		params: a.params,
 173  		isNTT:  a.isNTT,
 174  	}
 175  	copy(b.Coeffs, a.Coeffs)
 176  	return b
 177  }
 178  
 179  // Params returns the ring parameters.
 180  func (a *Poly) Params() Params { return a.params }
 181  
 182  // IsNTT reports whether the polynomial is in NTT form.
 183  func (a *Poly) IsNTT() bool { return a.isNTT }
 184  
 185  // addMod computes (a + b) mod q without overflow for uint32 inputs.
 186  func addMod(a, b, q uint32) uint32 {
 187  	sum := a + b
 188  	if sum >= q {
 189  		sum -= q
 190  	}
 191  	return sum
 192  }
 193  
 194  // subMod computes (a - b) mod q without overflow.
 195  func subMod(a, b, q uint32) uint32 {
 196  	if a >= b {
 197  		return a - b
 198  	}
 199  	return q - b + a
 200  }
 201  
 202  // mulMod computes (a * b) mod q using uint64 intermediate.
 203  func mulMod(a, b, q uint32) uint32 {
 204  	return uint32((uint64(a) * uint64(b)) % uint64(q))
 205  }
 206  
 207  // Add sets c = a + b (coefficient-wise mod Q). a and b must be same form.
 208  func Add(a, b *Poly) *Poly {
 209  	c := New(a.params)
 210  	c.isNTT = a.isNTT
 211  	q := a.params.Q
 212  	for i := range a.Coeffs {
 213  		c.Coeffs[i] = addMod(a.Coeffs[i], b.Coeffs[i], q)
 214  	}
 215  	return c
 216  }
 217  
 218  // Sub sets c = a - b (coefficient-wise mod Q).
 219  func Sub(a, b *Poly) *Poly {
 220  	c := New(a.params)
 221  	c.isNTT = a.isNTT
 222  	q := a.params.Q
 223  	for i := range a.Coeffs {
 224  		c.Coeffs[i] = subMod(a.Coeffs[i], b.Coeffs[i], q)
 225  	}
 226  	return c
 227  }
 228  
 229  // Neg sets c = -a (coefficient-wise mod Q).
 230  func Neg(a *Poly) *Poly {
 231  	c := New(a.params)
 232  	c.isNTT = a.isNTT
 233  	q := a.params.Q
 234  	for i := range a.Coeffs {
 235  		if a.Coeffs[i] != 0 {
 236  			c.Coeffs[i] = q - a.Coeffs[i]
 237  		}
 238  	}
 239  	return c
 240  }
 241  
 242  // MulPointwise computes c = a * b pointwise. Both must be in NTT form.
 243  func MulPointwise(a, b *Poly) *Poly {
 244  	c := New(a.params)
 245  	c.isNTT = true
 246  	q := a.params.Q
 247  	for i := range a.Coeffs {
 248  		c.Coeffs[i] = mulMod(a.Coeffs[i], b.Coeffs[i], q)
 249  	}
 250  	return c
 251  }
 252  
 253  // MulAdd computes c += a * b pointwise (fused multiply-accumulate in NTT domain).
 254  func MulAdd(c, a, b *Poly) {
 255  	q := a.params.Q
 256  	for i := range a.Coeffs {
 257  		c.Coeffs[i] = addMod(c.Coeffs[i], mulMod(a.Coeffs[i], b.Coeffs[i], q), q)
 258  	}
 259  }
 260  
 261  // ScalarMul computes c = s * a (each coefficient multiplied by scalar s mod Q).
 262  func ScalarMul(a *Poly, s uint32) *Poly {
 263  	c := New(a.params)
 264  	c.isNTT = a.isNTT
 265  	q := a.params.Q
 266  	s = s % q
 267  	for i := range a.Coeffs {
 268  		c.Coeffs[i] = mulMod(a.Coeffs[i], s, q)
 269  	}
 270  	return c
 271  }
 272  
 273  // Equal reports whether two polynomials have identical coefficients.
 274  func Equal(a, b *Poly) bool {
 275  	if len(a.Coeffs) != len(b.Coeffs) {
 276  		return false
 277  	}
 278  	for i := range a.Coeffs {
 279  		if a.Coeffs[i] != b.Coeffs[i] {
 280  			return false
 281  		}
 282  	}
 283  	return true
 284  }
 285