Tutorial: BigInt Hex Digits of Pi
This tutorial walks you through a GPU implementation of Bellard's formula for pi, built with cuPQC-BigInt. Instead of accumulating floating-point terms, every subterm of the series is evaluated in fixed-width big-integer arithmetic, so the program produces 256 exact fractional hex digits of pi and validates them against the known leading digits. Along the way you'll use modular exponentiation, division with remainder, shifts, and a warp-level reduction of 320-bit values.
This is the most advanced of the three BigInt examples: the series terms are spread across the whole GPU, each thread does real big-integer work, and the partial sums are reduced first inside each warp and then across the grid.
cuPQC-BigInt 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.
Unlike the other cuPQC libraries, cuPQC-BigInt is exposed entirely as __device__ functions. There are no host-side launch helpers: you call every operation from inside your own kernel, which is exactly what this example does.
Step 1: Project Setup
Clone the cuPQC repository:
This downloads all examples including the hex-digits-of-pi example. The Makefile in this directory will compile all examples, including example_bigint_pihex.
Step 2: Include Headers and Match the Architecture
Start with the cuPQC SDK headers:
#include <algorithm>
#include <array>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <bigint.hpp>
#include <cucheck.hpp>
using namespace cupqc;
using namespace std;
bigint.hpp provides the big-integer descriptor operators and device operations.
cucheck.hpp provides the CUDA_CHECK macro the host code wraps every CUDA call in.
A big-integer descriptor carries an SM<>() operator that selects the architecture code is generated for. The Makefile compiles with -arch=native, so the example derives SM<> from __CUDA_ARCH__ to keep the descriptor in step with whatever is actually being compiled:
#if defined(__CUDA_ARCH__) && (__CUDA_ARCH__ == 900)
#define PIHEX_SM() + SM<900>()
#elif defined(__CUDA_ARCH__) && (__CUDA_ARCH__ == 890)
#define PIHEX_SM() + SM<890>()
// ... 870, 860, 800 ...
#else
#define PIHEX_SM()
#endif
SM<>() is optional and defaults to SM<800> when omitted, which is what the host compilation pass and any architecture without a specialization use. The macro stops at 9.0 because that is the highest SM<> specialization the SDK provides, so a newer architecture falls through to the default rather than naming a specialization that does not exist.
Step 3: Define the Big-Integer Types
The example needs three widths, all owned by a single thread:
using bigint64_cfg = decltype(BitWidth<64>() + Thread() PIHEX_SM());
using bigint320_cfg = decltype(BitWidth<320>() + Thread() PIHEX_SM());
using bigint384_cfg = decltype(BitWidth<384>() + Thread() PIHEX_SM());
using bigint64_t = typename bigint64_cfg::bigint;
using bigint320_t = typename bigint320_cfg::bigint;
using bigint384_t = typename bigint384_cfg::bigint;
BitWidth<BW>() sets the width, which must be a multiple of 32; the limb count is BW / 32.
Thread() gives one thread ownership of the entire value, so no TPI<>() is set.
cfg::bigint is the value type produced by the descriptor.
The three widths each have a job. Denominators in Bellard's formula are always represented as 64-bit values, accumulators are 320 bits wide, and 384 bits gives room to shift a 64-bit residue left by 320 before dividing.
The limb bookkeeping derives the guard and result sizes from the descriptors rather than hard-coding them:
constexpr int BIGINT_WORDS = bigint320_cfg::num_limbs;
// Reserve 64 guard bits because denominators are limited to 64 bits.
constexpr int GUARD_WORDS = bigint64_cfg::num_limbs;
constexpr int RESULT_WORDS = BIGINT_WORDS - GUARD_WORDS;
static_assert(bigint384_cfg::num_limbs == BIGINT_WORDS + GUARD_WORDS);
constexpr int THREADS_PER_BLOCK = 64;
constexpr int REDUCE_THREADS = 1024;
constexpr int BELLARD_TERM_PAIRS = 7;
The calculation carries 320 fractional bits, of which the low 64 are guard bits that absorb rounding from the per-term divisions. They are dropped at the end, leaving RESULT_WORDS limbs — 256 bits, or 64 hex digits — of trustworthy output. Values are little-endian arrays of 32-bit limbs, so limb 0 is least significant.
Step 4: The Math, as Integer Arithmetic
Bellard's formula writes pi as seven alternating-sign subterms. The example's header comment carries the 2^-6 and (-1)^n 2^(-10n) factors into each subterm and splits n into its even and odd cases (n = 2m and n = 2m+1), which turns the seven alternating subterms into seven fixed-sign pairs — fourteen subterms in total, each a positive or negative 2^p / denom.
That restructuring is what makes the GPU version simple: there is no sign bookkeeping left, just seven pairs, each contributing one addition and one subtraction per iteration.
For a single positive subterm 2^p / denom, shifting pi by start_bit and retaining B = 320 fractional bits asks for
When pow >= 0, writing 2^pow = q*denom + r makes the q*2^B part vanish modulo 2^B, leaving floor(r * 2^B / denom) with r = 2^pow mod denom. So one modular exponentiation, one shift, and one division suffice. When -B <= pow < 0, the contribution is directly floor(2^(B+pow) / denom), and when pow < -B it lies below the retained precision and is zero.
Step 5: Evaluate One Subterm
pihex_term is the heart of the example — it turns the derivation above into three big-integer operations:
// term = floor(2^pow * 2^320 / denom) mod 2^320.
__device__ void pihex_term(bigint320_t& term, int64_t pow, uint64_t denom)
{
bigint64_t divisor(denom);
bigint384_t quotient;
bigint64_t remainder;
if (pow > 0) {
bigint64_t residue = bigint64_t(2u).pow_mod(static_cast<uint64_t>(pow), divisor);
bigint384_t dividend(residue.to_uint64());
dividend <<= 320;
dividend.div_rem(divisor, quotient, remainder);
term = narrow_to_320(quotient);
} else if (pow >= -320) {
// Negative powers occur in the tail (and at small start_bit values).
bigint384_t dividend(1u);
dividend <<= static_cast<uint32_t>(320 + pow);
dividend.div_rem(divisor, quotient, remainder);
term = narrow_to_320(quotient);
} else {
term = bigint320_t{};
}
}
bigint64_t(2u).pow_mod(pow, divisor) computes 2^pow mod denom entirely in 64-bit big-integer arithmetic.
residue.to_uint64() extracts that residue so it can seed a wider value.
dividend <<= 320 shifts the residue up by the retained precision, which is why the dividend is 384 bits wide.
dividend.div_rem(divisor, quotient, remainder) produces both the quotient and the remainder. div_rem is the one operation that accepts a narrower divisor than the dividend, which is what allows a 64-bit divisor against a 384-bit dividend.
bigint320_t{} is the zero contribution used when the subterm falls below the retained precision.
Note that the strict pow > 0 test sends pow == 0 into the second branch, where 2^(320+0) / denom computes the same value, so both branches agree at the boundary.
The quotient comes back at the dividend's width, so it is narrowed explicitly. cuPQC-BigInt never implicitly promotes or truncates between widths, so the copy is written out limb by limb:
// Explicitly narrow to the low 320 bits.
__device__ bigint320_t narrow_to_320(const bigint384_t& value)
{
bigint320_t result;
for (int i = 0; i < BIGINT_WORDS; ++i)
result[i] = value[i];
return result;
}
Step 6: Accumulate the Term Pairs
Each of the seven pairs contributes one positive and one negative subterm, which is a single add and a single subtract on the accumulator:
__device__ void accum_pihex_terms(bigint320_t& accum,
int64_t pow_a, uint64_t denom_a,
int64_t pow_b, uint64_t denom_b)
{
bigint320_t term_a;
bigint320_t term_b;
pihex_term(term_a, pow_a, denom_a);
pihex_term(term_b, pow_b, denom_b);
accum += term_a;
accum -= term_b;
}
The accumulator is unsigned and fixed-width, so += and -= wrap modulo 2^320. That is precisely the arithmetic the derivation asks for, which is why negative subterms need no special handling.
accum_bellard_pair then selects one of the seven pairs by term_id and supplies the powers and denominators for iteration n:
__device__ void accum_bellard_pair(bigint320_t& accum, int term_id,
int64_t start_bit, int64_t n)
{
switch (term_id) {
case 0:
accum_pihex_terms(accum, start_bit-20*n-16+5, 8*n+5,
start_bit-20*n-6+5, 8*n+1);
break;
// ... cases 1 through 5 ...
default:
accum_pihex_terms(accum, start_bit-20*n-6, 20*n+9,
start_bit-20*n-16, 20*n+19);
break;
}
}
Every case has the same shape: two pow expressions that decrease by 20 bits per iteration, and two denominators that grow linearly in n. The start_bit offset simply rides along inside each pow.
Step 7: The Calculation Kernel and Warp Reduction
calc_pihex splits the grid seven ways, one group of blocks per term pair, and walks the series with a grid-stride loop:
__global__ void calc_pihex(uint32_t* partials, int64_t start_bit,
int64_t n_max, int blocks_per_term)
{
const int term_id = blockIdx.x / blocks_per_term;
const int local_block = blockIdx.x % blocks_per_term;
const int64_t tid = static_cast<int64_t>(local_block) * blockDim.x + threadIdx.x;
const int64_t stride = static_cast<int64_t>(blocks_per_term) * blockDim.x;
bigint320_t accum;
for (int64_t n = tid; n < n_max; n += stride) {
accum_bellard_pair(accum, term_id, start_bit, n);
}
// ... warp reduction, then store ...
}
Each thread therefore owns a strided subset of the n values for one term pair, and holds a private 320-bit running sum in registers. Because the whole series is one unsigned sum modulo 2^320, the terms can be visited in any order and grouped arbitrarily.
The per-thread sums are then combined inside each warp with shuffles, one limb at a time:
bigint320_t other;
for (int offset = 16; offset > 0; offset >>= 1) {
for (int limb = 0; limb < BIGINT_WORDS; ++limb) {
other[limb] = __shfl_down_sync(0xffffffffu, accum[limb], offset);
}
accum += other;
}
if ((threadIdx.x & 31) == 0) {
const int partial = blockIdx.x * (THREADS_PER_BLOCK / 32) + threadIdx.x / 32;
accum.store(partials, partial);
}
Limb access with accum[limb] is what makes this work: the value's limbs are ordinary 32-bit registers, so a standard __shfl_down_sync tree moves a peer's whole 320-bit value into other, and one big-integer += adds it with correct carry propagation. After five halvings, lane 0 of each warp holds that warp's total and stores it with accum.store(partials, partial), where the second argument is the instance index into a packed array of 320-bit values.
With THREADS_PER_BLOCK of 64, each block produces two partials.
Step 8: Reduce the Partials to a Single Value
The second kernel is a single block of REDUCE_THREADS threads that collapses every partial into one 320-bit value:
__global__ void reduce_pihex(uint32_t* partials, int num_partials)
{
__shared__ uint32_t shared[REDUCE_THREADS * BIGINT_WORDS];
const int tid = threadIdx.x;
bigint320_t accum;
for (int i = tid; i < num_partials; i += REDUCE_THREADS) {
bigint320_t value(partials, i);
accum += value;
}
accum.store(shared, tid);
__syncthreads();
for (int stride = REDUCE_THREADS / 2; stride > 0; stride >>= 1) {
if (tid < stride) {
bigint320_t a(shared, tid);
bigint320_t b(shared, tid + stride);
a += b;
a.store(shared, tid);
}
__syncthreads();
}
if (tid == 0) {
bigint320_t result(shared, 0);
result.store(partials);
}
}
bigint320_t value(partials, i) constructs a value directly from instance i of a packed limb array — the same indexed form as store, which is what lets global and shared memory be treated as arrays of big integers.
a.store(shared, tid) writes back into shared memory at instance tid, so the classic halving tree needs no manual limb indexing.
result.store(partials) writes the final 320-bit sum to the front of the same buffer.
Step 9: Host Orchestration, Timing, and Validation
The host picks a launch geometry from the requested start_bit and the device's SM count:
const int64_t n_max = (start_bit + 320 + 19) / 20;
const int64_t wanted_blocks = (n_max + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
const int64_t max_blocks = static_cast<int64_t>(prop.multiProcessorCount) * 8;
const int blocks_per_term = static_cast<int>(min(wanted_blocks, max_blocks));
const int num_partials = blocks_per_term * BELLARD_TERM_PAIRS *
(THREADS_PER_BLOCK / 32);
n_max is how many iterations are needed to cover start_bit plus the 320 retained bits, since each iteration advances the series by 20 bits. blocks_per_term is one block per THREADS_PER_BLOCK iterations, capped at eight blocks per SM, and the grid is that many blocks for each of the seven term pairs.
The work is then issued on a stream with CUDA events between each phase:
CUDA_CHECK(cudaEventRecord(begin, stream));
calc_pihex<<<blocks_per_term * BELLARD_TERM_PAIRS, THREADS_PER_BLOCK, 0, stream>>>(
partials, start_bit, n_max, blocks_per_term);
CUDA_CHECK(cudaEventRecord(calculated, stream));
reduce_pihex<<<1, REDUCE_THREADS, 0, stream>>>(partials, num_partials);
CUDA_CHECK(cudaEventRecord(reduced, stream));
CUDA_CHECK(cudaMemcpyAsync(result.data(), partials + GUARD_WORDS,
result.size() * sizeof(uint32_t), cudaMemcpyDeviceToHost, stream));
CUDA_CHECK(cudaEventRecord(copied, stream));
CUDA_CHECK(cudaStreamSynchronize(stream));
Three intervals are timed and reported: the calculation kernel, the reduction kernel, and the copy back to the host. The destination array is registered with cudaHostRegister before the stream work so the copy is asynchronous, and unregistered afterwards.
Note the copy offset. Reading from partials + GUARD_WORDS skips the low GUARD_WORDS limbs of the reduced value, which is how the 64 guard bits are discarded — the host only ever sees the trustworthy RESULT_WORDS limbs.
Validation compares those limbs against a table of known digits, stored in cuPQC limb order:
// Known leading hex digits of pi (after the leading "3."), used to validate
// the start_bit=0 run. Stored in cuPQC limb order (LSW first); printed in
// reverse order to recover the natural human-readable hex string.
constexpr uint32_t PI_AT_ZERO[RESULT_WORDS] = {
0xec4e6c89u, 0x082efa98u, 0x299f31d0u, 0xa4093822u,
0x03707344u, 0x13198a2eu, 0x85a308d3u, 0x243f6a88u,
};
Because limb 0 is least significant, the digits are printed from the top limb down, and the check only applies when start_bit is 0:
for (int i = RESULT_WORDS - 1; i >= 0; --i) {
printf("%08x", result[i]);
}
// ...
if (start_bit == 0 && !equal(result.begin(), result.end(), PI_AT_ZERO)) {
fprintf(stderr, "FAIL: result does not match the leading hex digits of pi.\n");
return EXIT_FAILURE;
}
main parses the optional argument with strtoll and rejects anything that is not a fully-consumed nonnegative integer, printing a usage message instead.
Step 10: Build and Run
The Makefile expects the cuPQC SDK at /usr/local/cupqc-sdk, or at a path you supply through the CUPQC_SDK_DIR environment variable. It builds every example in the folder:
The program takes an optional nonnegative starting bit offset, defaulting to 0:
At offset 0 the program validates its output against the known leading hex digits of pi and reports GPU timings for the calculation, reduction, and copy-back phases.
cuPQC-BigInt is a static library built for a specific set of widths per TPI, and it is linked with link-time optimization: the Makefile compiles with -dlto and links -lcupqc-bigint. It also uses -arch=native, which is why the descriptors derive SM<> from __CUDA_ARCH__.
Step 11: Understanding the Output
Expected output at the default offset:
================================================================
BigInt pihex Example
================================================================
This example computes 256 fractional hex digits of pi starting at
bit offset 0, using Bellard's formula evaluated in 320- and
384-bit integer arithmetic across the whole GPU. It combines
pow_mod, div_rem, shifts and warp-level reduction of big integers.
Configuration: 320-bit accumulators, thread execution, 64 guard bits
pihex Example Program Results
pihex start_bit=0 GPU timings: calc=12.340ms reduce=1.230ms memcpy=0.040ms
pihex start_bit=0 (64 hex digits, most-significant first):
243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89
OK: matches the known leading hex digits of pi.
The timing values above are placeholders — actual numbers vary by GPU, and the calculation phase also grows with start_bit because a larger offset requires more series iterations. The digit string is the fractional part of pi in hex, after the leading 3., and the OK: line appears only when the offset is 0. A mismatch prints a FAIL: line on standard error and exits with a failure status.
Running with a nonzero offset prints the same banner and result lines with the offset substituted, but no validation line, since there is no reference table for arbitrary offsets.
Customization Tips
Compute Digits Further Out
The start_bit argument shifts the window of pi that gets computed, and the program adapts to it automatically: n_max grows with the offset, so more series iterations are needed and the calculation kernel does proportionally more work.
The 64 guard bits exist because denominators are limited to 64 bits. Since the largest denominator is approximately start_bit + 338, that reservation corresponds mathematically to start_bit up to roughly 2^64; the int64_t indexing used throughout imposes a lower practical limit of approximately 2^63.
Why This Example Stays on Thread()
All three widths use Thread(), giving each thread a complete value in its registers. That is the right fit here: the per-term work is serial, and the interesting parallelism is across series iterations rather than inside a single value. The warp-shuffle reduction in calc_pihex also assumes one value per thread, so moving to Warp() would mean restructuring the reduction as well as the descriptors.
See the BigInt vector addition tutorial for the trade-off between the two execution modes, and for why a change to BitWidth<>() or TPI<>() can surface as a link error rather than a compiler diagnostic.
Reuse the Big-Integer Reduction Pattern
The two-stage reduction is independent of Bellard's formula and transfers to any kernel that sums big integers. The warp stage shuffles limbs and adds with +=; the grid stage stores partials as a packed array and collapses them with the indexed constructor and store:
bigint320_t value(partials, i); // load instance i
accum += value;
accum.store(partials, instance); // store instance
Because the indexed constructor and store work on both global and shared memory, the same pattern serves for staging partials between kernels and for the in-block halving tree.
Learn More
- Example Source Code - Complete hex digits of pi example implementation
- cuPQC-BigInt User Guide - Usage guide with examples
- cuPQC-BigInt API Reference - Complete API documentation
- cuPQC-BigInt Features - Supported widths, operations, and reduction paths