MPFR with CUDA acceleration (CGBN backend)

This tree is GNU MPFR 4.2.1 with an optional CUDA backend built on NVIDIA's CGBN library, enabled with --enable-cgbn.

It adds two interfaces for evaluating MPFR arithmetic on a GPU — batches over ordinary arrays, and vectors that stay resident in device memory — plus an off-by-default hook that can divert a single scalar operation. Everything else is stock MPFR 4.2.1.

The results are identical to MPFR's own. Every operation produces the same correctly rounded value, the same ternary value and the same exception flags as the corresponding scalar function, for all five rounding modes and all special values. What changes is throughput, not semantics. tests/tcgbn is the proof: it cross-checks every operation against the scalar functions and nothing else.

Without --enable-cgbn the library is bit-for-bit the one MPFR builds today — not a single new symbol is exported, which is why -version-info is unchanged. This follows the convention already used for --enable-decimal-float.

---

mpfr-cgbn-design is a longer technical article on the same material: the design, the algorithms, the full measurements and the analysis of the 32768-bit boundary.

Contents

• Requirements
• Building
• Configure options
• Verifying the build
• Using the backend
• Performance
  • Full timings, all precisions
  • How to read them
• Why performance varies with bit width
• CGBN tuning parameters
• Would a bigger GPU help?
• Implementation notes
• Tools
• Limitations

---

Requirements

GMP	5.0.0 or later, with GMP_NUMB_BITS a multiple of 32 (checked by configure)
CUDA	a toolkit providing nvcc; developed and tested against 13.2
CGBN	the header tree, i.e. a directory containing cgbn/cgbn.h
GPU	compute capability 3.5 or later (CGBN needs warp shuffles and funnel shifts)
C++	nvcc compiles the device code as C++, so the C++ runtime is linked in

CGBN is header-only; nothing of it needs to be built or installed beyond unpacking the headers somewhere.

---

Building

Build out of tree, which keeps the source directory comparable to a pristine 4.2.1:

%>> mkdir -p build && cd build
%>> ../configure --enable-cgbn --with-cuda=/usr/local/cuda-13.2
%>> gmake -j$(nproc)
%>> gmake check

--with-cgbn defaults to /usr/local/include, so it can be omitted when the headers are there. If nvcc is already on $PATH, --with-cuda can be omitted too and the toolkit prefix is taken from the location of nvcc.

A minimal invocation, with everything in default places:

%>> ../configure --enable-cgbn

Static-only build, which is convenient under gdb and compute-sanitizer:

%>> ../configure --enable-cgbn --disable-shared

Cross-compiling, or building for a GPU that is not the build machine's:

%>> ../configure --enable-cgbn --with-cuda-arch=sm_90

Build cost

The device code is instantiated once per field width, one translation unit each, so make -j compiles them in parallel. On a 16-core host the whole library takes roughly 80 seconds, most of it nvcc; the 65536-bit unit alone is about a minute. Only src/mpfr-cgbn-*.cu are affected — an ordinary edit elsewhere in MPFR rebuilds as fast as it always did.

---

Configure options

Option	Meaning
--enable-cgbn	Build the CUDA/CGBN backend. Default: no.
--with-cuda=DIR	CUDA toolkit prefix, the directory containing bin/nvcc. Default: derived from the nvcc found on $PATH.
--with-cgbn=DIR	Directory containing the cgbn/ header tree. Default: /usr/local/include.
--with-cuda-arch=ARCH	Value passed to nvcc -arch, for example sm_90, sm_120 or all-major. Default: native, which targets the GPU of the build machine and therefore requires one to be present.

These compose with every existing MPFR option. Two worth knowing about:

Option	Interaction
--disable-shared	Fine. The device objects are compiled -fPIC and serve both the shared and the static library.
--enable-assert / =full	Fine, and recommended while modifying the backend: it turns on the MPFR_ASSERTD checks in the host-side packing code.

What configure does when --enable-cgbn is given:

1. locates nvcc, either under --with-cuda or on $PATH;
2. checks that $cgbn_path/cgbn/cgbn.h exists;
3. checks that GMP_NUMB_BITS is a multiple of 32 — the backend moves GMP limbs in and out of 32-bit CGBN fields and has no path for other limb sizes;
4. checks for pthread.h, and if present defines MPFR_CGBN_HAVE_PTHREAD so that device access is serialized on a mutex;
5. sets NVCCFLAGS to -O3 -arch=$cuda_arch -Xcompiler -fPIC and CUDA_LIBS to -lcudart -lstdc++ (plus -lpthread);
6. defines MPFR_WANT_CGBN, which is what guards every line of the backend and the four hook sites in mul.c, sqr.c, div.c, sqrt.c and fma.c.

Linking against the result

A program using the backend includes <mpfr-cgbn.h> and links as usual:

%>> gcc -o prog prog.c -lmpfr -lgmp

The shared libmpfr.so records its own dependency on libcudart and libstdc++, so nothing more is needed. Against a static libmpfr.a, add -lcudart -lstdc++ -lpthread and the CUDA library directory.

---

Verifying the build

%>> gmake check                       # 199 tests, including tcgbn
%>> cd tests && ./tcgbn              # just the backend, ~30 s
%>> GMP_CHECK_RANDOMIZE=1 ./tcgbn    # again with a random seed

tests/tcgbn compares the device against the scalar functions over every rounding mode, every width in the ladder and the boundaries just past each rung, the special values (NaN, ±Inf, ±0), exact ties, one-ulp cases, total cancellation, exponent differences that straddle the working width, in-place aliasing on every operand position, batches large enough to be split into chunks, and the device-vector path. It also checks that the batch functions ran where they were supposed to, and that the fused multiply-add test data still distinguishes a fused evaluation from a multiply followed by an add.

Two more checks worth running after touching the device code:

%>> cd tests && make cgbn-verify && ./cgbn-verify   # CGBN itself against GMP
%>> cd src && make check-exported-symbols check-gmp-symbols

---

Using the backend

examples/mpfr-cgbn-example.c is a complete worked example: it exercises every batch operation and every device-vector operation and checks each result against the corresponding scalar function. Build it with

%>> gcc -o mpfr-cgbn-example examples/mpfr-cgbn-example.c -lmpfr -lgmp

Status

#include <mpfr-cgbn.h>

int          mpfr_cgbn_available (void);        /* 1 if a device is usable   */
const char * mpfr_cgbn_device_name (void);      /* or NULL                   */
mpfr_prec_t  mpfr_cgbn_get_max_prec (void);     /* 65536                     */
void         mpfr_cgbn_free_cache (void);       /* release device resources  */

