This document describes the use and implementation of the modulars within the numerix package.

1.General modulars and moduli

Moduli are implemented in modulus.hpp. An object of type modulus is declared as follows:

modulus<C,V> p;

V is an implementation variant. The default variant, named modulus_naive, is implemented in modulus_naive.hpp and supports a type C with an Euclidean division. More variants are described below. Modulars are implemented in modular.hpp. An object of type modular can be used as follows:

modular<modulus<C, V>, W> x, y, z;
modular<modulus<C, V>, W>::set_modulus (p);
x = y * z;

The variant W is used to define the storage of the modulus. Default storage is modular_local that allows the current modulus to be changed at any time. The default modulus is 0. For the integer type the default variant is set to modulus_integer_naive, so that modular integers can be used as follows:

#include<numerix/integer.hpp>
#include<numerix/modular.hpp>
#include<numerix/modulus_integer_naive.hpp>

typedef modulus<integer> M;
M p (9973);
modular<M>::set_modulus (p);
modular<M> a (1), b (3);
c = a * b;
cout << c << "\n";

If the modulus is known to be fixed at compilation time then the following declaration is preferable for the sake of efficiency:

fixed_value <int, 3> w;
modulus<integer> p (w);

If the modulus is not intended to be changed and known at compilation time then the variant modular_fixed can be used.

2.Moduli for hardware integers

From now on a div p and a mod p represent the quotient and the remainder in the long division of a by p respectively. Throughout this section C denotes a genuine integer C type. We write n for the bit-size of C, and p for a modulus that fits in C.

2.1.Naive implementation

The default variant, modulus_int_naive<m>, is defined in modular_int.hpp. Each modular product performs one product and one integer division. The parameter m is the maximum bit-size of the modulus allowed by the variant. This parameter can be set to n. Best performances are obtained with m = n - 1. Note that we necessarily have m⩽n - 1 if C is a signed integer type. By convention the modulus 0 actually means 2^m.

2.2.Moduli with pre-computed inverse

In the variant modulus_int_preinverse<m> a suitable inverse of the modulus is pre-computed and stored as an element of type C, this yields faster computations when doing several operations with the same modulus. The behavior of the parameter m is as in the preceding naive variant. The best performances are attained for when m⩽n - 2.

Within the implementation of this variant the following auxiliary quantities are needed:

• a non-negative integer r is defined by 2^{r - 1}< p ≤ 2^r,
• an integer s such that 0 ≤s≤r - 1,
• an integer t such that t≥r and s + t - (r - 1) ≤n,
• an integer q = 2^{s + t} div p, that represents the numerical inverse of p to precision s + t - k.

By convention we let r = n if p is zero. Note that q<2ⁿ, hence fits in C.

2.2.1.Modular product

Let a be an integer such that 0 ≤a<p². The remainder a mod p can be obtained as follows:

1. Decompose a into a₀ and a₁ such that a = a₁2^s + a₀, with a₀<2^s, and compute b = a₁q div 2^t. From a<p² it follows that a₁q < p² q / 2^s≤p2^t. So a₁q has size at most r + t, and we have b<p.

2. Compute c = a - b p. From the definition of q we have:

   0≤

   2s + t

   p

   - q < 1 ⟹0≤ a - a ₁ q p /2 ^t < a ₁ p /2 ^t + a ₀ .

   It follows that 0≤c<p + a₁p/2^t + a₀.

3. Let h = 2^2r /2^{s + t} + 2^s/2^{r - 1}. From a₀<2^s we deduce that a₁p/2^t + a₀≤h p. It thus suffices to substract at most h times p from a to obtain a mod p.

In general we can always take t = r and s = r - 1, so that h = 3. For when r + 1≤n, it is better to take t = r + 1 and s = r - 1 so that h≤2. Finally, for when r + 2≤n one can take t = r + 3 and s = r - 2 so that h≤1. These settings are summarized in the following table:

In these three cases remark that the integer a₁q has size at most 2n.

It is clear that the case r + 2≤n leads to the fastest algorithm since step 3 reduces to one substraction/comparison that can be implemented as min(c,c - p) within type C alone. A special attention must be drawn to when r = 1, that corresponds to p = 2: we must suppose that s = r - 2 = - 1 is indeed computed within an unsigned int type, hence is sufficiently large to ensure a₁ = 0 and c = a, hence the correctness of the algorithm.

2.2.2.Numerical inverse

Let us now turn to the computation of the numerical inverse q of p. The strategy presented here consists in performing one long division in C followed by one Newton iteration. Let w = 2^{r - 1}/p be the numerical inverse of p, normalized so that 1/2 ≤w<1 holds. Let x₀ denote the initial approximation, and x₁ be its Newton iterate x₁ = 2x₀ - p x

Under the assumption that s₀≤n/2, we can compute a suitable x₀ as the floor part of 2^s₀w by the only use of operations within base type C. In other words we calculate q₀ = 2^{s₀ + r - 1} div p = 2^s₀x₀ as follows:

1. If r≤s₀ then q₀ is obtained by means of one division in C.
2. Otherwise r>s₀. We decompose p into p = p₁2^{r - s₀} + p₀, with p₀<2^{r - s₀}. Note that 2^{s₀ - 1}≤p₁≤2^s₀. Within C perform the long division 2^{2s₀ - 1} = a p₁ + b. By multiplying both side of the latter equality by 2^{r - s₀} we obtain that 2^{s₀ + r - 1} = a p - a p₀ + b2^{r - s₀}. From a and b, we easily deduce q₀ since b2^{r - s₀}<p and a p₀<2^{2s₀ - 1}2^{r - s₀}/2^{s₀ - 1} = 2^r.

The computation of q₁ = 2^{s₁ + r - 1} div p, that is the floor part of 2^s₁w, then continues as follows:

1. From the definition of x1, one can write: x12s11 = q02s0 - p q

   2

   0

   /2^r. We thus compute q₀2^s₀ and the ceil part of p q

   2

   0

   /2^r in order to deduce the floor part c of x₁2^s₁. Note that c fits C.
2. From the above inequalities it follows that 0≤2^s₁w - c<2. We thus obtain d = 2^{s₁ + r - 1} - c p < 2p. If d≥p then return q₁ = c + 1 else return q₁ = c.

Finally, q is deduced in this way:

1. Set s₀ as (s + t - (r - 1)) div 2, and compute q₁ = 2^{s₁ + r - 1} div p as just described.
2. Note that 1≤s + t - (r - 1) - s₁≤2. If s₁+1=s + t - (r - 1) then q = 2^{s₁ + r} div p otherwise q = 2^{s₁ + r + 1} div p.

2.2.3.Modular product by a constant

Let b be a modulo p integer that is to be involved in several modular products. Then one can compute β = 2^rb div p < 2^r just once and obtain a speed-up within the products by means of the follows algorithm:

1. Compute α = aβ div 2^r. It follows that 0≤a b - α p < 2p.
2. Compute c = a b - α p. If c≥p then c=c+1. Return c.

The computation of β can be done as follows from q:

1. Compute γ = q b div 2^{t + s - r}<2^r. Note that 2^rb/p - 2^rb/2^{t + s} - 1<γ≤2^rb/p.
2. Since 0≤2^rb - γp<g p, with g = 2 if r + 1≤n and g = 3 otherwise, deduce 2^rb div p from γ.

3.Further remarks

• On some processor architectures code to manipulate short int can be larger and slower than corresponding code which deals with int. This is particularly true on the Intel x86 processors executing 32 bit code. Every instruction which references a short int in such code is one byte larger and usually takes extra processor time to execute.
• Do not use char but signed char or unsigned char for the sake of portability.