🌐

A Test Code Example for Numerical Calculation Using a Residue Number System

As an example of numerical calculation based on a residue number system, this section introduces a program that actually performs the calculations a * b and a2 + b2 within a residue number system capable of handling 18-bit integers. The program is intended for verification of the processing flow in modular arithmetic.

This code is not intended as a tutorial on algorithms or programming techniques.
The reader is assumed to have undergraduate-level knowledge of mathematics, including the Chinese Remainder Theorem (CRT) and the extended Euclidean algorithm, as well as a working understanding of integer arithmetic, arrays, structures, and pointers in C99.

mod_calc -- modular calculation test program (in C99)

/* mod_calc -  modular calculation test program
 * Rev.1.61 (2026-01-04) (c) 2000 Takayuki HOSODA
 * SPDX-License-Identifier: BSD-3-Clause
 */
#include <stdio.h>
#include <inttypes.h>
#define N_VECT 6
typedef struct crt {
    int      n;
    int32_t  vect[N_VECT];
    int32_t  inv[N_VECT];
    int32_t  max;
} crt_t;
int32_t modinv(int32_t a, int32_t max);
int32_t mod2int(crt_t *v, int32_t *m);
void    init_crt(crt_t *v);
void    printv(char *s, int32_t *v);
char    *version = "\
mod_calc -- modular calculation test program\n\
Rev.1.61 (2000-11-01) (c) 2000, Takayuki HOSODA.\n";
char    *help = "\
Synopsys : mod_calc integer integer\n\
Example  : mod_calc 356 412\n";
int32_t modinv(int32_t a, int32_t m) {
    int32_t        q, s, t, u, v, w;
    s = a;
    t = m;
    u = 1;
    v = 0;
    while (s > 0) {
	q = t / s;
	w = t - q * s;
	t = s;
	s = w;
	w = v - q * u;
	v = u;
	u = w;
    }
    return (v < 0 ? v + m : v);
}
void init_crt(crt_t *v) {
    v->n = N_VECT;
    v->max = 1;
    for (int i = 0; i < v->n; i++)
        v->vect[i] = (const int32_t[]){5, 7, 9, 11, 13, 16}[i];
    for (int i = 0; i < v->n; i++)
        v->max *= v->vect[i];
    for (int i = 0; i < v->n; i++) 
	v->inv[i] = ((v->max / v->vect[i]) * modinv(v->max / v->vect[i], v->vect[i]));
}
int32_t mod2int(crt_t *v, int32_t *m) {
    int32_t  y = 0;
    for (int i = 0; i < v->n; i++) 
	y = (y + (m[i] * v->inv[i]) % v->max) % v->max;
    return y;
}
int32_t *int2mod(crt_t *v, int32_t a, int32_t *m) {
    for (int i = 0; i < v->n; i++) 
	m[i] = a % v->vect[i];
    return (m);
}
void printv(char *s, int32_t *m) {
    printf("%s\t", s);
    for (int i = 0; i < N_VECT; i++)
	printf("%"PRId32"\t", m[i]);
    printf("\n");
}
int main(int argc, char *argv[]) {
    crt_t              v;
    int32_t            a, b, y, z;
    int32_t            am[N_VECT];
    int32_t            bm[N_VECT];
    int32_t            zm[N_VECT];
    printf("%s\n", version);
    if (argc == 3) {
        sscanf(&argv[1][0], "%" SCNd32 "\n", &a);
        sscanf(&argv[2][0], "%" SCNd32 "\n", &b);
	init_crt(&v);
	printv("mod[] =", v.vect); 
	printv("inv[] =", v.inv);
	printf("max =\t%"PRId32"\n", v.max);
	int2mod(&v, a, am);
	int2mod(&v, b, bm);
	printf("\na * b calculation test:\n");
	y = (a * b) % v.max;
	printf("ref= %"PRId32" * %"PRId32"\n   = %"PRId32" (mod %"PRId32")\n", a, b, y, v.max);
	for (int i = 0; i < v.n; i++) zm[i] = (am[i] * bm[i]) % v.vect[i];
        for (int i = 0; i < 3; i++) 
            printv((char *[]){"am = ", "bm = ", "zm = "}[i], (int32_t *[]){am, bm, zm}[i]);
	z = mod2int(&v, zm);
	printf("z  = %"PRId32" --%s.\n", z, y == z ? "OK" : "NG");
	printf("\na^2 + b^2 calculation test:\n");
	y = (a * a + b * b) % v.max;
	printf("ref= %"PRId32"^2 + %"PRId32"^2\n   = %"PRId32" (mod %"PRId32")\n", a, b, y, v.max);
	for (int i = 0; i < v.n; i++)
	    zm[i] = ( ((am[i] * am[i]) % v.vect[i])
		    + ((bm[i] * bm[i]) % v.vect[i])) % v.vect[i];
        for (int i = 0; i < 3; i++) 
            printv((char *[]){"am = ", "bm = ", "zm = "}[i], (int32_t *[]){am, bm, zm}[i]);
	z = mod2int(&v, zm);
	printf("z  = %"PRId32" --%s.\n", z, y == z ? "OK" : "NG");
    } else {
        printf("%s", help);
    }
    return 0;
}

The set of pairwise coprime moduli {5, 7, 9, 11, 13, 16} was chosen so that their product exceeds 218 = 262144, the maximum value of an 18-bit unsigned integer, while keeping both the bit width and the number of moduli as small as possible. In this case, 5 �~ 7 �~ 9 �~ 11 �~ 13 �~ 16 = 720720 > 262144, which satisfies the requirement with sufficient margin. Alternatively, one could restrict the set to prime numbers; for example, choosing {19, 23, 29, 31} yields 19 �~ 23 �~ 31 = 392863 > 262144, which would also be acceptable.

Compilation

gcc -O -std=c99 mod_calc.c -o mod_calc

This code uses C99 features such as variable declarations within for statements and compound literals; therefore, a compiler that supports C99 is required.
When using GCC, version 5.0.0 or later is recommended.

Example of Execution Output

mod_calc -- modular calculation test program
Rev.1.61 (2000-11-01) (c) 2000, Takayuki HOSODA.

mod[] =	5	7	9	11	13	16	
inv[] =	576576	205920	320320	196560	277200	585585	
max =	720720

a * b calculation test:
ref= 356 * 412
   = 146672 (mod 720720)
am = 	1	6	5	4	5	4	
bm = 	2	6	7	5	9	12	
zm = 	2	1	8	9	6	0	
z  = 146672 --OK.

a^2 + b^2 calculation test:
ref= 356^2 + 412^2
   = 296480 (mod 720720)
am = 	1	6	5	4	5	4	
bm = 	2	6	7	5	9	12	
zm = 	0	2	2	8	2	0	
z  = 296480 --OK.

SEE ALSO

• An Example of Numerical Calculation Using Modular Arithmetic