Set MPFR_CGBN_DEVICE in the environment to choose among several GPUs. When no device is usable, every function below still works and computes on the CPU.

Batches over host arrays

long mpfr_cgbn_mul  (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b,
                     size_t n, mpfr_rnd_t rnd, int *tern);
long mpfr_cgbn_sqr  (mpfr_ptr r, mpfr_srcptr a, size_t n, mpfr_rnd_t, int *);
long mpfr_cgbn_add  (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, size_t, ...);
long mpfr_cgbn_sub  (...);
long mpfr_cgbn_div  (...);
long mpfr_cgbn_sqrt (mpfr_ptr r, mpfr_srcptr a, size_t n, mpfr_rnd_t, int *);
long mpfr_cgbn_fma  (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
                     size_t n, mpfr_rnd_t rnd, int *tern);

The arrays are contiguous, so an mpfr_t v[N] is passed as v[0]:

mpfr_t a[N], b[N], r[N];
int tern[N];
...
mpfr_cgbn_mul (r[0], a[0], b[0], N, MPFR_RNDN, tern);

All elements of r must share one precision; if they do not, or if that precision is above mpfr_cgbn_get_max_prec(), the work is done on the CPU instead. Operands may have any precisions. The return value is how many elements were evaluated on the device (0 when the CPU path was taken), which never affects the results, only tells you how they were obtained.

mpfr_cgbn_fma is a genuine fused multiply-add, like mpfr_fma: the product is formed exactly on the device and only the sum is rounded. With c close to -a*b a multiply followed by an add returns exactly zero, while the fused operation returns the tail of the exact product.

Device-resident vectors

The batch functions pay a PCIe round trip per call. When operations are chained, keeping the values on the device removes that cost entirely, and that is where the large factors are.

mpfr_cgbn_vec_ptr v = mpfr_cgbn_vec_init (n, prec);
mpfr_cgbn_vec_set (v, 0, host_array, n, MPFR_RNDN);      /* upload   */
mpfr_cgbn_vec_mul (r, a, b, MPFR_RNDN);                  /* on device */
mpfr_cgbn_vec_fma (r, a, b, c, MPFR_RNDN);
...
mpfr_cgbn_vec_get (v, 0, host_array, n, MPFR_RNDN, tern);/* download */
mpfr_cgbn_vec_clear (v);

vec_mul, vec_sqr, vec_add, vec_sub, vec_div, vec_sqrt and vec_fma all exist, the vectors must share length and precision, and r may alias a, b or c.

Arithmetic on vectors uses MPFR's extended exponent range, exactly like the body of a scalar MPFR function between MPFR_SAVE_EXPO_MARK and MPFR_SAVE_EXPO_FREE: an intermediate exponent is not checked against emin/emax and no exception flag is raised. Both happen in mpfr_cgbn_vec_get, where each value re-enters the ordinary MPFR world through mpfr_check_range.

Transparent offload of single operations

void        mpfr_cgbn_set_auto_threshold (mpfr_prec_t threshold);
mpfr_prec_t mpfr_cgbn_get_auto_threshold (void);

With a nonzero threshold, mpfr_mul, mpfr_sqr, mpfr_div, mpfr_sqrt and mpfr_fma divert to the device any call whose result and operands all have at least that precision. mpfr_fms and mpfr_hypot call mpfr_fma internally, so they follow it.

This is off by default, and it is not a win. A lone operation has to pay a kernel launch plus two transfers, and measured on the test system:

precision	hook off	hook on
1024	0.18 µs	17.8 µs	offloaded
32768	36.5 µs	116.5 µs	offloaded
65536	96.2 µs	435.6 µs	offloaded
131072	241.6 µs	238.6 µs	declined (above the maximum)

The results are identical in every case; only the time differs. The knob exists for measurement. Use the batch or vector interfaces for actual speedups.

When the threshold is zero the cost to the scalar path is one load and one predictable branch, and it does not show up in measurement — mpfr_mul runs at 7.8 / 8.1 / 24.2 / 117.0 ns at 53 / 113 / 256 / 1024 bits in a --enable-cgbn build, the same as without it.

---

Performance

Test system

GPU	NVIDIA RTX PRO 4500 Blackwell, sm_120, 82 SMs at 2.41 GHz
GPU memory	31.4 GiB, 256-bit bus at 14 GHz (≈896 GB/s), 64 MiB L2
Registers	65536 per SM, 255 per thread (an ISA limit, see below)
CUDA	13.2, driver 595.71.05
MPFR / GMP	4.2.1 / 6.3.0
CPU baseline	one core, the same machine, scalar MPFR

Reproduce with:

%>> cd tests && gmake bench-cgbn         # everything, about 6 minutes
%>> ./cgbnbench div sqrt                 # or just some operations

Three columns are measured, all wall clock, all per operation:

• CPU — the scalar MPFR function in a loop, one core;
• batch — mpfr_cgbn_*(), which packs the operands, sends them over PCIe, runs the kernel and brings the results back;
• resident — the same arithmetic on values already in device memory, which is what a chain of operations costs once the data is there.

The gap between the last two is the transfer, and it is the reason the device-vector interface exists. The last two columns are speedups over the CPU column; higher is better and anything below 1.0x is a loss.

Full timings, all precisions

mpfr_mul / mpfr_cgbn_mul

