mirror of
https://github.com/luxfi/math.git
synced 2026-07-27 03:38:49 +00:00
merge LP-107 Phase 3: NTT canonical body migration
This commit is contained in:
Hanzo AI
@@ -4,7 +4,6 @@ go 1.26.1
|
||||
|
||||
require (
|
||||
github.com/luxfi/ids v1.2.7
|
||||
github.com/luxfi/lattice/v7 v7.1.0
|
||||
github.com/luxfi/math/big v0.1.0
|
||||
github.com/luxfi/sampler v1.0.0
|
||||
github.com/stretchr/testify v1.11.1
|
||||
@@ -13,16 +12,15 @@ require (
|
||||
)
|
||||
|
||||
require (
|
||||
github.com/ALTree/bigfloat v0.2.0 // indirect
|
||||
github.com/davecgh/go-spew v1.1.2-0.20180830191138-d8f796af33cc // indirect
|
||||
github.com/google/go-cmp v0.7.0 // indirect
|
||||
github.com/klauspost/cpuid/v2 v2.3.0 // indirect
|
||||
github.com/kr/text v0.2.0 // indirect
|
||||
github.com/kr/pretty v0.3.1 // indirect
|
||||
github.com/luxfi/crypto v1.17.36 // indirect
|
||||
github.com/mr-tron/base58 v1.2.0 // indirect
|
||||
github.com/pmezard/go-difflib v1.0.1-0.20181226105442-5d4384ee4fb2 // indirect
|
||||
golang.org/x/crypto v0.49.0 // indirect
|
||||
golang.org/x/sys v0.42.0 // indirect
|
||||
gonum.org/v1/gonum v0.17.0 // indirect
|
||||
gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c // indirect
|
||||
gopkg.in/yaml.v3 v3.0.1 // indirect
|
||||
)
|
||||
|
||||
@@ -1,30 +1,29 @@
|
||||
github.com/ALTree/bigfloat v0.2.0 h1:AwNzawrpFuw55/YDVlcPw0F0cmmXrmngBHhVrvdXPvM=
|
||||
github.com/ALTree/bigfloat v0.2.0/go.mod h1:+NaH2gLeY6RPBPPQf4aRotPPStg+eXc8f9ZaE4vRfD4=
|
||||
github.com/creack/pty v1.1.9/go.mod h1:oKZEueFk5CKHvIhNR5MUki03XCEU+Q6VDXinZuGJ33E=
|
||||
github.com/davecgh/go-spew v1.1.2-0.20180830191138-d8f796af33cc h1:U9qPSI2PIWSS1VwoXQT9A3Wy9MM3WgvqSxFWenqJduM=
|
||||
github.com/davecgh/go-spew v1.1.2-0.20180830191138-d8f796af33cc/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=
|
||||
github.com/google/go-cmp v0.7.0 h1:wk8382ETsv4JYUZwIsn6YpYiWiBsYLSJiTsyBybVuN8=
|
||||
github.com/google/go-cmp v0.7.0/go.mod h1:pXiqmnSA92OHEEa9HXL2W4E7lf9JzCmGVUdgjX3N/iU=
|
||||
github.com/klauspost/cpuid/v2 v2.3.0 h1:S4CRMLnYUhGeDFDqkGriYKdfoFlDnMtqTiI/sFzhA9Y=
|
||||
github.com/klauspost/cpuid/v2 v2.3.0/go.mod h1:hqwkgyIinND0mEev00jJYCxPNVRVXFQeu1XKlok6oO0=
|
||||
github.com/kr/pretty v0.2.1/go.mod h1:ipq/a2n7PKx3OHsz4KJII5eveXtPO4qwEXGdVfWzfnI=
|
||||
github.com/kr/pretty v0.3.1 h1:flRD4NNwYAUpkphVc1HcthR4KEIFJ65n8Mw5qdRn3LE=
|
||||
github.com/kr/pretty v0.3.1/go.mod h1:hoEshYVHaxMs3cyo3Yncou5ZscifuDolrwPKZanG3xk=
|
||||
github.com/kr/pty v1.1.1/go.mod h1:pFQYn66WHrOpPYNljwOMqo10TkYh1fy3cYio2l3bCsQ=
|
||||
github.com/kr/text v0.1.0/go.mod h1:4Jbv+DJW3UT/LiOwJeYQe1efqtUx/iVham/4vfdArNI=
|
||||
github.com/kr/text v0.2.0 h1:5Nx0Ya0ZqY2ygV366QzturHI13Jq95ApcVaJBhpS+AY=
|
||||
github.com/kr/text v0.2.0/go.mod h1:eLer722TekiGuMkidMxC/pM04lWEeraHUUmBw8l2grE=
|
||||
github.com/luxfi/crypto v1.17.36 h1:fN/m5oF2kd44/5JXSILnBxjm3AmgCTw8qFwfkdJaAxg=
|
||||
github.com/luxfi/crypto v1.17.36/go.mod h1:UHzltsXPDF/tGcmTiL/DbkClAflr3YmvOaXlwWC6eXo=
|
||||
github.com/luxfi/ids v1.2.7 h1:3B42EbzR2cdY3veo1yOFOOIeEE+HULnCye2Ye32nXqI=
|
||||
github.com/luxfi/ids v1.2.7/go.mod h1:svLsj7e6ixJVfRYaQqv9RWjBIW11Vz738+d08BVpeCg=
|
||||
github.com/luxfi/lattice/v7 v7.1.0 h1:mr3HvN6olNTS2LT/xAW/JBhTqfvpsGmsopDMeR7BSJs=
|
||||
github.com/luxfi/lattice/v7 v7.1.0/go.mod h1:IaaUN+3ysnBG4BA8ILRYG0j80+qtYDP4C5lkaDb2pDE=
|
||||
github.com/luxfi/math/big v0.1.0 h1:Vz4c0RsZVPdIKPsHPgAJChH/R3p15WHRUz7LkLf+NIQ=
|
||||
github.com/luxfi/math/big v0.1.0/go.mod h1:BuxSu22RbO93xBLk5Eam5nldFponoJ73xDFz4uJ3Huk=
|
||||
github.com/luxfi/sampler v1.0.0 h1:k8Sf6otW83w4pQp0jXLA+g3J/joB7w7SqXQsWmNTOV0=
|
||||
github.com/luxfi/sampler v1.0.0/go.mod h1:f96/ozlj9vFfZj+akLtrHn4VpulQahwB+MQQhpeIekk=
|
||||
github.com/mr-tron/base58 v1.2.0 h1:T/HDJBh4ZCPbU39/+c3rRvE0uKBQlU27+QI8LJ4t64o=
|
||||
github.com/mr-tron/base58 v1.2.0/go.mod h1:BinMc/sQntlIE1frQmRFPUoPA1Zkr8VRgBdjWI2mNwc=
|
||||
github.com/pkg/diff v0.0.0-20210226163009-20ebb0f2a09e/go.mod h1:pJLUxLENpZxwdsKMEsNbx1VGcRFpLqf3715MtcvvzbA=
|
||||
github.com/pmezard/go-difflib v1.0.1-0.20181226105442-5d4384ee4fb2 h1:Jamvg5psRIccs7FGNTlIRMkT8wgtp5eCXdBlqhYGL6U=
|
||||
github.com/pmezard/go-difflib v1.0.1-0.20181226105442-5d4384ee4fb2/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4=
|
||||
github.com/rogpeppe/go-internal v1.9.0/go.mod h1:WtVeX8xhTBvf0smdhujwtBcq4Qrzq/fJaraNFVN+nFs=
|
||||
github.com/rogpeppe/go-internal v1.14.1 h1:UQB4HGPB6osV0SQTLymcB4TgvyWu6ZyliaW0tI/otEQ=
|
||||
github.com/rogpeppe/go-internal v1.14.1/go.mod h1:MaRKkUm5W0goXpeCfT7UZI6fk/L7L7so1lCWt35ZSgc=
|
||||
github.com/stretchr/testify v1.11.1 h1:7s2iGBzp5EwR7/aIZr8ao5+dra3wiQyKjjFuvgVKu7U=
|
||||
|
||||
@@ -0,0 +1,23 @@
|
||||
// Copyright (c) 2026 Lux Industries Inc.
|
||||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
// Package canonical owns the canonical Lattigo-derived Montgomery NTT
|
||||
// body — SubRing layout, root-table generator, scalar primitives
|
||||
// (MRed/MRedLazy/MForm/BRed/CRed/...), the loop-unrolled-by-16 NTT
|
||||
// kernel, and the SIMD AVX2 dispatch hooks.
|
||||
//
|
||||
// LP-107 Phase 3: this package is the single source of truth for the
|
||||
// Montgomery NTT body. luxfi/lattice/v7/ring re-exports it as a thin
|
||||
// shim so its downstream consumers (luxfi/pulsar, luxfi/fhe,
|
||||
// luxfi/threshold, ...) compile unchanged. Application callers should
|
||||
// route through one of:
|
||||
//
|
||||
// - github.com/luxfi/math/ntt — the canonical Service /
|
||||
// Backend interface (recommended for new code).
|
||||
// - github.com/luxfi/lattice/v7/ring — the historical SubRing API
|
||||
// for FHE-flavoured callers.
|
||||
//
|
||||
// This package is NOT a general-purpose import target: its surface is
|
||||
// optimized for the two consumers above and is not part of the public
|
||||
// luxfi/math semver contract for downstream code.
|
||||
package canonical
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,188 @@
|
||||
// Copyright (c) 2026 Lux Industries Inc.
|
||||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
//
|
||||
// Body migrated from github.com/luxfi/lattice/v7/ring/modular_reduction.go
|
||||
// (Lattigo-derived). LP-107 Phase 3 establishes luxfi/math as the
|
||||
// canonical owner of the Montgomery / Barrett scalar primitives; the
|
||||
// lattice copy now delegates to this package.
|
||||
|
||||
package canonical
|
||||
|
||||
import (
|
||||
"math/big"
|
||||
"math/bits"
|
||||
)
|
||||
|
||||
// MForm switches a to the Montgomery domain by computing
|
||||
// a*2^64 mod q.
|
||||
func MForm(a, q uint64, bredconstant [2]uint64) (r uint64) {
|
||||
mhi, _ := bits.Mul64(a, bredconstant[1])
|
||||
r = -(a*bredconstant[0] + mhi) * q
|
||||
if r >= q {
|
||||
r -= q
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// MFormLazy switches a to the Montgomery domain by computing
|
||||
// a*2^64 mod q in constant time.
|
||||
// The result is between 0 and 2*q-1.
|
||||
func MFormLazy(a, q uint64, bredconstant [2]uint64) (r uint64) {
|
||||
mhi, _ := bits.Mul64(a, bredconstant[1])
|
||||
r = -(a*bredconstant[0] + mhi) * q
|
||||
return
|
||||
}
|
||||
|
||||
// IMForm switches a from the Montgomery domain back to the
|
||||
// standard domain by computing a*(1/2^64) mod q.
|
||||
func IMForm(a, q, mredconstant uint64) (r uint64) {
|
||||
r, _ = bits.Mul64(a*mredconstant, q)
|
||||
r = q - r
|
||||
if r >= q {
|
||||
r -= q
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// IMFormLazy switches a from the Montgomery domain back to the
|
||||
// standard domain by computing a*(1/2^64) mod q in constant time.
|
||||
// The result is between 0 and 2*q-1.
|
||||
func IMFormLazy(a, q, mredconstant uint64) (r uint64) {
|
||||
r, _ = bits.Mul64(a*mredconstant, q)
|
||||
r = q - r
|
||||
return
|
||||
}
|
||||
|
||||
// GenMRedConstant computes the constant mredconstant = (q^-1) mod 2^64 required for MRed.
|
||||
func GenMRedConstant(q uint64) (mredconstant uint64) {
|
||||
mredconstant = 1
|
||||
for i := 0; i < 63; i++ {
|
||||
mredconstant *= q
|
||||
q *= q
|
||||
}
|
||||
return mredconstant
|
||||
}
|
||||
|
||||
// MRed computes x * y * (1/2^64) mod q.
|
||||
func MRed(x, y, q, mredconstant uint64) (r uint64) {
|
||||
mhi, mlo := bits.Mul64(x, y)
|
||||
hhi, _ := bits.Mul64(mlo*mredconstant, q)
|
||||
r = mhi - hhi + q
|
||||
if r >= q {
|
||||
r -= q
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// MRedLazy computes x * y * (1/2^64) mod q in constant time.
|
||||
// The result is between 0 and 2*q-1.
|
||||
func MRedLazy(x, y, q, mredconstant uint64) (r uint64) {
|
||||
ahi, alo := bits.Mul64(x, y)
|
||||
H, _ := bits.Mul64(alo*mredconstant, q)
|
||||
r = ahi - H + q
|
||||
return
|
||||
}
|
||||
|
||||
// GenBRedConstant computes the constant for the BRed algorithm.
|
||||
// Returns ((2^128)/q)/(2^64) and (2^128)/q mod 2^64.
|
||||
func GenBRedConstant(q uint64) [2]uint64 {
|
||||
bigR, _ := new(big.Int).SetString("100000000000000000000000000000000", 16)
|
||||
bigR.Quo(bigR, new(big.Int).SetUint64(q))
|
||||
|
||||
mlo := bigR.Uint64()
|
||||
mhi := bigR.Rsh(bigR, 64).Uint64()
|
||||
|
||||
return [2]uint64{mhi, mlo}
|
||||
}
|
||||
|
||||
// BRedAdd computes a mod q.
|
||||
func BRedAdd(a, q uint64, bredconstant [2]uint64) (r uint64) {
|
||||
mhi, _ := bits.Mul64(a, bredconstant[0])
|
||||
r = a - mhi*q
|
||||
if r >= q {
|
||||
r -= q
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// BRedAddLazy computes a mod q in constant time.
|
||||
// The result is between 0 and 2*q-1.
|
||||
func BRedAddLazy(x, q uint64, bredconstant [2]uint64) uint64 {
|
||||
s0, _ := bits.Mul64(x, bredconstant[0])
|
||||
return x - s0*q
|
||||
}
|
||||
|
||||
// BRed computes x*y mod q.
|
||||
func BRed(x, y, q uint64, bredconstant [2]uint64) (r uint64) {
|
||||
var mhi, mlo, lhi, hhi, hlo, s0, carry uint64
|
||||
|
||||
mhi, mlo = bits.Mul64(x, y)
|
||||
|
||||
r = mhi * bredconstant[0]
|
||||
|
||||
hhi, hlo = bits.Mul64(mlo, bredconstant[0])
|
||||
|
||||
r += hhi
|
||||
|
||||
lhi, _ = bits.Mul64(mlo, bredconstant[1])
|
||||
|
||||
s0, carry = bits.Add64(hlo, lhi, 0)
|
||||
|
||||
r += carry
|
||||
|
||||
hhi, hlo = bits.Mul64(mhi, bredconstant[1])
|
||||
|
||||
r += hhi
|
||||
|
||||
_, carry = bits.Add64(hlo, s0, 0)
|
||||
|
||||
r += carry
|
||||
|
||||
r = mlo - r*q
|
||||
|
||||
if r >= q {
|
||||
r -= q
|
||||
}
|
||||
|
||||
return
|
||||
}
|
||||
|
||||
// BRedLazy computes x*y mod q in constant time.
|
||||
// The result is between 0 and 2*q-1.
|
||||
func BRedLazy(x, y, q uint64, bredconstant [2]uint64) (r uint64) {
|
||||
var mhi, mlo, lhi, hhi, hlo, s0, carry uint64
|
||||
|
||||
mhi, mlo = bits.Mul64(x, y)
|
||||
|
||||
r = mhi * bredconstant[0]
|
||||
|
||||
hhi, hlo = bits.Mul64(mlo, bredconstant[0])
|
||||
|
||||
r += hhi
|
||||
|
||||
lhi, _ = bits.Mul64(mlo, bredconstant[1])
|
||||
|
||||
s0, carry = bits.Add64(hlo, lhi, 0)
|
||||
|
||||
r += carry
|
||||
|
||||
hhi, hlo = bits.Mul64(mhi, bredconstant[1])
|
||||
|
||||
r += hhi
|
||||
|
||||
_, carry = bits.Add64(hlo, s0, 0)
|
||||
|
||||
r += carry
|
||||
|
||||
r = mlo - r*q
|
||||
|
||||
return
|
||||
}
|
||||
|
||||
// CRed reduce returns a mod q where a is between 0 and 2*q-1.
|
||||
func CRed(a, q uint64) uint64 {
|
||||
if a >= q {
|
||||
return a - q
|
||||
}
|
||||
return a
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,405 @@
|
||||
//go:build goexperiment.simd && amd64
|
||||
|
||||
// SIMD-accelerated NTT body for AVX2 amd64 with GOEXPERIMENT=simd.
|
||||
// LP-107 Phase 3 — verbatim from luxfi/lattice/v7/ring/ntt_simd.go.
|
||||
|
||||
package canonical
|
||||
|
||||
import (
|
||||
"math/bits"
|
||||
"simd/archsimd"
|
||||
"unsafe"
|
||||
)
|
||||
|
||||
// simdAVX2 reports whether the CPU supports AVX2 (minimum for SIMD NTT).
|
||||
var simdAVX2 = archsimd.X86.AVX2()
|
||||
|
||||
// mredLazy4 computes 4 parallel MRedLazy operations using AVX2 VPMULUDQ.
|
||||
//
|
||||
// MRedLazy(x, y, Q, MRedConstant) computes:
|
||||
//
|
||||
// ahi, alo = x * y (128-bit multiply)
|
||||
// H, _ = (alo * MRedConstant) * Q (128-bit multiply, high 64 bits)
|
||||
// result = ahi - H + Q
|
||||
//
|
||||
// 64x64->128 bit multiply is decomposed into 32-bit widening multiplies:
|
||||
//
|
||||
// a*b = a_lo*b_lo + (a_lo*b_hi + a_hi*b_lo)<<32 + a_hi*b_hi<<64
|
||||
//
|
||||
// We use Uint32x8.MulEvenWiden (VPMULUDQ, AVX2) which multiplies even-indexed
|
||||
// 32-bit elements producing 64-bit results. By reinterpreting Uint64x4 as
|
||||
// Uint32x8, the even positions hold the low 32 bits of each uint64.
|
||||
//
|
||||
//go:nosplit
|
||||
func mredLazy4(xp, yp *[4]uint64, Q, MRedConstant uint64) archsimd.Uint64x4 {
|
||||
xv := archsimd.LoadUint64x4(xp)
|
||||
yv := archsimd.LoadUint64x4(yp)
|
||||
qv := archsimd.BroadcastUint64x4(Q)
|
||||
mredv := archsimd.BroadcastUint64x4(MRedConstant)
|
||||
|
||||
ahi, alo := mul64x4(xv, yv)
|
||||
|
||||
// t = alo * MRedConstant (low 64 bits only)
|
||||
t := mulLow64x4(alo, mredv)
|
||||
|
||||
// H = high 64 bits of (t * Q)
|
||||
H, _ := mul64x4(t, qv)
|
||||
|
||||
// result = ahi - H + Q
|
||||
return ahi.Sub(H).Add(qv)
|
||||
}
|
||||
|
||||
// mul64x4 computes 4 parallel 64x64->128 multiplies, returning high and low 64-bit halves.
|
||||
// Uses VPMULUDQ (AVX2) decomposition into 32-bit widening multiplies.
|
||||
//
|
||||
//go:nosplit
|
||||
func mul64x4(a, b archsimd.Uint64x4) (hi, lo archsimd.Uint64x4) {
|
||||
// Decompose: a = a_hi<<32 | a_lo, b = b_hi<<32 | b_lo
|
||||
a32 := a.AsUint32x8() // even indices = lo32, odd indices = hi32
|
||||
b32 := b.AsUint32x8()
|
||||
|
||||
// VPMULUDQ: multiplies even-indexed 32-bit elements -> 64-bit results
|
||||
// p0[i] = a_lo[i] * b_lo[i] (4 x 64-bit results)
|
||||
p0 := a32.MulEvenWiden(b32)
|
||||
|
||||
// For cross terms, shift to get hi32 into even positions
|
||||
aHi := a.ShiftAllRight(32).AsUint32x8()
|
||||
bHi := b.ShiftAllRight(32).AsUint32x8()
|
||||
|
||||
// p1[i] = a_lo[i] * b_hi[i]
|
||||
p1 := a32.MulEvenWiden(bHi)
|
||||
// p2[i] = a_hi[i] * b_lo[i]
|
||||
p2 := aHi.MulEvenWiden(b32)
|
||||
// p3[i] = a_hi[i] * b_hi[i]
|
||||
p3 := aHi.MulEvenWiden(bHi)
|
||||
|
||||
// Combine: result = p0 + (p1 + p2)<<32 + p3<<64
|
||||
// lo = p0 + (p1 + p2)<<32 (mod 2^64)
|
||||
// hi = p3 + (p1 + p2)>>32 + carries
|
||||
|
||||
mid := p1.Add(p2)
|
||||
|
||||
// Carry from mid addition: if mid < p1, overflow occurred (adds 1<<32 to hi)
|
||||
midCarry := mid.Less(p1)
|
||||
midCarryVal := archsimd.BroadcastUint64x4(1 << 32).Masked(midCarry)
|
||||
|
||||
lo = p0.Add(mid.ShiftAllLeft(32))
|
||||
|
||||
// Carry from lo addition: if lo < p0, overflow occurred (adds 1 to hi)
|
||||
loCarry := lo.Less(p0)
|
||||
loCarryVal := archsimd.BroadcastUint64x4(1).Masked(loCarry)
|
||||
|
||||
hi = p3.Add(mid.ShiftAllRight(32)).Add(midCarryVal).Add(loCarryVal)
|
||||
|
||||
return hi, lo
|
||||
}
|
||||
|
||||
// mulLow64x4 computes 4 parallel 64x64 multiplies, returning only the low 64 bits.
|
||||
// Uses VPMULUDQ (AVX2) decomposition.
|
||||
//
|
||||
//go:nosplit
|
||||
func mulLow64x4(a, b archsimd.Uint64x4) archsimd.Uint64x4 {
|
||||
a32 := a.AsUint32x8()
|
||||
b32 := b.AsUint32x8()
|
||||
|
||||
// p0 = a_lo * b_lo (full 64-bit, but we only need low 64 of final result)
|
||||
p0 := a32.MulEvenWiden(b32)
|
||||
|
||||
// Cross terms: only the low 32 bits of each contribute to the final low 64
|
||||
aHi := a.ShiftAllRight(32).AsUint32x8()
|
||||
bHi := b.ShiftAllRight(32).AsUint32x8()
|
||||
|
||||
p1 := a32.MulEvenWiden(bHi) // a_lo * b_hi
|
||||
p2 := aHi.MulEvenWiden(b32) // a_hi * b_lo
|
||||
|
||||
// lo = p0 + (p1 + p2) << 32 (mod 2^64)
|
||||
return p0.Add(p1.Add(p2).ShiftAllLeft(32))
|
||||
}
|
||||
|
||||
// nttCoreLazyAccel attempts SIMD-accelerated NTT.
|
||||
// Returns true if acceleration was applied, false to fall back to scalar.
|
||||
func nttCoreLazyAccel(p1, p2 []uint64, N int, Q, MRedConstant uint64, roots []uint64) bool {
|
||||
if !simdAVX2 || N < MinimumRingDegreeForLoopUnrolledNTT {
|
||||
return false
|
||||
}
|
||||
|
||||
var j1, j2, t int
|
||||
var F uint64
|
||||
|
||||
fourQ := 4 * Q
|
||||
twoQ := 2 * Q
|
||||
|
||||
// First round: copy p1 -> p2 with butterfly, using SIMD for 4-wide MRedLazy
|
||||
t = N >> 1
|
||||
F = roots[1]
|
||||
psiArr := [4]uint64{F, F, F, F}
|
||||
twoQv := archsimd.BroadcastUint64x4(twoQ)
|
||||
|
||||
for jx, jy := 0, t; jx < t; jx, jy = jx+4, jy+4 {
|
||||
/* #nosec G103 -- aligned 4-element access within bounds */
|
||||
xin := (*[4]uint64)(unsafe.Pointer(&p1[jx]))
|
||||
yin := (*[4]uint64)(unsafe.Pointer(&p1[jy]))
|
||||
|
||||
V := mredLazy4(yin, &psiArr, Q, MRedConstant)
|
||||
xv := archsimd.LoadUint64x4(xin)
|
||||
|
||||
xv.Add(V).StoreSlice(p2[jx:])
|
||||
xv.Add(twoQv).Sub(V).StoreSlice(p2[jy:])
|
||||
}
|
||||
|
||||
// Remaining stages
|
||||
var reduce bool
|
||||
for m := 2; m < N; m <<= 1 {
|
||||
/* #nosec G115 -- m cannot be negative */
|
||||
reduce = (bits.Len64(uint64(m))&1 == 1)
|
||||
t >>= 1
|
||||
|
||||
if t >= 8 {
|
||||
for i := 0; i < m; i++ {
|
||||
j1 = (i * t) << 1
|
||||
j2 = j1 + t
|
||||
F = roots[m+i]
|
||||
psiArr = [4]uint64{F, F, F, F}
|
||||
fourQv := archsimd.BroadcastUint64x4(fourQ)
|
||||
|
||||
if reduce {
|
||||
for jx, jy := j1, j1+t; jx < j2; jx, jy = jx+4, jy+4 {
|
||||
xp := (*[4]uint64)(unsafe.Pointer(&p2[jx]))
|
||||
yp := (*[4]uint64)(unsafe.Pointer(&p2[jy]))
|
||||
|
||||
xv := archsimd.LoadUint64x4(xp)
|
||||
// Conditional reduction: if U >= fourQ, U -= fourQ
|
||||
needReduce := xv.GreaterEqual(fourQv)
|
||||
xv = xv.Sub(fourQv.Masked(needReduce))
|
||||
|
||||
V := mredLazy4(yp, &psiArr, Q, MRedConstant)
|
||||
xout := xv.Add(V)
|
||||
yout := xv.Add(twoQv).Sub(V)
|
||||
|
||||
xout.Store(xp)
|
||||
yout.Store(yp)
|
||||
}
|
||||
} else {
|
||||
for jx, jy := j1, j1+t; jx < j2; jx, jy = jx+4, jy+4 {
|
||||
xp := (*[4]uint64)(unsafe.Pointer(&p2[jx]))
|
||||
yp := (*[4]uint64)(unsafe.Pointer(&p2[jy]))
|
||||
|
||||
xv := archsimd.LoadUint64x4(xp)
|
||||
V := mredLazy4(yp, &psiArr, Q, MRedConstant)
|
||||
xout := xv.Add(V)
|
||||
yout := xv.Add(twoQv).Sub(V)
|
||||
|
||||
xout.Store(xp)
|
||||
yout.Store(yp)
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Small t stages: fall back to scalar (different memory access patterns)
|
||||
nttSmallTStagesScalar(p2, m, t, Q, MRedConstant, twoQ, fourQ, roots, reduce)
|
||||
}
|
||||
}
|
||||
|
||||
return true
|
||||
}
|
||||
|
||||
// inttCoreLazyAccel attempts SIMD-accelerated INTT.
|
||||
// Returns true if acceleration was applied, false to fall back to scalar.
|
||||
func inttCoreLazyAccel(p1, p2 []uint64, N int, Q, MRedConstant uint64, roots []uint64) bool {
|
||||
if !simdAVX2 || N < MinimumRingDegreeForLoopUnrolledNTT {
|
||||
return false
|
||||
}
|
||||
|
||||
var h, t int
|
||||
var F uint64
|
||||
|
||||
t = 1
|
||||
h = N >> 1
|
||||
twoQ := Q << 1
|
||||
fourQ := Q << 2
|
||||
|
||||
// First round (t=1): scalar, different access pattern per butterfly
|
||||
for i, j := h, 0; i < 2*h; i, j = i+8, j+16 {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[8]uint64)(unsafe.Pointer(&roots[i]))
|
||||
xin := (*[16]uint64)(unsafe.Pointer(&p1[j]))
|
||||
xout := (*[16]uint64)(unsafe.Pointer(&p2[j]))
|
||||
|
||||
xout[0], xout[1] = invbutterfly(xin[0], xin[1], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[2], xout[3] = invbutterfly(xin[2], xin[3], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[4], xout[5] = invbutterfly(xin[4], xin[5], psi[2], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[6], xout[7] = invbutterfly(xin[6], xin[7], psi[3], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[8], xout[9] = invbutterfly(xin[8], xin[9], psi[4], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[10], xout[11] = invbutterfly(xin[10], xin[11], psi[5], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[12], xout[13] = invbutterfly(xin[12], xin[13], psi[6], twoQ, fourQ, Q, MRedConstant)
|
||||
xout[14], xout[15] = invbutterfly(xin[14], xin[15], psi[7], twoQ, fourQ, Q, MRedConstant)
|
||||
}
|
||||
|
||||
t <<= 1
|
||||
twoQv := archsimd.BroadcastUint64x4(twoQ)
|
||||
fourQv := archsimd.BroadcastUint64x4(fourQ)
|
||||
for m := N >> 1; m > 1; m >>= 1 {
|
||||
h = m >> 1
|
||||
|
||||
if t >= 8 {
|
||||
for i, j1, j2 := 0, 0, t; i < h; i, j1, j2 = i+1, j1+2*t, j2+2*t {
|
||||
F = roots[h+i]
|
||||
psiArr := [4]uint64{F, F, F, F}
|
||||
|
||||
for jx, jy := j1, j1+t; jx < j2; jx, jy = jx+4, jy+4 {
|
||||
xp := (*[4]uint64)(unsafe.Pointer(&p2[jx]))
|
||||
yp := (*[4]uint64)(unsafe.Pointer(&p2[jy]))
|
||||
|
||||
xv := archsimd.LoadUint64x4(xp)
|
||||
yv := archsimd.LoadUint64x4(yp)
|
||||
|
||||
// X = U + V; if X >= twoQ, X -= twoQ
|
||||
sum := xv.Add(yv)
|
||||
needReduce := sum.GreaterEqual(twoQv)
|
||||
xout := sum.Sub(twoQv.Masked(needReduce))
|
||||
|
||||
// Y = MRedLazy(U + fourQ - V, Psi, Q, MRedConstant)
|
||||
diff := xv.Add(fourQv).Sub(yv)
|
||||
var diffArr [4]uint64
|
||||
diff.Store(&diffArr)
|
||||
yout := mredLazy4(&diffArr, &psiArr, Q, MRedConstant)
|
||||
|
||||
xout.Store(xp)
|
||||
yout.Store(yp)
|
||||
}
|
||||
}
|
||||
} else if t == 4 {
|
||||
for i, j1 := h, 0; i < 2*h; i, j1 = i+2, j1+4*t {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[2]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
x[0], x[4] = invbutterfly(x[0], x[4], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[1], x[5] = invbutterfly(x[1], x[5], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[2], x[6] = invbutterfly(x[2], x[6], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[3], x[7] = invbutterfly(x[3], x[7], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[8], x[12] = invbutterfly(x[8], x[12], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[9], x[13] = invbutterfly(x[9], x[13], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[10], x[14] = invbutterfly(x[10], x[14], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[11], x[15] = invbutterfly(x[11], x[15], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
}
|
||||
} else {
|
||||
for i, j1 := h, 0; i < 2*h; i, j1 = i+4, j1+8*t {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[4]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
x[0], x[2] = invbutterfly(x[0], x[2], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[1], x[3] = invbutterfly(x[1], x[3], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[4], x[6] = invbutterfly(x[4], x[6], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[5], x[7] = invbutterfly(x[5], x[7], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[8], x[10] = invbutterfly(x[8], x[10], psi[2], twoQ, fourQ, Q, MRedConstant)
|
||||
x[9], x[11] = invbutterfly(x[9], x[11], psi[2], twoQ, fourQ, Q, MRedConstant)
|
||||
x[12], x[14] = invbutterfly(x[12], x[14], psi[3], twoQ, fourQ, Q, MRedConstant)
|
||||
x[13], x[15] = invbutterfly(x[13], x[15], psi[3], twoQ, fourQ, Q, MRedConstant)
|
||||
}
|
||||
}
|
||||
|
||||
t <<= 1
|
||||
}
|
||||
|
||||
return true
|
||||
}
|
||||
|
||||
// nttSmallTStagesScalar handles NTT butterfly stages where t < 8.
|
||||
// These stages use non-contiguous memory access patterns (stride != 1)
|
||||
// that don't benefit from 4-wide contiguous SIMD loads/stores.
|
||||
func nttSmallTStagesScalar(p2 []uint64, m, t int, Q, MRedConstant, twoQ, fourQ uint64, roots []uint64, reduce bool) {
|
||||
var V uint64
|
||||
|
||||
if t == 4 {
|
||||
if reduce {
|
||||
for i, j1 := m, 0; i < 2*m; i, j1 = i+2, j1+4*t {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[2]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
x[0], x[4] = butterfly(x[0], x[4], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[1], x[5] = butterfly(x[1], x[5], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[2], x[6] = butterfly(x[2], x[6], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[3], x[7] = butterfly(x[3], x[7], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[8], x[12] = butterfly(x[8], x[12], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[9], x[13] = butterfly(x[9], x[13], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[10], x[14] = butterfly(x[10], x[14], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[11], x[15] = butterfly(x[11], x[15], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
}
|
||||
} else {
|
||||
for i, j1 := m, 0; i < 2*m; i, j1 = i+2, j1+4*t {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[2]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
V = MRedLazy(x[4], psi[0], Q, MRedConstant)
|
||||
x[0], x[4] = x[0]+V, x[0]+twoQ-V
|
||||
V = MRedLazy(x[5], psi[0], Q, MRedConstant)
|
||||
x[1], x[5] = x[1]+V, x[1]+twoQ-V
|
||||
V = MRedLazy(x[6], psi[0], Q, MRedConstant)
|
||||
x[2], x[6] = x[2]+V, x[2]+twoQ-V
|
||||
V = MRedLazy(x[7], psi[0], Q, MRedConstant)
|
||||
x[3], x[7] = x[3]+V, x[3]+twoQ-V
|
||||
V = MRedLazy(x[12], psi[1], Q, MRedConstant)
|
||||
x[8], x[12] = x[8]+V, x[8]+twoQ-V
|
||||
V = MRedLazy(x[13], psi[1], Q, MRedConstant)
|
||||
x[9], x[13] = x[9]+V, x[9]+twoQ-V
|
||||
V = MRedLazy(x[14], psi[1], Q, MRedConstant)
|
||||
x[10], x[14] = x[10]+V, x[10]+twoQ-V
|
||||
V = MRedLazy(x[15], psi[1], Q, MRedConstant)
|
||||
x[11], x[15] = x[11]+V, x[11]+twoQ-V
|
||||
}
|
||||
}
|
||||
} else if t == 2 {
|
||||
if reduce {
|
||||
for i, j1 := m, 0; i < 2*m; i, j1 = i+4, j1+8*t {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[4]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
x[0], x[2] = butterfly(x[0], x[2], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[1], x[3] = butterfly(x[1], x[3], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[4], x[6] = butterfly(x[4], x[6], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[5], x[7] = butterfly(x[5], x[7], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[8], x[10] = butterfly(x[8], x[10], psi[2], twoQ, fourQ, Q, MRedConstant)
|
||||
x[9], x[11] = butterfly(x[9], x[11], psi[2], twoQ, fourQ, Q, MRedConstant)
|
||||
x[12], x[14] = butterfly(x[12], x[14], psi[3], twoQ, fourQ, Q, MRedConstant)
|
||||
x[13], x[15] = butterfly(x[13], x[15], psi[3], twoQ, fourQ, Q, MRedConstant)
|
||||
}
|
||||
} else {
|
||||
for i, j1 := m, 0; i < 2*m; i, j1 = i+4, j1+8*t {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[4]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
V = MRedLazy(x[2], psi[0], Q, MRedConstant)
|
||||
x[0], x[2] = x[0]+V, x[0]+twoQ-V
|
||||
V = MRedLazy(x[3], psi[0], Q, MRedConstant)
|
||||
x[1], x[3] = x[1]+V, x[1]+twoQ-V
|
||||
V = MRedLazy(x[6], psi[1], Q, MRedConstant)
|
||||
x[4], x[6] = x[4]+V, x[4]+twoQ-V
|
||||
V = MRedLazy(x[7], psi[1], Q, MRedConstant)
|
||||
x[5], x[7] = x[5]+V, x[5]+twoQ-V
|
||||
V = MRedLazy(x[10], psi[2], Q, MRedConstant)
|
||||
x[8], x[10] = x[8]+V, x[8]+twoQ-V
|
||||
V = MRedLazy(x[11], psi[2], Q, MRedConstant)
|
||||
x[9], x[11] = x[9]+V, x[9]+twoQ-V
|
||||
V = MRedLazy(x[14], psi[3], Q, MRedConstant)
|
||||
x[12], x[14] = x[12]+V, x[12]+twoQ-V
|
||||
V = MRedLazy(x[15], psi[3], Q, MRedConstant)
|
||||
x[13], x[15] = x[13]+V, x[13]+twoQ-V
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// t == 1
|
||||
for i, j1 := m, 0; i < 2*m; i, j1 = i+8, j1+16 {
|
||||
/* #nosec G103 -- behavior and consequences well understood */
|
||||
psi := (*[8]uint64)(unsafe.Pointer(&roots[i]))
|
||||
x := (*[16]uint64)(unsafe.Pointer(&p2[j1]))
|
||||
x[0], x[1] = butterfly(x[0], x[1], psi[0], twoQ, fourQ, Q, MRedConstant)
|
||||
x[2], x[3] = butterfly(x[2], x[3], psi[1], twoQ, fourQ, Q, MRedConstant)
|
||||
x[4], x[5] = butterfly(x[4], x[5], psi[2], twoQ, fourQ, Q, MRedConstant)
|
||||
x[6], x[7] = butterfly(x[6], x[7], psi[3], twoQ, fourQ, Q, MRedConstant)
|
||||
x[8], x[9] = butterfly(x[8], x[9], psi[4], twoQ, fourQ, Q, MRedConstant)
|
||||
x[10], x[11] = butterfly(x[10], x[11], psi[5], twoQ, fourQ, Q, MRedConstant)
|
||||
x[12], x[13] = butterfly(x[12], x[13], psi[6], twoQ, fourQ, Q, MRedConstant)
|
||||
x[14], x[15] = butterfly(x[14], x[15], psi[7], twoQ, fourQ, Q, MRedConstant)
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,18 @@
|
||||
//go:build !(goexperiment.simd && amd64)
|
||||
|
||||
// SIMD stub for non-AVX2 / non-experimental-simd builds.
|
||||
// LP-107 Phase 3 — verbatim from luxfi/lattice/v7/ring/ntt_simd_stub.go.
|
||||
|
||||
package canonical
|
||||
|
||||
// nttCoreLazyAccel is the stub for non-SIMD builds.
|
||||
// Returns false to indicate no acceleration was applied.
|
||||
func nttCoreLazyAccel(p1, p2 []uint64, N int, Q, MRedConstant uint64, roots []uint64) bool {
|
||||
return false
|
||||
}
|
||||
|
||||
// inttCoreLazyAccel is the stub for non-SIMD builds.
|
||||
// Returns false to indicate no acceleration was applied.
|
||||
func inttCoreLazyAccel(p1, p2 []uint64, N int, Q, MRedConstant uint64, roots []uint64) bool {
|
||||
return false
|
||||
}
|
||||
@@ -0,0 +1,269 @@
|
||||
// Copyright (c) 2026 Lux Industries Inc.
|
||||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
//
|
||||
// Body migrated from luxfi/lattice/v7/ring/subring.go.
|
||||
// LP-107 Phase 3 — luxfi/math is the canonical owner of the Lattigo-derived
|
||||
// Montgomery NTT body and SubRing layout. luxfi/lattice/ring re-exports.
|
||||
//
|
||||
// Changes from upstream:
|
||||
// - factorization.GetFactors -> GetFactors (same package)
|
||||
// - utils.BitReverse64 -> BitReverse64 (same package)
|
||||
// - bignum.NewInt(uint64) -> new(big.Int).SetUint64(...)
|
||||
|
||||
package canonical
|
||||
|
||||
import (
|
||||
"fmt"
|
||||
"math/big"
|
||||
"math/bits"
|
||||
)
|
||||
|
||||
// MinimumRingDegreeForLoopUnrolledOperations is the minimum ring degree required to
|
||||
// safely perform loop-unrolled operations.
|
||||
const MinimumRingDegreeForLoopUnrolledOperations = 8
|
||||
|
||||
// Type is the type of ring used by the cryptographic scheme.
|
||||
type Type int
|
||||
|
||||
// Standard and ConjugateInvariant are two types of Rings.
|
||||
const (
|
||||
Standard = Type(0) // Z[X]/(X^N + 1) (Default)
|
||||
ConjugateInvariant = Type(1) // Z[X+X^-1]/(X^2N + 1)
|
||||
)
|
||||
|
||||
// String returns the string representation of the ring Type.
|
||||
func (rt Type) String() string {
|
||||
switch rt {
|
||||
case Standard:
|
||||
return "Standard"
|
||||
case ConjugateInvariant:
|
||||
return "ConjugateInvariant"
|
||||
default:
|
||||
return "Invalid"
|
||||
}
|
||||
}
|
||||
|
||||
// SubRing is a struct storing precomputation
|
||||
// for fast modular reduction and NTT for
|
||||
// a given modulus.
|
||||
type SubRing struct {
|
||||
ntt NumberTheoreticTransformer
|
||||
|
||||
// Polynomial nb.Coefficients
|
||||
N int
|
||||
|
||||
// Modulus
|
||||
Modulus uint64
|
||||
|
||||
// Unique factors of Modulus-1
|
||||
Factors []uint64
|
||||
|
||||
// 2^bit_length(Modulus) - 1
|
||||
Mask uint64
|
||||
|
||||
// Fast reduction constants
|
||||
BRedConstant [2]uint64 // Barrett Reduction
|
||||
MRedConstant uint64 // Montgomery Reduction
|
||||
|
||||
*NTTTable // NTT related constants
|
||||
}
|
||||
|
||||
// NewSubRing creates a new SubRing with the standard NTT.
|
||||
// NTT constants still need to be generated using .GenNTTConstants(NthRoot uint64).
|
||||
func NewSubRing(N int, Modulus uint64) (s *SubRing, err error) {
|
||||
return NewSubRingWithCustomNTT(N, Modulus, NewNumberTheoreticTransformerStandard, 2*N)
|
||||
}
|
||||
|
||||
// NewSubRingWithCustomNTT creates a new SubRing with degree N and modulus Modulus with user-defined NTT transform and primitive Nth root of unity.
|
||||
// Modulus should be equal to 1 modulo the root of unity.
|
||||
// N must be a power of two larger than 8. An error is returned with a nil *SubRing in the case of non NTT-enabling parameters.
|
||||
func NewSubRingWithCustomNTT(N int, Modulus uint64, ntt func(*SubRing, int) NumberTheoreticTransformer, NthRoot int) (s *SubRing, err error) {
|
||||
|
||||
// Checks if N is a power of 2
|
||||
if N < MinimumRingDegreeForLoopUnrolledOperations || (N&(N-1)) != 0 && N != 0 {
|
||||
return nil, fmt.Errorf("invalid ring degree: must be a power of 2 greater than %d", MinimumRingDegreeForLoopUnrolledOperations)
|
||||
}
|
||||
|
||||
if NthRoot <= 0 {
|
||||
panic(fmt.Errorf("invalid NthRoot: NthRoot=%d should be greater than 0", NthRoot))
|
||||
}
|
||||
|
||||
s = &SubRing{}
|
||||
|
||||
s.N = N
|
||||
|
||||
s.Modulus = Modulus
|
||||
/* #nosec G115 -- Modulus is ensured to be greater than 0 */
|
||||
s.Mask = (1 << uint64(bits.Len64(Modulus-1))) - 1
|
||||
|
||||
// Computes the fast modular reduction constants for the Ring
|
||||
s.BRedConstant = GenBRedConstant(Modulus)
|
||||
|
||||
// If qi is not a power of 2, we can compute the MRed (otherwise, it
|
||||
// would return an error as there is no valid Montgomery form mod a power of 2)
|
||||
if (Modulus&(Modulus-1)) != 0 && Modulus != 0 {
|
||||
s.MRedConstant = GenMRedConstant(Modulus)
|
||||
}
|
||||
|
||||
s.NTTTable = new(NTTTable)
|
||||
s.NthRoot = uint64(NthRoot)
|
||||
|
||||
s.ntt = ntt(s, N)
|
||||
|
||||
return
|
||||
}
|
||||
|
||||
// Type returns the Type of subring which might be either `Standard` or `ConjugateInvariant`.
|
||||
func (s *SubRing) Type() Type {
|
||||
switch s.ntt.(type) {
|
||||
case NumberTheoreticTransformerStandard:
|
||||
return Standard
|
||||
case NumberTheoreticTransformerConjugateInvariant:
|
||||
return ConjugateInvariant
|
||||
default:
|
||||
// Sanity check
|
||||
panic(fmt.Errorf("invalid NumberTheoreticTransformer type: %T", s.ntt))
|
||||
}
|
||||
}
|
||||
|
||||
// GenerateNTTConstants generates the NTT constant for the target SubRing.
|
||||
// The fields `PrimitiveRoot` and `Factors` can be set manually to
|
||||
// bypass the search for the primitive root (which requires to
|
||||
// factor Modulus-1) and speedup the generation of the constants.
|
||||
func (s *SubRing) GenerateNTTConstants() (err error) {
|
||||
|
||||
if s.N == 0 || s.Modulus == 0 {
|
||||
return fmt.Errorf("invalid t parameters (missing)")
|
||||
}
|
||||
|
||||
Modulus := s.Modulus
|
||||
NthRoot := s.NthRoot
|
||||
|
||||
// Checks if each qi is prime and equal to 1 mod NthRoot
|
||||
if !IsPrime(Modulus) {
|
||||
return fmt.Errorf("invalid modulus: %d is not prime)", Modulus)
|
||||
}
|
||||
|
||||
if Modulus&(NthRoot-1) != 1 {
|
||||
return fmt.Errorf("invalid modulus: %d != 1 mod NthRoot)", Modulus)
|
||||
}
|
||||
|
||||
// It is possible to manually set the primitive root along with the factors of q-1.
|
||||
// This is notably useful when marshalling the SubRing, to avoid re-factoring q-1.
|
||||
// If both are set, then checks that that the root is indeed primitive.
|
||||
// Else, factorize q-1 and finds a primitive root.
|
||||
if s.PrimitiveRoot != 0 && s.Factors != nil {
|
||||
if err = CheckPrimitiveRoot(s.PrimitiveRoot, s.Modulus, s.Factors); err != nil {
|
||||
return
|
||||
}
|
||||
} else {
|
||||
if s.PrimitiveRoot, s.Factors, err = PrimitiveRoot(Modulus, s.Factors); err != nil {
|
||||
return
|
||||
}
|
||||
}
|
||||
|
||||
logNthRoot := int(bits.Len64(NthRoot>>1) - 1)
|
||||
|
||||
// 1.1 Computes N^(-1) mod Q in Montgomery form
|
||||
s.NInv = MForm(ModExp(NthRoot>>1, Modulus-2, Modulus), Modulus, s.BRedConstant)
|
||||
|
||||
// 1.2 Computes Psi and PsiInv in Montgomery form
|
||||
|
||||
// Computes Psi and PsiInv in Montgomery form
|
||||
PsiMont := MForm(ModExp(s.PrimitiveRoot, (Modulus-1)/NthRoot, Modulus), Modulus, s.BRedConstant)
|
||||
PsiInvMont := MForm(ModExp(s.PrimitiveRoot, Modulus-((Modulus-1)/NthRoot)-1, Modulus), Modulus, s.BRedConstant)
|
||||
|
||||
s.RootsForward = make([]uint64, NthRoot>>1)
|
||||
s.RootsBackward = make([]uint64, NthRoot>>1)
|
||||
|
||||
s.RootsForward[0] = MForm(1, Modulus, s.BRedConstant)
|
||||
s.RootsBackward[0] = MForm(1, Modulus, s.BRedConstant)
|
||||
|
||||
// Computes nttPsi[j] = nttPsi[j-1]*Psi and RootsBackward[j] = RootsBackward[j-1]*PsiInv
|
||||
for j := uint64(1); j < NthRoot>>1; j++ {
|
||||
|
||||
indexReversePrev := BitReverse64(j-1, logNthRoot)
|
||||
indexReverseNext := BitReverse64(j, logNthRoot)
|
||||
|
||||
s.RootsForward[indexReverseNext] = MRed(s.RootsForward[indexReversePrev], PsiMont, Modulus, s.MRedConstant)
|
||||
s.RootsBackward[indexReverseNext] = MRed(s.RootsBackward[indexReversePrev], PsiInvMont, Modulus, s.MRedConstant)
|
||||
}
|
||||
|
||||
return
|
||||
}
|
||||
|
||||
// PrimitiveRoot computes the smallest primitive root of the given prime q
|
||||
// The unique factors of q-1 can be given to speed up the search for the root.
|
||||
func PrimitiveRoot(q uint64, factors []uint64) (uint64, []uint64, error) {
|
||||
|
||||
if factors != nil {
|
||||
if err := CheckFactors(q-1, factors); err != nil {
|
||||
return 0, factors, err
|
||||
}
|
||||
} else {
|
||||
|
||||
factorsBig := GetFactors(new(big.Int).SetUint64(q - 1)) //Factor q-1, might be slow
|
||||
|
||||
factors = make([]uint64, len(factorsBig))
|
||||
for i := range factors {
|
||||
factors[i] = factorsBig[i].Uint64()
|
||||
}
|
||||
}
|
||||
|
||||
notFoundPrimitiveRoot := true
|
||||
|
||||
var g uint64 = 2
|
||||
|
||||
for notFoundPrimitiveRoot {
|
||||
g++
|
||||
for _, factor := range factors {
|
||||
// if for any factor of q-1, g^(q-1)/factor = 1 mod q, g is not a primitive root
|
||||
if ModExp(g, (q-1)/factor, q) == 1 {
|
||||
notFoundPrimitiveRoot = true
|
||||
break
|
||||
}
|
||||
notFoundPrimitiveRoot = false
|
||||
}
|
||||
}
|
||||
|
||||
return g, factors, nil
|
||||
}
|
||||
|
||||
// CheckFactors checks that the given list of factors contains
|
||||
// all the unique primes of m.
|
||||
func CheckFactors(m uint64, factors []uint64) (err error) {
|
||||
|
||||
for _, factor := range factors {
|
||||
|
||||
if !IsPrime(factor) {
|
||||
return fmt.Errorf("composite factor")
|
||||
}
|
||||
|
||||
for m%factor == 0 {
|
||||
m /= factor
|
||||
}
|
||||
}
|
||||
|
||||
if m != 1 {
|
||||
return fmt.Errorf("incomplete factor list")
|
||||
}
|
||||
|
||||
return
|
||||
}
|
||||
|
||||
// CheckPrimitiveRoot checks that g is a valid primitive root mod q,
|
||||
// given the factors of q-1.
|
||||
func CheckPrimitiveRoot(g, q uint64, factors []uint64) (err error) {
|
||||
|
||||
if err = CheckFactors(q-1, factors); err != nil {
|
||||
return
|
||||
}
|
||||
|
||||
for _, factor := range factors {
|
||||
if ModExp(g, (q-1)/factor, q) == 1 {
|
||||
return fmt.Errorf("invalid primitive root")
|
||||
}
|
||||
}
|
||||
|
||||
return
|
||||
}
|
||||
@@ -0,0 +1,147 @@
|
||||
// Copyright (c) 2026 Lux Industries Inc.
|
||||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
package canonical
|
||||
|
||||
import (
|
||||
"math/rand/v2"
|
||||
"testing"
|
||||
)
|
||||
|
||||
// pulsarN256 is the Pulsar/LP-073 NTT instance: N=256, Q=2^48 + 4*N + 1.
|
||||
const (
|
||||
pulsarN = 256
|
||||
pulsarQ = uint64(0x1000000004A01)
|
||||
)
|
||||
|
||||
// TestSubRing_RoundTrip exercises the full canonical path:
|
||||
// constructor + GenerateNTTConstants + forward NTT + inverse NTT.
|
||||
func TestSubRing_RoundTrip(t *testing.T) {
|
||||
sr, err := NewSubRing(pulsarN, pulsarQ)
|
||||
if err != nil {
|
||||
t.Fatalf("NewSubRing: %v", err)
|
||||
}
|
||||
if err := sr.GenerateNTTConstants(); err != nil {
|
||||
t.Fatalf("GenerateNTTConstants: %v", err)
|
||||
}
|
||||
|
||||
r := rand.New(rand.NewPCG(0xdeadbeef, 0x12345678))
|
||||
a := make([]uint64, pulsarN)
|
||||
for i := range a {
|
||||
a[i] = r.Uint64() % pulsarQ
|
||||
}
|
||||
saved := make([]uint64, pulsarN)
|
||||
copy(saved, a)
|
||||
|
||||
sr.NTT(a, a)
|
||||
sr.INTT(a, a)
|
||||
|
||||
for i := range a {
|
||||
if a[i] != saved[i] {
|
||||
t.Fatalf("round-trip [%d]: %d != %d", i, a[i], saved[i])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// TestSubRing_LazyVariants exercises NTTLazy / INTTLazy.
|
||||
func TestSubRing_LazyVariants(t *testing.T) {
|
||||
sr, err := NewSubRing(pulsarN, pulsarQ)
|
||||
if err != nil {
|
||||
t.Fatalf("NewSubRing: %v", err)
|
||||
}
|
||||
if err := sr.GenerateNTTConstants(); err != nil {
|
||||
t.Fatalf("GenerateNTTConstants: %v", err)
|
||||
}
|
||||
|
||||
a := make([]uint64, pulsarN)
|
||||
r := rand.New(rand.NewPCG(0xa5a5, 0))
|
||||
for i := range a {
|
||||
a[i] = r.Uint64() % pulsarQ
|
||||
}
|
||||
b := make([]uint64, pulsarN)
|
||||
copy(b, a)
|
||||
|
||||
sr.NTTLazy(a, a)
|
||||
for i := range a {
|
||||
// Lazy bound: a[i] in [0, 6q-2]
|
||||
if a[i] >= 6*pulsarQ {
|
||||
t.Fatalf("NTTLazy [%d]=%d exceeds 6q-2", i, a[i])
|
||||
}
|
||||
}
|
||||
sr.INTTLazy(a, a)
|
||||
// INTTLazy(NTTLazy(x)) is not identity (different reduction
|
||||
// regimes); run NTT/INTT on b for the round-trip check.
|
||||
sr.NTT(b, b)
|
||||
sr.INTT(b, b)
|
||||
rng := rand.New(rand.NewPCG(0xa5a5, 0))
|
||||
for i := range b {
|
||||
want := rng.Uint64() % pulsarQ
|
||||
if b[i] != want {
|
||||
t.Fatalf("NTT/INTT round-trip [%d]: %d != %d", i, b[i], want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// TestModularReduction_Primitives smoke-tests the scalar primitives so
|
||||
// a regression in any of MRed / MRedLazy / MForm / IMForm / BRed /
|
||||
// CRed surfaces here, not just behind the NTT body.
|
||||
func TestModularReduction_Primitives(t *testing.T) {
|
||||
q := pulsarQ
|
||||
brc := GenBRedConstant(q)
|
||||
mrc := GenMRedConstant(q)
|
||||
|
||||
x, y := uint64(0xC0FFEE), uint64(0xBEEF)
|
||||
|
||||
// MForm <-> IMForm round-trip
|
||||
xm := MForm(x, q, brc)
|
||||
if got := IMForm(xm, q, mrc); got != x {
|
||||
t.Errorf("IMForm(MForm(x))=%d, want %d", got, x)
|
||||
}
|
||||
|
||||
// MRed should match (x*y) mod q after demontgomery-ization
|
||||
ym := MForm(y, q, brc)
|
||||
xmy := MRed(xm, ym, q, mrc)
|
||||
xy := IMForm(xmy, q, mrc)
|
||||
want := MulModSlow(x, y, q)
|
||||
if xy != want {
|
||||
t.Errorf("MRed(MForm(x), MForm(y))=%d (de-mont), want %d", xy, want)
|
||||
}
|
||||
|
||||
// CRed bound
|
||||
if got := CRed(2*q-1, q); got >= q {
|
||||
t.Errorf("CRed(2q-1)=%d, want < q", got)
|
||||
}
|
||||
}
|
||||
|
||||
// MulModSlow is reference (x*y) mod q via 128-bit multiply, used only
|
||||
// in tests.
|
||||
func MulModSlow(x, y, q uint64) uint64 {
|
||||
hi, lo := mul64(x, y)
|
||||
return mod128(hi, lo, q)
|
||||
}
|
||||
|
||||
func mul64(x, y uint64) (hi, lo uint64) {
|
||||
xLo, xHi := x&0xffffffff, x>>32
|
||||
yLo, yHi := y&0xffffffff, y>>32
|
||||
t := xLo*yLo
|
||||
lo = t & 0xffffffff
|
||||
t = (t >> 32) + xHi*yLo
|
||||
w1 := t & 0xffffffff
|
||||
w2 := t >> 32
|
||||
t = xLo*yHi + w1
|
||||
lo |= (t & 0xffffffff) << 32
|
||||
hi = xHi*yHi + w2 + (t >> 32)
|
||||
return
|
||||
}
|
||||
|
||||
func mod128(hi, lo, q uint64) uint64 {
|
||||
// Slow but obviously-correct: divmod via long division.
|
||||
for i := 0; i < 64; i++ {
|
||||
hi = (hi << 1) | (lo >> 63)
|
||||
lo <<= 1
|
||||
if hi >= q {
|
||||
hi -= q
|
||||
}
|
||||
}
|
||||
return hi
|
||||
}
|
||||
@@ -0,0 +1,63 @@
|
||||
// Copyright (c) 2026 Lux Industries Inc.
|
||||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
//
|
||||
// Helpers migrated from luxfi/lattice/v7/ring (utils.go, primes.go) and
|
||||
// luxfi/lattice/v7/utils (BitReverse64). LP-107 Phase 3.
|
||||
|
||||
package canonical
|
||||
|
||||
import (
|
||||
"math/big"
|
||||
"math/bits"
|
||||
)
|
||||
|
||||
// BitReverse64 returns the bit-reverse value of the input value, within a context of 2^bitLen.
|
||||
func BitReverse64(index uint64, bitLen int) uint64 {
|
||||
return bits.Reverse64(index) >> (64 - bitLen)
|
||||
}
|
||||
|
||||
// IsPrime applies the Baillie-PSW, which is 100% accurate for numbers below 2^64.
|
||||
func IsPrime(x uint64) bool {
|
||||
return new(big.Int).SetUint64(x).ProbablyPrime(0)
|
||||
}
|
||||
|
||||
// ModExp performs the modular exponentiation x^e mod p,
|
||||
// x and p are required to be at most 64 bits to avoid an overflow.
|
||||
func ModExp(x, e, p uint64) (result uint64) {
|
||||
brc := GenBRedConstant(p)
|
||||
result = 1
|
||||
for i := e; i > 0; i >>= 1 {
|
||||
if i&1 == 1 {
|
||||
result = BRed(result, x, p, brc)
|
||||
}
|
||||
x = BRed(x, x, p, brc)
|
||||
}
|
||||
return result
|
||||
}
|
||||
|
||||
// ModExpPow2 performs the modular exponentiation x^e mod p, where p is a power of two,
|
||||
// x and p are required to be at most 64 bits to avoid an overflow.
|
||||
func ModExpPow2(x, e, p uint64) (result uint64) {
|
||||
result = 1
|
||||
for i := e; i > 0; i >>= 1 {
|
||||
if i&1 == 1 {
|
||||
result *= x
|
||||
}
|
||||
x *= x
|
||||
}
|
||||
return result & (p - 1)
|
||||
}
|
||||
|
||||
// ModexpMontgomery performs the modular exponentiation x^e mod p,
|
||||
// where x is in Montgomery form, and returns x^e in Montgomery form.
|
||||
func ModexpMontgomery(x uint64, e int, q, mredconstant uint64, bredconstant [2]uint64) (result uint64) {
|
||||
result = MForm(1, q, bredconstant)
|
||||
|
||||
for i := e; i > 0; i >>= 1 {
|
||||
if i&1 == 1 {
|
||||
result = MRed(result, x, q, mredconstant)
|
||||
}
|
||||
x = MRed(x, x, q, mredconstant)
|
||||
}
|
||||
return result
|
||||
}
|
||||
@@ -0,0 +1,70 @@
|
||||
// Copyright (c) 2026 Lux Industries Inc.
|
||||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
//
|
||||
// Body migrated from luxfi/lattice/v7/ring/vec_ops.go.
|
||||
// LP-107 Phase 3 — only the helpers required by the NTT body
|
||||
// (reducevec, mulscalarmontgomeryvec, mulscalarmontgomerylazyvec).
|
||||
|
||||
package canonical
|
||||
|
||||
import "unsafe"
|
||||
|
||||
// reducevec applies BRedAdd to every coefficient of p1 -> p2 (size multiple of 8).
|
||||
func reducevec(p1, p2 []uint64, modulus uint64, bredconstant [2]uint64) {
|
||||
N := len(p1)
|
||||
for j := 0; j < N; j = j + 8 {
|
||||
/* #nosec G103 -- behavior and consequences well understood, possible buffer overflow if len(p1)%8 */
|
||||
x := (*[8]uint64)(unsafe.Pointer(&p1[j]))
|
||||
/* #nosec G103 -- behavior and consequences well understood, possible buffer overflow if len(p2)%8 */
|
||||
z := (*[8]uint64)(unsafe.Pointer(&p2[j]))
|
||||
|
||||
z[0] = BRedAdd(x[0], modulus, bredconstant)
|
||||
z[1] = BRedAdd(x[1], modulus, bredconstant)
|
||||
z[2] = BRedAdd(x[2], modulus, bredconstant)
|
||||
z[3] = BRedAdd(x[3], modulus, bredconstant)
|
||||
z[4] = BRedAdd(x[4], modulus, bredconstant)
|
||||
z[5] = BRedAdd(x[5], modulus, bredconstant)
|
||||
z[6] = BRedAdd(x[6], modulus, bredconstant)
|
||||
z[7] = BRedAdd(x[7], modulus, bredconstant)
|
||||
}
|
||||
}
|
||||
|
||||
// mulscalarmontgomeryvec multiplies p1 by scalarMont (Montgomery form) into p2.
|
||||
func mulscalarmontgomeryvec(p1 []uint64, scalarMont uint64, p2 []uint64, modulus, mredconstant uint64) {
|
||||
N := len(p1)
|
||||
for j := 0; j < N; j = j + 8 {
|
||||
/* #nosec G103 -- behavior and consequences well understood, possible buffer overflow if len(p1)%8 */
|
||||
x := (*[8]uint64)(unsafe.Pointer(&p1[j]))
|
||||
/* #nosec G103 -- behavior and consequences well understood, possible buffer overflow if len(p2)%8 */
|
||||
z := (*[8]uint64)(unsafe.Pointer(&p2[j]))
|
||||
|
||||
z[0] = MRed(x[0], scalarMont, modulus, mredconstant)
|
||||
z[1] = MRed(x[1], scalarMont, modulus, mredconstant)
|
||||
z[2] = MRed(x[2], scalarMont, modulus, mredconstant)
|
||||
z[3] = MRed(x[3], scalarMont, modulus, mredconstant)
|
||||
z[4] = MRed(x[4], scalarMont, modulus, mredconstant)
|
||||
z[5] = MRed(x[5], scalarMont, modulus, mredconstant)
|
||||
z[6] = MRed(x[6], scalarMont, modulus, mredconstant)
|
||||
z[7] = MRed(x[7], scalarMont, modulus, mredconstant)
|
||||
}
|
||||
}
|
||||
|
||||
// mulscalarmontgomerylazyvec is the lazy variant of mulscalarmontgomeryvec.
|
||||
func mulscalarmontgomerylazyvec(p1 []uint64, scalarMont uint64, p2 []uint64, modulus, mredconstant uint64) {
|
||||
N := len(p1)
|
||||
for j := 0; j < N; j = j + 8 {
|
||||
/* #nosec G103 -- behavior and consequences well understood, possible buffer overflow if len(p1)%8 */
|
||||
x := (*[8]uint64)(unsafe.Pointer(&p1[j]))
|
||||
/* #nosec G103 -- behavior and consequences well understood, possible buffer overflow if len(p2)%8 */
|
||||
z := (*[8]uint64)(unsafe.Pointer(&p2[j]))
|
||||
|
||||
z[0] = MRedLazy(x[0], scalarMont, modulus, mredconstant)
|
||||
z[1] = MRedLazy(x[1], scalarMont, modulus, mredconstant)
|
||||
z[2] = MRedLazy(x[2], scalarMont, modulus, mredconstant)
|
||||
z[3] = MRedLazy(x[3], scalarMont, modulus, mredconstant)
|
||||
z[4] = MRedLazy(x[4], scalarMont, modulus, mredconstant)
|
||||
z[5] = MRedLazy(x[5], scalarMont, modulus, mredconstant)
|
||||
z[6] = MRedLazy(x[6], scalarMont, modulus, mredconstant)
|
||||
z[7] = MRedLazy(x[7], scalarMont, modulus, mredconstant)
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,149 @@
|
||||
// Body migrated from luxfi/lattice/v7/utils/factorization/weierstrass.go.
|
||||
// LP-107 Phase 3 — randInt inlined to drop the lattice/v7/utils/sampling
|
||||
// import so this package has zero lattice dependencies.
|
||||
|
||||
package canonical
|
||||
|
||||
import (
|
||||
"crypto/rand"
|
||||
"math/big"
|
||||
)
|
||||
|
||||
// randInt generates a random Int in [0, max-1].
|
||||
// Inlined from luxfi/lattice/v7/utils/randInt.
|
||||
func randInt(max *big.Int) (n *big.Int) {
|
||||
var err error
|
||||
if n, err = rand.Int(rand.Reader, max); err != nil {
|
||||
panic(err)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// Weierstrass is an elliptic curve y^2 = x^3 + ax + b mod N.
|
||||
type Weierstrass struct {
|
||||
A, B, N *big.Int
|
||||
}
|
||||
|
||||
// Point represents an elliptic curve point in standard coordinates.
|
||||
type Point struct {
|
||||
X, Y *big.Int
|
||||
}
|
||||
|
||||
// Add adds two Weierstrass points together with respect
|
||||
// to the underlying Weierstrass curve.
|
||||
// This method does not check if the points lie on
|
||||
// the underlying curve.
|
||||
func (w *Weierstrass) Add(P, Q Point) Point {
|
||||
|
||||
tmp := new(big.Int)
|
||||
|
||||
xR, yR := new(big.Int), new(big.Int)
|
||||
|
||||
if P.X.Cmp(tmp.SetUint64(0)) == 0 && P.Y.Cmp(tmp.SetUint64(1)) == 0 {
|
||||
|
||||
return Point{xR.Set(Q.X), yR.Set(Q.Y)}
|
||||
}
|
||||
|
||||
if Q.X.Cmp(tmp.SetUint64(0)) == 0 && Q.Y.Cmp(tmp.SetUint64(1)) == 0 {
|
||||
return Point{xR.Set(P.X), yR.Set(P.Y)}
|
||||
}
|
||||
|
||||
xP, yP := P.X, P.Y
|
||||
xQ, yQ := Q.X, Q.Y
|
||||
|
||||
N := w.N
|
||||
|
||||
if xP.Cmp(xQ) == 0 && yP.Cmp(new(big.Int).Sub(N, yQ)) == 0 {
|
||||
return Point{xR.SetUint64(0), yR.SetUint64(0)}
|
||||
}
|
||||
|
||||
S := new(big.Int) // slope
|
||||
|
||||
if xP != xQ {
|
||||
|
||||
// S = (yQ-yP)/(xQ-xP)
|
||||
S.Sub(yQ, yP)
|
||||
tmp.Sub(xQ, xP)
|
||||
tmp.ModInverse(tmp, N)
|
||||
S.Mul(S, tmp)
|
||||
S.Mod(S, N)
|
||||
|
||||
} else {
|
||||
|
||||
// S = (3*(xP^2) + a)/(2*yP)
|
||||
S.Mul(xP, xP)
|
||||
S.Mod(S, N)
|
||||
S.Mul(S, new(big.Int).SetUint64(3))
|
||||
S.Add(S, w.A)
|
||||
S.Mod(S, N)
|
||||
tmp.Add(yP, yP)
|
||||
tmp.ModInverse(tmp, N)
|
||||
S.Mul(S, tmp)
|
||||
S.Mod(S, N)
|
||||
}
|
||||
|
||||
// s^2 - xP - xQ
|
||||
xR.Mul(S, S)
|
||||
xR.Mod(xR, N)
|
||||
xR.Sub(xR, xP)
|
||||
xR.Sub(xR, xQ)
|
||||
xR.Mod(xR, N)
|
||||
|
||||
// s*(xP-xR)-yP
|
||||
yR.Sub(xP, xR)
|
||||
yR.Mul(yR, S)
|
||||
yR.Mod(yR, N)
|
||||
yR.Sub(yR, yP)
|
||||
yR.Mod(yR, N)
|
||||
|
||||
return Point{X: xR, Y: yR}
|
||||
}
|
||||
|
||||
// NewRandomWeierstrassCurve generates a new random Weierstrass curve modulo N,
|
||||
// along with a random point that lies on the curve.
|
||||
func NewRandomWeierstrassCurve(N *big.Int) (Weierstrass, Point) {
|
||||
|
||||
var A, B, xG, yG *big.Int
|
||||
for {
|
||||
|
||||
// Select random values for A, xG and yG
|
||||
A = randInt(N)
|
||||
xG = randInt(N)
|
||||
yG = randInt(N)
|
||||
|
||||
// Deduces B from Y^2 = X^3 + A * X + B evaluated at point (xG, yG)
|
||||
yGpow2 := new(big.Int).Mul(yG, yG)
|
||||
yGpow2.Mod(yGpow2, N)
|
||||
|
||||
xGpow3 := new(big.Int).Mul(xG, xG)
|
||||
xGpow3.Mod(xGpow3, N)
|
||||
xGpow3.Sub(xGpow3, A)
|
||||
xGpow3.Mul(xGpow3, xG)
|
||||
xGpow3.Mod(xGpow3, N)
|
||||
|
||||
B = new(big.Int).Sub(yGpow2, xGpow3) // B = yG^2 - xG*(xG^2 - A)
|
||||
B.Mod(B, N)
|
||||
|
||||
// Checks that 4A^3 + 27B^2 != 0
|
||||
fourACube := new(big.Int).Add(A, A)
|
||||
fourACube.Mul(fourACube, fourACube)
|
||||
fourACube.Mod(fourACube, N)
|
||||
fourACube.Mul(fourACube, A)
|
||||
|
||||
twentySevenBSquare := new(big.Int).Mul(B, B)
|
||||
twentySevenBSquare.Mod(twentySevenBSquare, N)
|
||||
twentySevenBSquare.Mul(twentySevenBSquare, new(big.Int).SetUint64(27))
|
||||
twentySevenBSquare.Mod(twentySevenBSquare, N)
|
||||
|
||||
jInvariantQuotient := new(big.Int).Add(fourACube, twentySevenBSquare)
|
||||
jInvariantQuotient.Mod(jInvariantQuotient, N)
|
||||
|
||||
if jInvariantQuotient.Cmp(new(big.Int).SetUint64(0)) != 0 && new(big.Int).GCD(nil, nil, N, jInvariantQuotient).Cmp(new(big.Int).SetUint64(1)) == 0 {
|
||||
return Weierstrass{
|
||||
A: A,
|
||||
B: B,
|
||||
N: N,
|
||||
}, Point{X: xG, Y: yG}
|
||||
}
|
||||
}
|
||||
}
|
||||
+28
-26
@@ -7,22 +7,22 @@ import (
|
||||
"fmt"
|
||||
"sync"
|
||||
|
||||
"github.com/luxfi/math/ntt/canonical"
|
||||
"github.com/luxfi/math/params"
|
||||
|
||||
"github.com/luxfi/lattice/v7/ring"
|
||||
)
|
||||
|
||||
// pureGoBackend is the canonical pure-Go NTT realization. It delegates
|
||||
// to github.com/luxfi/lattice/v7/ring's SubRing.NTT / INTT — the
|
||||
// canonical Lattigo-derived Montgomery NTT — so callers see no
|
||||
// behavior change vs the v0.1.x lattice path.
|
||||
// pureGoBackend is the canonical pure-Go NTT realization. It owns the
|
||||
// Lattigo-derived Montgomery NTT body via the internal canonical
|
||||
// package; lattice/v7/ring re-exports the same body for downstream
|
||||
// consumers, so callers see no behavior change vs the v0.1.x lattice
|
||||
// path.
|
||||
//
|
||||
// LP-107 Phase 3 will invert this dependency: the canonical kernel
|
||||
// body will live in this package, and luxfi/lattice will import
|
||||
// luxfi/math/ntt to expose ring.SubRing.NTT.
|
||||
// LP-107 Phase 3 (this file): the canonical kernel body now lives in
|
||||
// luxfi/math/ntt/internal/canonical. luxfi/lattice/ring is a thin shim
|
||||
// that delegates to it.
|
||||
type pureGoBackend struct {
|
||||
mu sync.RWMutex
|
||||
rings map[params.NTTParamID]*ring.Ring
|
||||
mu sync.RWMutex
|
||||
subRings map[params.NTTParamID]*canonical.SubRing
|
||||
}
|
||||
|
||||
// PureGoBackend returns the singleton pure-Go NTT backend. Always
|
||||
@@ -32,7 +32,7 @@ func PureGoBackend() Backend {
|
||||
}
|
||||
|
||||
var thePureGo = pureGoBackend{
|
||||
rings: make(map[params.NTTParamID]*ring.Ring),
|
||||
subRings: make(map[params.NTTParamID]*canonical.SubRing),
|
||||
}
|
||||
|
||||
func init() {
|
||||
@@ -49,26 +49,30 @@ func (b *pureGoBackend) Supports(p *Params) bool {
|
||||
return p != nil && p.Validate() == nil
|
||||
}
|
||||
|
||||
// resolveRing returns or builds the cached *ring.Ring for p.
|
||||
func (b *pureGoBackend) resolveRing(p *Params) (*ring.Ring, error) {
|
||||
// resolveSubRing returns or builds the cached *canonical.SubRing for p.
|
||||
func (b *pureGoBackend) resolveSubRing(p *Params) (*canonical.SubRing, error) {
|
||||
b.mu.RLock()
|
||||
r, ok := b.rings[p.ID]
|
||||
sr, ok := b.subRings[p.ID]
|
||||
b.mu.RUnlock()
|
||||
if ok {
|
||||
return r, nil
|
||||
return sr, nil
|
||||
}
|
||||
b.mu.Lock()
|
||||
defer b.mu.Unlock()
|
||||
if r, ok := b.rings[p.ID]; ok {
|
||||
return r, nil
|
||||
if sr, ok := b.subRings[p.ID]; ok {
|
||||
return sr, nil
|
||||
}
|
||||
rr, err := ring.NewRing(int(p.N), []uint64{p.Q})
|
||||
sr, err := canonical.NewSubRing(int(p.N), p.Q)
|
||||
if err != nil {
|
||||
return nil, fmt.Errorf("ntt(pure-go): ring.NewRing(N=%d, Q=%d): %w",
|
||||
return nil, fmt.Errorf("ntt(pure-go): canonical.NewSubRing(N=%d, Q=%d): %w",
|
||||
p.N, p.Q, err)
|
||||
}
|
||||
b.rings[p.ID] = rr
|
||||
return rr, nil
|
||||
if err := sr.GenerateNTTConstants(); err != nil {
|
||||
return nil, fmt.Errorf("ntt(pure-go): GenerateNTTConstants(N=%d, Q=%d): %w",
|
||||
p.N, p.Q, err)
|
||||
}
|
||||
b.subRings[p.ID] = sr
|
||||
return sr, nil
|
||||
}
|
||||
|
||||
// Forward implements Backend.
|
||||
@@ -81,11 +85,10 @@ func (b *pureGoBackend) Forward(dst []uint64, p *Params, batch uint32) error {
|
||||
return fmt.Errorf("ntt(pure-go): buffer too small: need %d got %d",
|
||||
int(batch)*N, len(dst))
|
||||
}
|
||||
r, err := b.resolveRing(p)
|
||||
sr, err := b.resolveSubRing(p)
|
||||
if err != nil {
|
||||
return err
|
||||
}
|
||||
sr := r.SubRings[0]
|
||||
for i := uint32(0); i < batch; i++ {
|
||||
off := int(i) * N
|
||||
sr.NTT(dst[off:off+N], dst[off:off+N])
|
||||
@@ -103,11 +106,10 @@ func (b *pureGoBackend) Inverse(dst []uint64, p *Params, batch uint32) error {
|
||||
return fmt.Errorf("ntt(pure-go): buffer too small: need %d got %d",
|
||||
int(batch)*N, len(dst))
|
||||
}
|
||||
r, err := b.resolveRing(p)
|
||||
sr, err := b.resolveSubRing(p)
|
||||
if err != nil {
|
||||
return err
|
||||
}
|
||||
sr := r.SubRings[0]
|
||||
for i := uint32(0); i < batch; i++ {
|
||||
off := int(i) * N
|
||||
sr.INTT(dst[off:off+N], dst[off:off+N])
|
||||
|
||||
Reference in New Issue
Block a user