Tutorial: NTT Forward and Inverse Transform
This tutorial walks you through building a batched Number Theoretic Transform (NTT) application using cuPQC-NTT. The NTT is the modular-arithmetic analogue of the FFT: it turns polynomial multiplication into cheap pointwise multiplication, which is why it sits on the critical path of lattice-based post-quantum schemes, fully homomorphic encryption, and zero-knowledge proof systems.
You'll build a forward and inverse NTT round-trip at N = 1024 over the KoalaBear prime field, processing one polynomial per CUDA block, and verify that INTT(NTT(f)) == f for every coefficient of every polynomial in the batch.
cuPQC-NTT requires cuPQC SDK 0.6.0 or newer, the release that introduced it. cuPQC 0.6 targets compute capabilities 8.0, 8.6, 8.7, 8.9, and 9.0. The library implements the cyclic transform, so the arithmetic is in the ring modulo x^N - 1.
Step 1: Project Setup
Clone the cuPQC repository:
This will download all examples including the NTT example. The Makefile in this directory will compile all examples, including example_ntt.
Step 2: Include Required Headers
Start by including the necessary cuPQC SDK headers:
The ntt.hpp header provides the NTT descriptors and operators, the built-in prime fields and their precomputed roots of unity, and the shared-memory sizing helpers.
Step 3: Define the Forward and Inverse NTT Types
An NTT type is composed from operators, in the same style as the other cuPQC libraries. You need two descriptors, one per direction:
constexpr uint32_t NTT_N = 1024;
using ForwardNTT = decltype(Algorithm<algorithm::NTT>()
+ Direction<nttDirection::FORWARD>()
+ Precision<uint32_t>()
+ Size<NTT_N>()
+ Block()
+ BlockDim<128>());
using InverseNTT = decltype(Algorithm<algorithm::NTT>()
+ Direction<nttDirection::INVERSE>()
+ Precision<uint32_t>()
+ Size<NTT_N>()
+ Block()
+ BlockDim<128>());
Algorithm<algorithm::NTT>(): selects the Number Theoretic Transform algorithm
Direction<nttDirection::FORWARD>() / Direction<nttDirection::INVERSE>(): selects the transform direction; the two descriptors are otherwise identical
Precision<uint32_t>(): selects the coefficient and modulus width
Size<NTT_N>(): sets the transform length, here 1024 = 2^10
Block(): requests block execution, meaning one whole transform is performed cooperatively by a CUDA block
BlockDim<128>(): sets the thread count that cooperates on the transform
The composition is resolved at compile time, and the resulting type exposes the configuration back to your code as ForwardNTT::Size and ForwardNTT::BlockDim, which you'll use for the launch geometry.
Step 4: Generate the Twiddle Tables
Twiddle factors are generated once on the device and then reused across the entire batch. Each direction needs its own table:
template<class NTT>
__global__ void make_twiddles_kernel(uint32_t* twiddles, const uint32_t p, const uint32_t g) {
NTT().make_twiddles(twiddles, p, g);
}
template<class NTT>
__global__ void transform_twiddles_to_mont_kernel(uint32_t* twiddles, const uint32_t p) {
NTT().transform_twiddles_to_mont(twiddles, p);
}
make_twiddles(twiddles, p, g): fills the table for the prime p and generator g. The forward table uses a primitive N-th root of unity g; the inverse table uses g^-1 mod p.
transform_twiddles_to_mont(twiddles, p): converts the table into the Montgomery domain the device transforms operate in
make_twiddles runs single-threaded (<<<1, 1>>>) and only needs to be called once per prime and generator pair, while the Montgomery conversion is parallel over the table.
The inverse generator is derived on the host by Fermat's little theorem, g^-1 mod p = g^(p-2) mod p:
// Host-side modular exponentiation: base^exp mod m.
static uint32_t host_modpow(uint64_t base, uint64_t exp, uint64_t mod) {
uint64_t result = 1;
base %= mod;
while (exp > 0) {
if (exp & 1) { result = result * base % mod; }
base = base * base % mod;
exp >>= 1;
}
return static_cast<uint32_t>(result);
}
Step 5: Create the Forward NTT Kernel
Device transforms operate on Montgomery-domain data, so conversion happens at the memory boundary rather than per operation. A kernel loads into shared memory, executes, and stores back:
__global__ void forward_ntt_kernel(uint32_t* polys, const uint32_t* twiddles, const nttConst<uint32_t> ntt_const) {
uint32_t* poly = polys + blockIdx.x * ForwardNTT::Size;
extern __shared__ uint32_t sdata[];
ForwardNTT().load_to_mont(sdata, poly, ntt_const);
__syncthreads();
ForwardNTT().execute(sdata, twiddles, ntt_const.p);
__syncthreads();
ForwardNTT().store_from_mont(sdata, poly, ntt_const);
}
load_to_mont(sdata, poly, ntt_const): reads the polynomial from global memory into the shared-memory workspace, converting into the Montgomery domain on the way in
execute(sdata, twiddles, ntt_const.p): performs the transform in place on the shared-memory workspace
store_from_mont(sdata, poly, ntt_const): writes the workspace back to global memory, converting out of the Montgomery domain
nttConst<uint32_t> ntt_const: carries the prime and the derived Montgomery constants; construct it once on the host with nttConst<uint32_t>(p) and pass it by value
The __syncthreads() calls between stages are required because the whole block cooperates on one shared buffer: every thread must finish loading before any thread starts transforming, and the transform must complete before any thread stores.
blockIdx.x * ForwardNTT::Size is what makes the kernel batched. Each block picks a different polynomial out of the contiguous device array and works on it independently.
Step 6: Create the Inverse NTT Kernel
The inverse kernel follows the same load-execute-store pattern, with one extra argument:
__global__ void inverse_ntt_kernel(uint32_t* polys, const uint32_t* inv_twiddles,
const nttConst<uint32_t> ntt_const, const uint32_t N_inv) {
uint32_t* poly = polys + blockIdx.x * InverseNTT::Size;
extern __shared__ uint32_t sdata[];
InverseNTT().load_to_mont(sdata, poly, ntt_const);
__syncthreads();
InverseNTT().execute(sdata, inv_twiddles, ntt_const.p, N_inv);
__syncthreads();
InverseNTT().store_from_mont(sdata, poly, ntt_const);
}
The inverse execute additionally takes N_inv, which is N^-1 mod p. The transform applies it to normalise its output, so that the round-trip returns the original coefficients rather than a scaled copy of them. Obtain it with the n_inv<N>(p) helper.
Step 7: Implement the Host Round-Trip
The host function sizes the shared-memory workspace, allocates device memory, generates the twiddle tables, launches the two kernels back to back, and verifies the result:
bool ntt_round_trip(const unsigned int batch,
const uint32_t p,
const uint32_t g,
const uint32_t g_inv,
const uint32_t N_inv)
{
constexpr uint32_t N = ForwardNTT::Size;
constexpr size_t smem = ntt_shared_workspace_size<N, uint32_t>();
// Initialise batch polynomials with values 0, 1, 2, ... (mod p).
std::vector<uint32_t> h_data(N * batch);
// ...
const std::vector<uint32_t> reference = h_data;
// Allocate device memory for polynomials and twiddle tables, and upload the input.
// ...
// Generate twiddle tables once - shared across all polynomials in the batch.
make_twiddles_kernel<ForwardNTT><<<1, 1>>>(d_twiddles, p, g);
make_twiddles_kernel<InverseNTT><<<1, 1>>>(d_inv_twiddles, p, g_inv);
transform_twiddles_to_mont_kernel<ForwardNTT><<<1, ForwardNTT::BlockDim>>>(d_twiddles, p);
transform_twiddles_to_mont_kernel<InverseNTT><<<1, InverseNTT::BlockDim>>>(d_inv_twiddles, p);
const nttConst<uint32_t> ntt_const(p);
// Forward NTT: transform each polynomial into its NTT representation.
forward_ntt_kernel<<<batch, ForwardNTT::BlockDim, smem>>>(d_data, d_twiddles, ntt_const);
// Inverse NTT: recover the original polynomials.
inverse_ntt_kernel<<<batch, InverseNTT::BlockDim, smem>>>(d_data, d_inv_twiddles, ntt_const, N_inv);
// Transfer results back to host and free device memory.
// ...
// Verify: INTT(NTT(f)) == f for every coefficient of every polynomial.
unsigned int errors = 0;
for (size_t i = 0; i < reference.size(); i++) {
if (h_data[i] != reference[i]) {
errors++;
}
}
if (errors == 0) {
std::cout << "Round-trip OK: " << batch << " polynomial(s) of degree " << N << " verified." << std::endl;
return true;
}
std::cout << "Round-trip FAILED: " << errors << " coefficient(s) did not match." << std::endl;
return false;
}
ntt_shared_workspace_size<N, uint32_t>() returns the number of bytes one block needs for its workspace. It is a constexpr, so you can pass it straight through as the dynamic shared-memory argument of the launch configuration.
The launch geometry is <<<batch, ForwardNTT::BlockDim, smem>>>: one block per polynomial, BlockDim threads per block cooperating on that polynomial, and smem bytes of shared memory per block. Both kernels operate in place on d_data, so the inverse pass reads exactly what the forward pass wrote.
Step 8: Create the Main Application
The main function picks the field, derives the two scalars the inverse transform needs, and runs the round-trip over a batch of 8 polynomials:
int main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) {
std::cout << "================================================================\n";
std::cout << "Number Theoretic Transform (NTT) Example\n";
std::cout << "================================================================\n\n";
// ...
// KoalaBear prime: p = 2^31 - 2^24 + 1 = 2130706433.
// KoalaBear_primitive_root_10 is a primitive 2^10-th root of unity mod p.
constexpr uint32_t p = cupqc::KoalaBear;
constexpr uint32_t g = cupqc::KoalaBear_primitive_root_10;
// g^{-1} mod p = g^{p-2} mod p (Fermat's little theorem)
const uint32_t g_inv = host_modpow(g, p - 2, p);
// N^{-1} mod p, applied by the inverse NTT to normalise the output
const uint32_t N_inv = n_inv<NTT_N>(p);
const unsigned int batch = 8;
const bool ok = ntt_round_trip(batch, p, g, g_inv, N_inv);
std::cout << "\n" << (ok ? "Example completed successfully.\n"
: "Example failed.\n");
return ok ? EXIT_SUCCESS : EXIT_FAILURE;
}
BabyBear and KoalaBear are built in, along with precomputed roots named *_primitive_root_S for transform sizes 2^S where 10 <= S <= 24. Because this example uses Size<1024> and 1024 = 2^10, the matching root is KoalaBear_primitive_root_10. Custom primes up to 62 bits are also supported.
Step 9: Build and Run
The Makefile will build all examples in the folder. Run the NTT example:
The Makefile expects the cuPQC SDK at /usr/local/cupqc-sdk, or at a path of your choosing by setting the CUPQC_SDK_DIR environment variable or editing the Makefile.
cuPQC-NTT is a static library linked with LTO (-dlto -lcupqc-ntt), built for a specific set of transform sizes and sub-sizes. A configuration the type system accepts but the shipped library does not provide compiles cleanly and then fails to link, so a link error on an unusual Size or SubSize is the expected signal that the combination is not available.
Step 10: Understanding the Output
Expected output:
================================================================
Number Theoretic Transform (NTT) Example
================================================================
This example demonstrates a batched NTT round-trip using cuPQC SDK.
The NTT is the modular-arithmetic analogue of the FFT, and is the
standard way to turn polynomial multiplication into cheap pointwise
multiplication for lattice cryptography, FHE, and ZK proof systems.
Configuration: N=1024, KoalaBear field, block execution
Layout: one polynomial per CUDA block, shared-memory workspace
Round-trip OK: 8 polynomial(s) of degree 1024 verified.
Example completed successfully.
The verification is exact, not approximate: all 1024 * 8 coefficients must match the reference bit for bit. If any coefficient does not match, the program prints Round-trip FAILED with the number of mismatches followed by Example failed., and exits with a non-zero status.
Customization Tips
Change the Precision
Precision selects the coefficient and modulus width. Besides uint32_t, the library supports uint16_t and uint64_t:
using ForwardNTT64 = decltype(Algorithm<algorithm::NTT>()
+ Direction<nttDirection::FORWARD>()
+ Precision<uint64_t>()
+ Size<NTT_N>()
+ Block()
+ BlockDim<128>());
The precision has to be consistent across the descriptor, the twiddle table, the nttConst, and the workspace sizing helper, so switch the element type everywhere at once: nttConst<uint64_t>, ntt_shared_workspace_size<N, uint64_t>(), and the polynomial and twiddle buffers.
Use a Custom Prime
The example uses the built-in cupqc::KoalaBear field and its precomputed root, but custom primes up to 62 bits are also supported. In that case you supply the prime and a primitive N-th root of unity yourself, pass them to make_twiddles, and build the nttConst from your prime. The inverse generator and the normalisation factor come from the same helpers as before: host_modpow(g, p - 2, p) for g^-1 mod p and n_inv<N>(p) for N^-1 mod p.
Tune Size and BlockDim
Size<N> sets the transform length and BlockDim<...> sets how many threads cooperate on it. Both are compile-time choices, and both feed the launch: the shared-memory request follows from Size through ntt_shared_workspace_size, and the thread count comes straight from ForwardNTT::BlockDim. Reading the geometry back off the type rather than hard-coding it, as this example does, means a change to either operator needs no edits at the launch sites.
Move to a Staged Transform for Large N
For large N, a single block's shared memory cannot hold the full workspace. Adding SubSize<M> splits the transform into two kernel passes, which bounds each pass's per-block working set well below the full N. The staged NTT tutorial works through a round-trip at N = 16384 and the rules for choosing M.
Learn More
- Example Source Code - Complete NTT round-trip example implementation
- cuPQC-NTT User Guide - Usage guide with examples
- cuPQC-NTT API Reference - Complete API documentation
- cuPQC-NTT Features - Supported transform sizes and prime fields