precision	n	CPU µs/op	batch µs/op	resident µs/op	batch	resident
128	1	0.012	15.506	6.3803	0.0x	0.0x
128	64	0.015	0.298	0.1694	0.0x	0.1x
128	1024	0.020	0.046	0.0081	0.4x	2.5x
128	16384	0.019	0.019	0.0006	1.0x	32.8x
256	1	0.027	15.412	6.4090	0.0x	0.0x
256	64	0.029	0.293	0.1183	0.1x	0.2x
256	1024	0.032	0.060	0.0074	0.5x	4.4x
256	16384	0.033	0.025	0.0006	1.3x	51.9x
512	1	0.055	15.612	6.4793	0.0x	0.0x
512	64	0.063	0.304	0.1205	0.2x	0.5x
512	1024	0.062	0.064	0.0076	1.0x	8.1x
512	16384	0.062	0.034	0.0007	1.8x	90.0x
1024	1	0.120	15.631	6.7815	0.0x	0.0x
1024	64	0.129	0.332	0.1191	0.4x	1.1x
1024	1024	0.128	0.069	0.0077	1.9x	16.6x
1024	16384	0.128	0.057	0.0011	2.3x	117.5x
2048	1	0.361	16.556	7.4993	0.0x	0.0x
2048	64	0.332	0.454	0.1228	0.7x	2.7x
2048	1024	0.340	0.114	0.0093	3.0x	36.5x
2048	16384	0.355	0.118	0.0026	3.0x	134.8x
4096	1	1.159	17.545	8.3906	0.1x	0.1x
4096	64	1.068	0.615	0.1383	1.7x	7.7x
4096	1024	1.124	0.229	0.0131	4.9x	85.6x
4096	16384	1.223	0.265	0.0053	4.6x	230.3x
8192	1	3.623	23.605	14.2229	0.2x	0.3x
8192	64	3.285	0.965	0.2386	3.4x	13.8x
8192	1024	3.478	0.377	0.0290	9.2x	119.7x
8192	16384	3.828	0.566	0.0167	6.8x	229.5x
16384	1	10.332	38.775	28.8575	0.3x	0.4x
16384	64	9.499	1.272	0.4885	7.5x	19.4x
16384	1024	10.079	0.717	0.0789	14.1x	127.8x
16384	16384	10.504	1.017	0.0545	10.3x	192.6x
32768	1	31.471	104.361	93.6821	0.3x	0.3x
32768	64	29.173	2.961	1.5540	9.9x	18.8x
32768	1024	30.902	1.514	0.2770	20.4x	111.5x
32768	16384	31.803	1.958	0.2068	16.2x	153.8x
40960	1	40.927	154.182	143.0526	0.3x	0.3x
40960	64	37.880	4.093	2.3476	9.3x	16.1x
40960	1024	40.556	1.985	0.4107	20.4x	98.8x
40960	16384	41.376	2.451	0.3128	16.9x	132.3x
49152	1	54.777	238.432	226.4823	0.2x	0.2x
49152	64	50.663	5.754	3.6863	8.8x	13.7x
49152	1024	54.063	2.554	0.6189	21.2x	87.4x
49152	16384	55.686	2.967	0.4760	18.8x	117.0x
57344	1	68.039	334.602	322.1672	0.2x	0.2x
57344	64	62.443	7.627	5.1941	8.2x	12.0x
57344	1024	66.575	3.177	0.8506	21.0x	78.3x
57344	16384	68.079	3.645	0.6584	18.7x	103.4x
65536	1	81.976	423.920	411.4637	0.2x	0.2x
65536	64	75.667	9.438	6.6999	8.0x	11.3x
65536	1024	80.755	3.874	1.1538	20.8x	70.0x
65536	16384	82.395	4.415	0.9152	18.7x	90.0x

mpfr_add / mpfr_cgbn_add

precision	n	CPU µs/op	batch µs/op	resident µs/op	batch	resident
128	1	0.010	15.713	6.8178	0.0x	0.0x
128	64	0.009	0.318	0.1521	0.0x	0.1x
128	1024	0.008	0.047	0.0096	0.2x	0.9x
128	16384	0.013	0.023	0.0007	0.5x	18.6x
256	1	0.015	15.852	6.8664	0.0x	0.0x
256	64	0.013	0.305	0.1322	0.0x	0.1x
256	1024	0.013	0.061	0.0083	0.2x	1.6x
256	16384	0.019	0.028	0.0008	0.7x	24.1x
512	1	0.020	15.810	6.8822	0.0x	0.0x
512	64	0.015	0.321	0.1351	0.0x	0.1x
512	1024	0.015	0.067	0.0085	0.2x	1.8x
512	16384	0.021	0.037	0.0008	0.6x	26.3x
1024	1	0.032	15.820	6.8631	0.0x	0.0x
1024	64	0.022	0.337	0.1225	0.1x	0.2x
1024	1024	0.024	0.072	0.0077	0.3x	3.1x
1024	16384	0.028	0.058	0.0009	0.5x	33.1x
2048	1	0.056	15.568	6.5108	0.0x	0.0x
2048	64	0.038	0.443	0.1079	0.1x	0.3x
2048	1024	0.042	0.115	0.0070	0.4x	6.0x
2048	16384	0.045	0.104	0.0009	0.4x	49.0x
4096	1	0.106	15.834	6.7752	0.0x	0.0x
4096	64	0.075	0.592	0.1103	0.1x	0.7x
4096	1024	0.077	0.196	0.0074	0.4x	10.5x
4096	16384	0.102	0.261	0.0011	0.4x	90.5x
8192	1	0.222	16.166	6.7883	0.0x	0.0x
8192	64	0.144	0.838	0.1180	0.2x	1.2x
8192	1024	0.147	0.357	0.0085	0.4x	17.2x
8192	16384	0.228	0.547	0.0019	0.4x	120.1x
16384	1	0.422	17.376	7.3702	0.0x	0.1x
16384	64	0.277	0.949	0.1597	0.3x	1.7x
16384	1024	0.292	0.651	0.0158	0.4x	18.5x
16384	16384	0.498	0.964	0.0070	0.5x	71.0x
32768	1	0.824	22.259	11.5445	0.0x	0.1x
32768	64	0.539	1.695	0.2881	0.3x	1.9x
32768	1024	0.577	1.260	0.0388	0.5x	14.9x
32768	16384	0.890	1.742	0.0219	0.5x	40.6x
40960	1	1.034	21.127	9.6001	0.0x	0.1x
40960	64	0.665	1.977	0.2457	0.3x	2.7x
40960	1024	0.720	1.590	0.0399	0.5x	18.0x
40960	16384	1.085	2.103	0.0268	0.5x	40.4x
49152	1	1.231	22.846	11.0553	0.1x	0.1x
49152	64	0.792	2.329	0.2934	0.3x	2.7x
49152	1024	0.881	1.957	0.0489	0.5x	18.0x
49152	16384	1.279	2.483	0.0322	0.5x	39.7x
57344	1	1.442	24.014	11.5503	0.1x	0.1x
57344	64	0.922	2.683	0.3156	0.3x	2.9x
57344	1024	1.077	2.348	0.0545	0.5x	19.8x
57344	16384	1.481	2.852	0.0370	0.5x	40.0x
65536	1	1.600	31.656	18.9194	0.1x	0.1x
65536	64	1.060	3.165	0.5002	0.3x	2.1x
65536	1024	1.278	2.794	0.0766	0.5x	16.7x
65536	16384	1.697	3.267	0.0478	0.5x	35.5x

mpfr_div / mpfr_cgbn_div

precision	n	CPU µs/op	batch µs/op	resident µs/op	batch	resident
128	1	0.031	16.592	7.7731	0.0x	0.0x
128	64	0.031	0.362	0.1965	0.1x	0.2x
128	1024	0.037	0.049	0.0125	0.7x	3.0x
128	16384	0.042	0.020	0.0009	2.1x	45.0x
256	1	0.059	17.613	8.5245	0.0x	0.0x
256	64	0.058	0.361	0.1886	0.2x	0.3x
256	1024	0.062	0.065	0.0120	1.0x	5.2x
256	16384	0.063	0.026	0.0014	2.4x	44.7x
512	1	0.107	20.459	11.5524	0.0x	0.0x
512	64	0.102	0.551	0.3634	0.2x	0.3x
512	1024	0.104	0.080	0.0228	1.3x	4.6x
512	16384	0.105	0.037	0.0031	2.8x	33.2x
1024	1	0.240	24.534	15.6315	0.0x	0.0x
1024	64	0.226	0.632	0.4138	0.4x	0.5x
1024	1024	0.225	0.089	0.0268	2.5x	8.4x
1024	16384	0.224	0.063	0.0056	3.6x	40.0x
2048	1	0.671	28.004	19.0298	0.0x	0.0x
2048	64	0.622	0.637	0.3006	1.0x	2.1x
2048	1024	0.617	0.127	0.0229	4.9x	27.0x
2048	16384	0.617	0.115	0.0103	5.4x	59.7x
4096	1	2.022	33.132	24.0035	0.1x	0.1x
4096	64	1.876	0.867	0.3817	2.2x	4.9x
4096	1024	1.869	0.223	0.0330	8.4x	56.6x
4096	16384	1.927	0.279	0.0174	6.9x	110.9x
8192	1	6.184	48.556	39.1302	0.1x	0.2x
8192	64	5.937	1.352	0.6187	4.4x	9.6x
8192	1024	5.921	0.411	0.0644	14.4x	92.0x
8192	16384	6.064	0.580	0.0398	10.5x	152.2x
16384	1	18.366	84.088	74.1717	0.2x	0.2x
16384	64	18.112	2.006	1.2056	9.0x	15.0x
16384	1024	18.105	0.827	0.1825	21.9x	99.2x
16384	16384	18.419	1.079	0.1173	17.1x	157.1x
32768	1	54.896	230.791	219.9671	0.2x	0.2x
32768	64	54.379	4.958	3.5323	11.0x	15.4x
32768	1024	54.500	1.792	0.5488	30.4x	99.3x
32768	16384	54.918	2.222	0.4263	24.7x	128.8x
40960	1	77.949	8541.106	8529.5441	0.0x	0.0x
40960	64	77.827	169.954	168.1291	0.5x	0.5x
40960	1024	78.182	36.675	35.0021	2.1x	2.2x
40960	16384	78.556	34.155	23.4965	2.3x	3.3x
49152	1	103.476	11821.067	13585.8539	0.0x	0.0x
49152	64	103.349	273.920	271.8688	0.4x	0.4x
49152	1024	103.949	67.813	61.9392	1.5x	1.7x
49152	16384	104.275	52.917	42.6265	2.0x	2.4x
57344	1	133.738	17602.872	18365.6716	0.0x	0.0x
57344	64	132.710	420.388	426.4542	0.3x	0.3x
57344	1024	133.188	106.911	97.4560	1.2x	1.4x
57344	16384	133.361	90.594	69.3402	1.5x	1.9x
65536	1	168.602	24043.487	26995.4980	0.0x	0.0x
65536	64	168.645	571.694	568.9138	0.3x	0.3x
65536	1024	169.187	167.406	150.0392	1.0x	1.1x
65536	16384	169.543	140.713	115.4445	1.2x	1.5x

mpfr_sqrt / mpfr_cgbn_sqrt

precision	n	CPU µs/op	batch µs/op	resident µs/op	batch	resident
128	1	0.088	15.859	6.8705	0.0x	0.0x
128	64	0.089	0.371	0.2166	0.2x	0.4x
128	1024	0.107	0.042	0.0142	2.5x	7.6x
128	16384	0.114	0.015	0.0011	7.8x	106.3x
256	1	0.122	17.562	8.7010	0.0x	0.0x
256	64	0.117	0.412	0.2451	0.3x	0.5x
256	1024	0.147	0.053	0.0155	2.8x	9.5x
256	16384	0.155	0.019	0.0019	8.2x	81.4x
512	1	0.204	21.061	12.0132	0.0x	0.0x
512	64	0.209	0.561	0.3862	0.4x	0.5x
512	1024	0.244	0.074	0.0248	3.3x	9.8x
512	16384	0.251	0.027	0.0034	9.3x	73.8x
1024	1	0.349	25.831	16.8799	0.0x	0.0x
1024	64	0.365	0.658	0.4614	0.6x	0.8x
1024	1024	0.407	0.075	0.0305	5.4x	13.3x
1024	16384	0.411	0.045	0.0064	9.2x	64.5x
2048	1	0.671	31.097	22.0461	0.0x	0.0x
2048	64	0.682	0.608	0.3637	1.1x	1.9x
2048	1024	0.722	0.100	0.0265	7.2x	27.2x
2048	16384	0.723	0.081	0.0117	9.0x	61.6x
4096	1	1.587	37.422	28.3685	0.0x	0.1x
4096	64	1.551	0.842	0.4600	1.8x	3.4x
4096	1024	1.584	0.172	0.0372	9.2x	42.6x
4096	16384	1.625	0.179	0.0192	9.1x	84.6x
8192	1	4.323	54.122	44.8320	0.1x	0.1x
8192	64	4.270	1.319	0.7441	3.2x	5.7x
8192	1024	4.274	0.309	0.0698	13.8x	61.2x
8192	16384	4.404	0.379	0.0409	11.6x	107.6x
16384	1	12.869	96.581	86.8228	0.1x	0.1x
16384	64	12.883	1.993	1.4099	6.5x	9.1x
16384	1024	12.839	0.609	0.1687	21.1x	76.1x
16384	16384	12.999	0.740	0.1132	17.6x	114.8x
32768	1	39.244	232.850	222.3122	0.2x	0.2x
32768	64	39.736	4.572	3.5777	8.7x	11.1x
32768	1024	39.455	1.398	0.5718	28.2x	69.0x
32768	16384	39.759	1.621	0.4439	24.5x	89.6x
40960	1	55.002	18349.903	18336.0468	0.0x	0.0x
40960	64	55.451	322.089	320.2901	0.2x	0.2x
40960	1024	55.231	81.422	75.6311	0.7x	0.7x
40960	16384	55.577	59.465	54.1651	0.9x	1.0x
49152	1	74.142	26391.501	26317.2796	0.0x	0.0x
49152	64	74.356	459.697	460.5828	0.2x	0.2x
49152	1024	74.151	147.171	132.7680	0.5x	0.6x
49152	16384	74.433	110.315	100.3821	0.7x	0.7x
57344	1	92.622	35670.410	34572.5443	0.0x	0.0x
57344	64	93.490	681.721	645.3638	0.1x	0.1x
57344	1024	93.281	211.995	202.9699	0.4x	0.5x
57344	16384	93.438	179.977	157.0540	0.5x	0.6x
65536	1	117.296	44241.387	45837.8247	0.0x	0.0x
65536	64	117.862	1091.891	1015.8000	0.1x	0.1x
65536	1024	117.281	302.672	295.1039	0.4x	0.4x
65536	16384	117.810	285.203	231.9795	0.4x	0.5x

mpfr_fma / mpfr_cgbn_fma

precision	n	CPU µs/op	batch µs/op	resident µs/op	batch	resident
128	1	0.040	15.681	6.8885	0.0x	0.0x
128	64	0.037	0.343	0.1773	0.1x	0.2x
128	1024	0.038	0.057	0.0123	0.7x	3.1x
128	16384	0.038	0.026	0.0009	1.5x	41.0x
256	1	0.042	16.603	7.6495	0.0x	0.0x
256	64	0.046	0.341	0.1621	0.1x	0.3x
256	1024	0.048	0.071	0.0103	0.7x	4.6x
256	16384	0.047	0.033	0.0011	1.4x	44.9x
512	1	0.071	16.603	7.7560	0.0x	0.0x
512	64	0.080	0.365	0.1673	0.2x	0.5x
512	1024	0.078	0.063	0.0105	1.2x	7.4x
512	16384	0.078	0.046	0.0011	1.7x	69.0x
1024	1	0.182	16.971	7.9783	0.0x	0.0x
1024	64	0.193	0.433	0.1523	0.4x	1.3x
1024	1024	0.189	0.092	0.0100	2.0x	18.8x
1024	16384	0.196	0.086	0.0016	2.3x	123.3x
2048	1	0.525	18.230	9.0814	0.0x	0.1x
2048	64	0.542	0.626	0.1496	0.9x	3.6x
2048	1024	0.536	0.153	0.0106	3.5x	50.6x
2048	16384	0.541	0.164	0.0030	3.3x	180.1x
4096	1	1.621	18.862	9.7597	0.1x	0.2x
4096	64	1.501	0.787	0.1627	1.9x	9.2x
4096	1024	1.575	0.280	0.0150	5.6x	104.7x
4096	16384	1.628	0.438	0.0062	3.7x	260.9x
8192	1	4.864	23.609	13.9656	0.2x	0.3x
8192	64	4.411	0.846	0.2342	5.2x	18.8x
8192	1024	4.642	0.526	0.0289	8.8x	160.4x
8192	16384	4.757	0.728	0.0167	6.5x	284.1x
16384	1	13.397	42.625	32.4762	0.3x	0.4x
16384	64	12.492	1.566	0.5708	8.0x	21.9x
16384	1024	13.214	0.942	0.0836	14.0x	158.0x
16384	16384	13.613	1.310	0.0558	10.4x	244.0x
32768	1	34.991	115.275	104.1351	0.3x	0.3x
32768	64	32.656	3.603	1.7383	9.1x	18.8x
32768	1024	34.619	1.986	0.3008	17.4x	115.1x
32768	16384	35.594	2.494	0.2244	14.3x	158.6x
40960	1	48.272	174.598	162.8546	0.3x	0.3x
40960	64	44.768	4.930	2.6512	9.1x	16.9x
40960	1024	47.766	2.604	0.4446	18.3x	107.4x
40960	16384	48.818	3.186	0.3444	15.3x	141.8x
49152	1	61.866	248.973	240.9318	0.2x	0.3x
49152	64	57.710	6.538	3.8701	8.8x	14.9x
49152	1024	61.281	3.350	0.6469	18.3x	94.7x
49152	16384	62.366	3.908	0.5006	16.0x	124.6x
57344	1	75.777	348.892	336.0542	0.2x	0.2x
57344	64	70.179	8.552	5.4477	8.2x	12.9x
57344	1024	75.057	4.112	0.8872	18.3x	84.6x
57344	16384	76.471	4.518	0.6897	16.9x	110.9x
65536	1	93.811	435.352	420.0096	0.2x	0.2x
65536	64	87.181	10.429	6.9569	8.4x	12.5x
65536	1024	92.819	5.382	1.1814	17.2x	78.6x
65536	16384	94.388	5.410	0.9264	17.4x	101.9x

How to read them

A single operation never pays. Every n = 1 row is a loss, from 15 µs at the low end to 44 ms for a 65536-bit square root. The floor is the kernel launch plus the synchronisation, about 6.5 µs resident and 16 µs with transfers, and no precision CGBN can represent takes the CPU that long. This is why the auto-threshold ships disabled.

Batch crosses over between n = 64 and n = 1024, earlier at higher precision and for the more expensive operations. Peak batch figures: div 30.4x and sqrt 28.2x at 32768 bits, mul 21.2x at 49152, fma 18.3x at 40960, all at n = 1024.

Resident is where the backend earns its keep, because it does no transfer at all: fma 284x at 8192 bits, mul 230x at 4096, div 157x at 16384, add 120x at 8192, sqrt 115x at 16384, all at n = 16384.

Addition never wins on transfer. Every add cell in the batch column is ≤ 0.7x. A CPU addition costs 0.02–1.9 µs and a PCIe round trip cannot be repaid at that price. Its 120x resident column is the whole story for add and sub: offload them between other operations, never on their own.

Peak resident speedup declines past 8192–16384 bits. At the wide end one instance needs a full warp, so 16384 instances no longer fit the 82 SMs in a single wave, and each instance is doing O(n²) schoolbook work. The sweet spot for this backend is 4096–16384 bits with batches of a few thousand.

Division and square root fall off a cliff above 32768 bits. That is the subject of the next section.

---

Why performance varies with bit width

The width ladder

CGBN fixes the width of its arithmetic at compile time, so the backend carries a ladder of widths and rounds each precision up to the next rung:

128  256  512  1024  2048  4096  8192  16384  32768  40960  49152  57344  65536

Powers of two up to 32768, then steps of 8192. Cost follows the field a precision lands in, not the precision asked for. A 33000-bit operand does the work of a 40960-bit one; one bit more than 40960 moves it to 49152.

Each rung must satisfy two constraints:

• (bits/32) % TPI == 0, so CGBN inserts no padding limbs. The kernels rely on this: it is what lets them treat cgbn_clz() == BITS as "value is zero" and read a cgbn_wide_t as a flat 2*BITS-bit integer.
• bits % GMP_NUMB_BITS == 0, so a significand of exactly that precision fits (mpfr_cgbn_pack asserts words >= xn*WPL).

With TPI = 32 that means any multiple of 1024 bits.

The cliff, and what causes it

CGBN chooses the reciprocal algorithm behind its divide and square root from the ratio of limbs per thread to threads per instance (core/core.cu):

DLIMBS_ALG = (LIMBS <= TPI/2) ? dlimbs_algs_half
           : (LIMBS <= TPI)   ? dlimbs_algs_full
                              : dlimbs_algs_multi;

LIMBS = bits/(32*TPI), so the fast paths need bits <= 32*TPI². TPI cannot exceed a warp, so the boundary is exactly 32768 bits — which is precisely the maximum size CGBN documents. It is a real algorithmic boundary, not a conservative one.

Measured per instance at TPI = 32 (tests/cgbn-tune), 512 instances:

bits	limbs/thread	mul µs	div µs	sqrt µs
16384	16	0.061	0.144	0.155
24576	24	0.127	0.268	0.282
32768	32	0.225	0.462	0.460
40960	40	0.332	32.68	74.91
49152	48	0.472	56.90	138.94
57344	56	0.625	89.34	214.93
65536	64	0.855	156.34	319.17

One step of 8192 bits — a 25% wider operand — costs division 71x and square root 163x, while multiplication pays the 1.5x its extra width accounts for and keeps scaling quadratically right through. That is the moment LIMBS exceeds TPI.

It is not register pressure

This matters for the hardware question below, so it is worth establishing carefully. Registers and spill per kernel, from cuobjdump -res-usage on the built objects:

bits	tpi	limbs/thr	mul reg (spill)	div reg (spill)	sqrt reg (spill)	fma reg (spill)
128	4	1	35 (none)	40 (none)	36 (none)	41 (none)
256	8	1	42 (none)	40 (none)	36 (none)	44 (none)
512	8	2	48 (none)	47 (none)	45 (none)	56 (none)
1024	16	2	46 (none)	47 (none)	45 (none)	52 (none)
2048	32	2	44 (none)	48 (none)	40 (none)	47 (none)
4096	32	4	52 (none)	64 (none)	48 (none)	64 (none)
8192	32	8	72 (none)	88 (none)	72 (none)	96 (none)
16384	32	16	128 (none)	150 (none)	112 (none)	168 (none)
32768	32	32	224 (none)	255 (72 B)	214 (none)	255 (24 B)
40960	32	40	255 (72 B)	255 (1328 B)	255 (1312 B)	255 (184 B)
49152	32	48	255 (168 B)	255 (1776 B)	255 (1744 B)	255 (544 B)
57344	32	56	255 (416 B)	255 (2336 B)	255 (2304 B)	255 (936 B)
65536	32	64	255 (760 B)	255 (2880 B)	255 (2944 B)	255 (1136 B)

Spilling is real above 32768 bits. It is not what causes the cliff:

1. At 32768 bits div is already at 255 registers and already spilling 72 B — and it runs at 0.46 µs. High register usage plus spilling is not by itself slow.
2. The time jumps discontinuously while the spill grows smoothly. From 32768 to 40960 the div spill grows 18x (72 → 1328 B) but the time grows 71x. From 40960 to 65536 the spill grows 2.2x while the time grows 4.8x. If spilling drove the cost, time would track spill volume; it does not.
3. mul at 65536 spills 760 B and pays essentially nothing for it. Its cost from 32768 to 65536 goes 0.225 → 0.855 µs, a factor 3.8 against the 4.0 that doubling a schoolbook multiply costs anyway. So ~760 B of spill is worth a few percent, not a factor of 70.
4. The discontinuity lands exactly at LIMBS > TPI, the documented dispatch above.

Register pressure costs a few percent. The algorithm change costs two orders of magnitude.

---

CGBN tuning parameters

CGBN takes a parameter class (cgbn_default_parameters_t) with four knobs, and a threads-per-instance count. All were swept at 65536 bits, 512 instances (tests/cgbn-tune):

configuration	mul µs	div µs	sqrt µs
TPI 4 (512 limbs/thread)	38.44	11121	26600
TPI 8 (256 limbs/thread)	15.30	2948	6499
TPI 16 (128 limbs/thread)	3.61	579	1455
TPI 32 (64 limbs/thread)	0.854	156.3	317.6
TPI 32, SHM_LIMIT 32768	0.854	156.4	321.9
TPI 32, MAX_ROTATION 8	0.854	156.5	321.1

TPI — threads per instance

The only knob that matters, and it is already at its optimum everywhere.

CGBN accepts 4, 8, 16 or 32 and validates it (tpi != 32 && tpi != 16 && tpi != 8 && tpi != 4 reports cgbn_unsupported_threads_per_instance). 32 is the ceiling because an instance is a cooperative group inside one warp.

LIMBS = bits/(32*TPI), so raising TPI is doubly good: it reduces the registers each thread holds and it raises the LIMBS <= TPI threshold that keeps divide and square root on their fast path. Lowering TPI is doubly bad, which is why the sweep degrades so sharply: at 65536 bits, TPI 16 is 3.7x worse than TPI 32 and TPI 4 is 71x worse.

The only reason to use less than 32 is at small widths, where a whole warp per instance would leave most lanes with nothing to do and the cross-lane shuffles would dominate. The recommended value is the smallest TPI that keeps at least one 32-bit limb per thread, capped at 32:

Recommended values per rung

These are what the ladder uses (MPFR_CGBN_SIZES_DECL in src/mpfr-cgbn-impl.h), and the sweeps confirm each one:

rung (bits)	TPI	limbs/thread	SHM_LIMIT	MAX_ROTATION	TPB	CONSTANT_TIME
128	4	1	0	4	0	false
256	8	1	0	4	0	false
512	8	2	0	4	0	false
1024	16	2	0	4	0	false
2048	32	2	0	4	0	false
4096	32	4	0	4	0	false
8192	32	8	0	4	0	false
16384	32	16	0	4	0	false
32768	32	32	0	4	0	false
40960	32	40	0	4	0	false
49152	32	48	0	4	0	false
57344	32	56	0	4	0	false
65536	32	64	0	4	0	false

In other words: TPI 32 everywhere from 2048 bits up, CGBN's defaults for everything else. There is no width at which any of the other parameters is worth changing on this hardware.

The other parameters, and why they are left alone

Parameter	Default	Verdict
SHM_LIMIT	0	Lets CGBN stage cross-lane rotations through shared memory instead of registers, which sounds like exactly the fix for the spilling above 32768 bits. Measured: no effect at all (156.4 µs against 156.3). The spill is in the arithmetic's working set, not in the rotation path. Setting it also requires a matching TPB, which would fix the block size. Left at 0.
MAX_ROTATION	4	Controls how far CGBN unrolls its rotation network. Tried at 2 and 8: no effect (156.5 and 144.5 µs, within noise). Left at 4.
TPB	0	Threads per block, used only to validate blockDim.x when error checking is on. The backend launches with 128 threads per block, a multiple of every supported TPI. Leaving it 0 disables a check that costs nothing to skip, since the launcher chooses the block size itself.
CONSTANT_TIME	false	CGBN reports cgbn_unsupported_operation for it. Not applicable.

The backend also uses cgbn_no_checks rather than cgbn_report_monitor: every precondition CGBN would verify (divisor nonzero, num._high < denom for the wide divide, no padding) is established by construction in the kernels, and the monitor costs a branch per operation.

Adding a rung

One more .cu file and about a minute of nvcc. Satisfy both constraints above, add the width to MPFR_CGBN_SIZES_DECL and MPFR_CGBN_NSIZES, declare mpfr_cgbn_launch_<bits>, add the dispatch case in mpfr-cgbn-cuda.cu, add the file to src/Makefile.am (twice — EXTRA_DIST and MPFR_CGBN_CU), and add the width and width+1 to big_precs[] in tests/tcgbn.c.

Add rungs where they pay. Below 32768 bits the cost curve is gentle and powers of two are fine. Above it the curve is steep: rounding a 40000-bit operand up to 65536, as a bare doubling would, costs 2.9x on a multiply and 5.0x on a divide compared with rounding it up to 40960.

---

Would a bigger GPU help?

The question worth answering precisely is: is the performance loss above 32768 bits caused by register spills and register pressure, such that a GPU with more registers, more CUDA cores and more RAM would fix it?

No, on all three counts. Do not buy hardware to fix this.

Registers per thread are fixed by the ISA, not by the GPU

NVIDIA's own table of per-compute-capability limits gives, identically for compute capability 7.5, 8.0, 8.6, 8.7, 8.9, 9.0, 10.x, 11.0 and 12.x:

	value
Maximum number of 32-bit registers per thread	255
Number of 32-bit registers per SM	64 K
Maximum local memory per thread	512 KB

The 255-register ceiling is a property of the instruction set, not of the chip. Every NVIDIA GPU from Turing to Blackwell has exactly the same budget, and so does the next one. A "bigger" GPU has more SMs, each with the same 64 K register file and the same 255-per-thread cap. There is no GPU on which these kernels would spill less.

The loss is not caused by spilling anyway

Established in the section above: at 32768 bits div already runs at 255 registers with 72 B of spill and takes 0.46 µs; at 40960 bits it takes 32.7 µs. The spill grew 18x and the time grew 71x, and mul at 65536 bits spills 760 B while paying only the factor its extra width accounts for. The discontinuity is CGBN switching divide and square root to dlimbs_algs_multi when LIMBS exceeds TPI. More registers, if they existed, would not move that boundary — it is a compile-time algorithm selection driven by the ratio of limbs to threads.

More RAM would not help either

The working set is tiny. One 65536-bit value is 8 KiB; the largest configuration benchmarked here, 16384 elements across four resident vectors at 65536 bits, is about 0.5 GiB against the 31.4 GiB on this card. Host staging is capped at 32 MiB by design. Memory capacity is nowhere near a constraint, and nothing in the backend is bandwidth-bound at these sizes — a 65536-bit multiply moves 24 KiB and spends 0.85 µs on it, which is 28 GB/s against 896 available.

What a bigger GPU would actually buy

More SMs, and that is worth something — just not for this problem:

• Aggregate throughput scales roughly with SM count for the operations that are already parallel: multiply, fused multiply-add, add, sub, and divide and square root below 32768 bits. At large batch sizes this card is SM-limited, so twice the SMs is close to twice the elements per second.
• It does not change any crossover. The n = 1 floor is launch latency, not throughput, so single operations stay a loss. The batch-versus-resident gap is PCIe, which an SM count does not touch.
• It does not touch the cliff. Divide and square root above 32768 bits are serial work inside one instance; more instances in flight does not make one instance faster.

There is one genuine register effect worth naming, and it is also not purchasable: at 255 registers per thread only 256 threads — 8 warps of the 48 an SM can hold — are resident, so occupancy at the widest rungs is about 17% and there is little left to hide latency with. This is part of why peak resident speedup falls from 230x at 4096 bits to 90x at 65536. But the per-SM register file is 64 K on every current capability too, so no GPU improves it.

What actually fixes it

In decreasing order of practicality:

1. Stay at or below 32768 bits for division and square root. This is the single most effective change available and costs nothing: at 32768 bits div is 30.4x and sqrt 28.2x faster than the CPU in batch; one rung higher they are 2.3x and 0.9x.
2. Use the intermediate rungs. If a precision above 32768 is unavoidable, 40960 and 49152 are 4.9x and 2.7x better for division than letting it round up to 65536.
3. Keep the data resident. The transfer, not the arithmetic, is what caps the batch column at a single-digit factor for the cheap operations.
4. Replace CGBN's wide divide and square root with Newton–Raphson or Barrett reduction built on cgbn_mul_wide, which is the operation that scales cleanly past the boundary. This is a real project, but it is where the factor of 70 is, and it is a software change rather than a hardware purchase.

---

Implementation notes

Everything is behind #ifdef MPFR_WANT_CGBN; a build without --enable-cgbn compiles none of it and exports nothing.

Files

File	Role
src/mpfr-cgbn.h	The installed public header: batch functions, device vectors, the auto-threshold knob.
src/mpfr-cgbn-impl.h	The C/C++ ABI shared with the device code: packed layout, launch descriptor, width ladder, the thin CUDA runtime wrapper. Must stay in the common subset of C and C++.
src/mpfr-cgbn.c	Host glue: mpfr_t ⇄ device layout, chunked staging through page-locked buffers, the batch driver, the CPU fallback, the hook entry points. Plain C, includes mpfr-impl.h.
src/mpfr-cgbn-kernels.h	All the device arithmetic, templated on the CGBN width. Not compiled on its own.
src/mpfr-cgbn-s<BITS>.cu	One per rung; defines MPFR_CGBN_BITS / MPFR_CGBN_TPI and includes the kernels. Separate translation units so make -j parallelises them.
src/mpfr-cgbn-cuda.cu	The only file that touches the CUDA runtime API.
tests/tcgbn.c	The correctness test.
examples/mpfr-cgbn-example.c	Worked example of every operation.

The five hook sites are in src/mul.c, sqr.c, div.c, sqrt.c and fma.c, each a single #ifdef-guarded branch just after MPFR_LOG_FUNC.

How a number is represented on the device

A batch is a structure of two arrays: a significand array of n * (bits/32) little-endian 32-bit words, exactly the layout of cgbn_mem_t<bits>, and a metadata array of {int64 exp, int32 sign, int32 kind}.

The significand is left-aligned in the field, so a regular number has bit bits-1 set and its value is sign * 0.significand * 2^exp — MPFR's own convention, zero-padded on the right out to the field width. A number of precision p has its low bits-p bits zero.

Special values are a separate kind tag rather than an encoded exponent as in mpfr_t, which keeps the device hot path branch-free.

Things that are easy to get wrong

These are the traps the implementation had to navigate; they are worth knowing before modifying it.

• Division and square root produce their significand exactly but only one further bit of the infinite tail, and that bit lives in the remainder, not in the quotient limbs. regular_div and regular_sqrt extract it and wide_set_round_bit parks it where round_wide expects the first bit below the significand. Collapsing it into the sticky bit gives wrong results whenever the target precision equals the field width — it produced an off-by-one-ulp result with an inverted ternary value for a/b with b = a + 1 ulp.
• The fused multiply-add needs a 3*BITS accumulator (triple_t), where everything else rounds out of a 2*BITS wide value. Its exact product is already 2*BITS, and the addend can either sit a full field below it or cancel it down to its last bit. Anything falling out of that window is provably below the round bit, so it only feeds the sticky bit — that proof is what bounds the accumulator, so re-derive it before changing the placement.
• MPFR raises its NaN flag for every NaN result, including the plain propagation of a NaN operand (see MPFR_RET_NAN). This is not IEEE 754 behaviour, and the kernel does it in one place at the store: label.
• Device arithmetic runs in the extended exponent range, like the body of a scalar MPFR function between MPFR_SAVE_EXPO_MARK and MPFR_SAVE_EXPO_FREE. emin/emax and the exception flags are applied on the host, in mpfr_cgbn_unpack, through mpfr_check_range.
• The CPU fallback must not be reachable from the hooks. mpfr_cgbn_batch falls back to the scalar functions, which are hooked, which would recurse. The hooks call mpfr_cgbn_device_batch instead, which does only what the device can and returns the count; when it declines, mpfr_mul and friends simply carry on with their own code.
• nvcc rejects the -fPIC libtool would pass a C compiler, so src/Makefile.am writes the .lo descriptors itself. The first four lines of those files must keep matching libtool's ^# Generated by .*libtool, or libtool rejects them as "not a valid libtool object".
• CGBN's host-side fallback is selected by __GMP_H__, so in the host compilation pass cgbn_env_t is the mpz emulation, which lacks members such as PADDING. Device code must stick to the API that exists in both.

Testing

tests/tcgbn is the safety net and compares against the scalar functions only — never against hard-coded expectations. It covers every rounding mode, the width-ladder boundaries and the precisions just past each rung, the special values, exact ties, one-ulp neighbours, total cancellation, in-place aliasing on every operand position, chunked batches, and the device-vector path. Run it with GMP_CHECK_RANDOMIZE=1 after touching the kernels.

check_fma_is_fused guards the multiply-add specifically: it fails if the generated data ever stops distinguishing a fused evaluation from a multiply followed by an add, so that check cannot quietly go vacuous.

ABI

libmpfr's -version-info is deliberately unchanged. Without --enable-cgbn no new symbol is exported, so the default ABI is identical; with it, the added symbols are a local build-time choice, the same convention MPFR already uses for --enable-decimal-float. make check-exported-symbols and make check-gmp-symbols both pass.

---

Tools

All four live in tests/ and are EXTRA_PROGRAMS: built only on demand, never by make check, because they report rather than assert.

Program	Answers	Build time
cgbnbench	How fast is the backend against the scalar functions? Name operations on the command line (./cgbnbench div sqrt) to measure only some. Repetition counts are chosen from a timed pass, so a full run is minutes rather than half an hour.	seconds
cgbn-diag	Which precisions reach the device; what the auto-threshold hook costs; is the fused multiply-add really fused; what does the backend cost the scalar path when idle.	seconds
cgbn-verify	Is CGBN itself correct against GMP at every width in the ladder, and do the primitive semantics the kernels rely on still hold (clz of zero, the mask_and sign convention, shift edges, carry-out)?	~50 s
cgbn-tune	How does cost vary with threads per instance and with width? This is what locates the cliff.	~9 min

%>> cd tests
%>> make bench-cgbn            # cgbnbench, everything
%>> make cgbn-diag  && ./cgbn-diag
%>> make cgbn-verify && ./cgbn-verify
%>> make cgbn-tune   && ./cgbn-tune

cgbnbench and cgbn-diag are ordinary C linked against libmpfr and build without the backend too, printing a notice and exiting. cgbn-verify and cgbn-tune talk to CGBN directly and are compiled by nvcc.

---

Limitations

• Precision above 65536 bits runs on the CPU. The ladder stops there.
• Division and square root above 32768 bits are slower than the CPU on this hardware, for the structural reason set out above. They are still sent to the device; the numbers are published here so the choice is an informed one.
• A single operation is never faster on the device. The auto-threshold exists for measurement.
• Addition and subtraction do not repay a transfer. Use them on resident data, between other operations.
• All results in one batch must share a precision; otherwise the call falls back to the CPU. Operand precisions may differ freely.
• Device access is serialized on an internal mutex, so concurrent calls from several threads are safe but do not overlap on the GPU. Exception flags and the exponent range remain thread-local, as everywhere else in MPFR.
• GMP_NUMB_BITS must be a multiple of 32. Configure checks this.
• MPFR_RNDF (faithful rounding) is evaluated as MPFR_RNDZ, which is a faithful result; its ternary value is not meaningful, as MPFR documents.
• CGBN is used past its documented maximum size. Widths above 32768 bits were added only after checking CGBN's multiply, divide and square root against GMP at each of them, and tests/cgbn-verify re-checks that at any time.