Faster Chinese remaindering

 Joris van der Hoeven

 Laboratoire d'informatique, UMR 7161 CNRS Campus de l'École polytechnique 1, rue Honoré d'Estienne d'Orves Bâtiment Alan Turing, CS35003 91120 Palaiseau

 December 16, 2016

 The Chinese remainder theorem is a key tool for the design of efficient multi-modular algorithms. In this paper, we study the case when the moduli are fixed and can even be chosen by the user. Through an appropriate use of the technique of FFT-trading, we will show that this assumption allows for the gain of an asymptotic factor in the complexity of “Chinese remaindering”. For small , we will also show how to choose “gentle moduli” that allow for further gains at the other end. The multiplication of integer matrices is one typical application where we expect practical gains for various common matrix dimensions and integer bitsizes. Keywords: Chinese remainder theorem, algorithm, complexity, integer matrix multiplication

1.Introduction

Modular reduction is an important tool in computer algebra and elsewhere for speeding up computations. The technique allows to reduce a problem that involves large integer or polynomial coefficients to one or more similar problems that only involve small modular coefficients. Depending on the application, the solution to the initial problem is reconstructed via the Chinese remainder theorem or Hensel's lemma. We refer to [10, chapter 5] for a gentle introduction to this topic.

In this paper, we will mainly be concerned with multi-modular algorithms over the integers that rely on the Chinese remainder theorem. The archetype of such an algorithm works as follows. We start with a polynomial function . For any modulus , reduction of modulo yields a new function such that

for all . Given an algorithm to compute that only uses ring operations on integers, it suffices to replace each ring operations by its reduction modulo in order to obtain an algorithm that computes . Now given integers and , assume that we know a bound with for and for . Then the following multi-modular algorithm provides with an alternative way to compute :

Select moduli with that are mutually coprime.

1. For , compute for .

2. For , compute .

3. For , reconstruct from the values with .

Step 1 consists of multi-modular reductions (finding the as a function of ) and step 3 of multi-modular reconstructions (finding as a function the ); this is where the Chinese remainder theorem comes in. For a more detailed example with an application to integer matrix multiplication, we refer to section 4.5.

In favourable cases, the cost of steps 0, 1 and 3 is negligible with respect to the cost of step 2. In such situations, the multi-modular algorithm to compute is usually much faster than the original algorithm. In less favourable cases, the cost of steps 1 and 3 can no longer be neglected. This raises the question whether it is possible to reduce the cost of these steps as much as possible.

Two observations are crucial here. First of all, the moduli are the same for all multi-modular reductions and multi-modular reconstructions. If is large, then this means that we can essentially assume that were fixed once and for all. Secondly, we are free to choose in any way that suits us. By making each fit into a machine word, one may ensure that every modular operation only takes a few cycles. Special “FFT-moduli” are often used as well for speeding up polynomial arithmetic.

In this paper, we will show how to exploit both of the above observations. For fixed moduli, we will show in section 3 how to make Chinese remaindering asymptotically more efficient by a factor when gets large. In section 4, we show that it is possible to construct “gentle modulo” that allow for speed-ups when is small (). Both results can be combined in order to accelerate Chinese remaindering for all possible values of .

The new asymptotic complexity bounds make heavy use of discrete Fourier transforms. For our purposes, it is crucial to avoid “synthetic” FFT schemes that require the adjunction of artificial roots of unity as in Schönhage–Strassen multiplication [24]. Instead, one should use “inborn” FFT schemes that work with approximate roots of unity in or roots of unity with high smooth orders in finite fields; see [24, section 3] and [23, 15]. Basic complexity properties of integer multiplication and division based on fast Fourier techniques are recalled in section 2.

Let be the bit complexity for multiplying two -bit numbers. Given pairwise comprime moduli of bit-size , it is well known that multi-modular reduction and reconstruction can be carried out in time using so called remainder trees [8, 20, 3]. Recent improvements of this technique can be found in [4, 2]. The main goal of section 3 is to show that this complexity essentially drops down to in the case when all moduli are fixed; see Theorems 6 and 10 for more precise statements. The main idea is to increase the arities of nodes in the remainder tree, while performing the bulk of the computations at each node using Fourier representations. This technique of trading faster algorithms against faster representations was also used in [16], where we called it FFT-trading; see also [1]. The same approach can also be applied to the problem of base conversion (see section 3.8) and for univariate polynomials instead of integers (see section 3.9).

Having obtained a non trivial asymptotic speed-up for large , we next turn our attention to the case when is small (say ). The main goal of section 4 there is to exhibit the existence of gentle moduli for which Chinese remaindering becomes more efficient than usual. The first idea is to pick moduli of the form , where is somewhat smaller than the hardware word size, is even, and . In section 4.1, we will show that multi-modular reduction and reconstruction both become a lot simpler for such moduli. Secondly, each can be factored as and, if we are lucky, then both and can be factored into moduli that fit into machine words. If we are very lucky, then this allows us to obtain moduli of bitsize that are mutually coprime and for which Chinese remaindering can be implemented efficiently. Gentle moduli can be regarded as the integer analogue of “special sets of points” that allowed for speed-ups of multi-point evaluation and interpolation in [5].

Acknowledgments. We would like to thank Grégoire Lecerf for pointing us to Bernstein's work [2] on the topic of this paper.

2.Preliminaries

2.1.Integer multiplication

Throughout this paper we will assuming the deterministic multitape Turing model [21] in order to analyze the “bit complexity” of our algorithms. We will denote by the cost of -bit integer multiplication. The best current bound [15] for is

where is called the iterator of the logarithm.

For large , it is well known that the fastest algorithms for integer multiplication [23, 24, 9, 15] are all based on the discrete Fourier transform [7]: denoting by the cost of a “suitable Fourier transform” of bitsize and by the cost of the “inner multiplications” for this bitsize, one has

For the best current algorithm from [15], we have

One always has . The actual size of Fourier transforms is usually somewhat restricted: for efficiency reasons, it should be the product of powers of small prime numbers only, such as , and . Fortunately, for large numbers , it is always possible to find with and .

It is also well known that fast Fourier transforms allow for several tricks. For instance, if one of the multiplicands of an -bit integer product is fixed, then its Fourier transform can be precomputed. This means that the cost of the multiplication drops down to

In particular, the complexity of multiplying an -bit integer with an -bit one (for ) satisfies

Squares of -bit numbers can be computed in time for the same reason. Yet another example is the multiplication of two matrices with -bit integer entries: such multiplications can be done in time by transforming the input matrices, multiplying the transformed matrices in the “Fourier model”, and then transforming the result back.

In the remainder of this paper, we will systematically assume that asymptotically fast integer multiplication is based on fast Fourier transforms. In particular, we have (1) for certain functions and . We will also assume that the functions and are (not necessarily strictly) increasing and that . These additional conditions are satisfied for (2) and (3). The first condition is violated whenever the FFT scheme requires the adjunction of artificial roots of unity. This happens for Schönhage–Strassen multiplication, in which case we have and ). We will say that an FFT-scheme is inborn if it satisfies our additional requirements.

For a hypothetical integer multiplication that runs in time , but for which is (not necessarily strictly) increasing, we also notice that it is possible to design an inborn FFT-based integer multiplication method that runs in time ; this is for instance used in [13].

2.2.Euclidean division of integers

Let denote the cost of euclidean division with remainder of an -bit integer by an -bit one. In [16, section 3.2], we gave an algorithm divide for the euclidean division of a polynomial of degree by another polynomial of degree . This algorithm is based on FFT trading, a technique that consists of doing as much work as possible in the FFT-model even at the expense of using naive, suboptimal algorithms in the FFT-model.

The straightforward adaptation of this division to integer arithmetic yields the asymptotic bound

Furthermore, the discrete Fourier transforms for the dividend make up for roughly one fifth of the total amount of transforms. For , the cost of the transforms for the dividend does not grow with , which leads to the refinement

Similarly, if , then

since the bulk of the computation consists of multiplying the approximate -bit quotient with the -bit dividend. If the dividend is fixed, then we even get

since the Fourier transforms for the dividend can be precomputed.

2.3.Approximate products modulo one

Let us start with a few definitions and notations. Given and , we define

be the set of dyadic fixed point numbers of bitsize and with exponent . Given and , we denote by

the remainder of the euclidean division of by . Given and , we say that is an -approximation of if . We also say that is a circular -approximation of if . Here denotes the circular distance between and .

Let , and . Mutatis mutandis, Bernstein observed in [2] that we may compute a circular -approximation for as follows. We first compute the product of and modulo and then take .

Let us show that indeed satisfies the required property. By construction, there exists an with such that . Therefore, and , whence .

More generally, if we only have a circular -approximation of a number (instead of a number as above), then the algorithm computes a circular -approximation of .

Bernstein's trick is illustrated in Figure 1: we are only interested in the highlighted portion of the product. We notice that these kinds of truncated products are reminiscent of the “middle product” in the case of polynomials [12]; in our setting, we also have to cope with carries and rounding errors. When using FFT-multiplication, products modulo can be computed using three discrete Fourier transforms of bitsize , so the cost is essentially the same as the cost of an -bit integer multiplication. If one of the arguments is fixed, then the cost becomes asymptotic to .

 Figure 1. Product modulo one of and with significant bits.

More generally, for , let be a circular -approximation of a number and let . Then we may compute a circular approximation of as follows. Consider the expansions

 (4)

with ; see Figure 2. By what precedes, for , we may compute circular -approximations for

Setting

it follows that , whereas . Taking , it follows that and

When using FFT-multiplication, we notice that the sum can be computed in the FFT-model, before being transformed back. In that way, we only need instead of transforms of size , for a total cost of .

 Figure 2. Product modulo one of and with significant bits.

2.4.Machine arithmetic

For actual machine implementations of large integer arithmetic, it is customary to choose a base of the form and to perform all computations with respect to that base. We will call the soft word size. For processors that have good hardware support for integer arithmetic, taking or is usually most efficient. The GMP package [11] uses this approach.

However, modern processors are often better at floating point arithmetic. General purpose processors generally provide double precision IEEE-768 compliant instructions, whereas GPUs are frequently limited to single precision. The respective choices and are most adequate in these cases. It is good to pick slightly below the maximum bitsize of the mantissa in order to accelerate carry handling. We refer to [17] for more details on the implementation of multiple precision arithmetic based on floating point arithmetic.

Another particularity of recent processors is the availability of ever wider SIMD (Single Instruction Multiple Data) instructions. For modern implementations of large integer arithmetic, we therefore recommend to focus on the problem of multiplying several () large integers of a given bitsize instead of a single one. We again refer to [17] for more details.

In what follows, when using integer arithmetic, we will denote by the maximal bitsize such that we have a hardware instruction to multiply two integers of bits (e.g. or ). When using floating point arithmetic, we let be the bitsize of a mantissa (i.e. or ). We will call the machine word size. For implementations of multiple precision arithmetic, we always have , but it can be interesting to take .

For moduli that fit into a machine word, arithmetic modulo can be implemented efficiently using hardware instructions. In this case, the available algorithms again tend to be somewhat faster if the size of is a few bit smaller than . We refer to [19] for a survey of the best currently available algorithms in various cases and how to exploit SIMD instructions.

When using arithmetic modulo , it is often possible to delay the reductions modulo as much as possible. One typical example is modular matrix multiplication. Assume that we have two matrices with coefficients modulo (represented by integers between and , say). If fits into a machine word, then we may multiply the matrices using integer or floating point arithmetic and reduce the result modulo . This has the advantage that we may use standard, highly optimized implementations of matrix multiplication. One drawback is that the intermediate results before reduction require at least twice as much space. Also, the bitsize of the modulus is at least twice as small as .

3.Asymptotically fast Chinese remaindering

3.1.The Chinese remainder theorem

For any integer , we will write . We recall:

Chinese Remainder Theorem. Let be positive integers that are mutually coprime and denote . Given , there exists a unique with for .

We will prove the following more constructive version of this theorem.

Theorem 1. Let be positive integers that are mutually coprime and denote . There exist such that for any , the number

satisfies for .

Notation. We call the cofactors for in and also denote these numbers by .

Proof. If , then it suffices to take . If , then the extended Euclidean algorithm allows us to compute a Bezout relation

where and . Let us show that we may take

Indeed, given and , we have

In particular, divides . Since (5) implies , it follows that . Similarly, .

For , we will use induction. Let , and . By induction, we may compute , , and . We claim that we may take

Indeed, for , we get

whence . For we obtain in a similar way.

3.2.Naive multi-modular reduction and reconstruction

Let , , and be as in the Chinese remainder theorem. We will refer to the computation of as a function of as the problem of multi-modular reduction. The inverse problem is called multi-modular reconstruction. In what follows, we assume that have been fixed once and for all.

The simplest way to perform multi-modular reduction is to simply take

Inversely, Theorem 1 provides us with a formula for multi-modular reconstruction:

Since are fixed, the computation of the cofactors can be regarded as a pre-computation.

Let us analyze the cost of the above methods in terms of the complexity of -bit integer multiplication. Assume that for . Then multi-modular reduction can clearly be done in time .

As to multi-modular reconstruction, assume that for , where is such that . Cutting the cofactors in chunks of bits as in (4), we precompute the Fourier transform of all obtained chunks. The Fourier transforms of can be computed in time . The sum can be computed in the Fourier model in time and transformed back in time . Our assumption that for ensures that the computed sum is correct. The remainder can be computed in time . The total time of the reconstruction is therefore bounded by

If we only assume that , then we need to increase the bitsize by . If , then this means that we have to multiply the right-hand side of (8) by .

3.3.Scaled remainders

The above complexity analysis shows that naive multi-modular recomposition can be done faster than naive multi-modular reduction. In order to make the latter operation more efficient, one may work with scaled remainders that were introduced in [2]. The idea is that each remainder of the form is replaced by . The quotient is regarded as a real number and its remainder modulo as a number in the interval .

If we allow ourselves to compute with exact real numbers, then this leads us to replace the relation (6) by

and (7) by

For actual machine computations, we rather work with fixed point approximations of the scaled remainders. In order to retrieve the actual remainder from the scaled one , we need a circular -approximation of .

Now assume that with

Given a circular -approximation of in with

the algorithm at the end of section 2.3 allows us to compute a circular -approximation modulo of , by applying the formula (9). Since , we may then recover the number using one final -bit multiplication. Moreover, in the FFT-model, the transforms for need to be computed only once and the transforms for the numbers can be precomputed. In summary, given an approximation for the scaled remainder , we may thus compute approximations for the scaled remainders in time

From this, we may recover the actual remainders in time .

Scaled remainders can also be used for multi-modular reconstruction, but carry handling requires more care and the overall payoff is less when compared to the algorithm from the previous subsection.

3.4.Remainder trees

It is well-known that Chinese remaindering can be accelerated using a similar dichotomic technique as in the proof of Theorem 1. This time, we subdivide into parts with for . We denote and assume that (if , then we apply the native algorithms from the previous subsections).

Fast multi-modular reduction proceeds as follows. We first compute

using the algorithm from the previous subsection. Next, we recursively apply fast multi-modular reduction to obtain

The computation process can be represented in the form of a so called remainder tree; see Figure 3. The root of the tree is labeled by . The children of the root are the remainder trees for modulo , where . If needed, then the arity can be adjusted as a function of the bitsize of the moduli and .

Fast multi-modular reconstruction is done in a similar way, following the opposite direction. We first reconstruct

followed by

The computation flow can again be represented using a tree, but this time the computations are done in a bottom-up way.

Following the approach from subsection 3.3, it is also possible to systematically work with fixed point approximations of scaled remainders instead of usual remainders . In that case, the computation process gives rise to a scaled remainder tree as in Figure 4. Of course, the precision of the approximations has to be chosen with care. Before we come to this, let us first show how to choose .

 Figure 3. Example of a remainder tree with arities and at the two recursion levels. In the case of a reduction, the remainders are computed top-down. In the case of a reconstruction, they are reconstructed in a bottom-up fashion.
 Figure 4. The scaled remainder tree corresponding to Example 3, while using fixed point approximations for the scaled remainders.

3.5.Specification of the arities of nodes in the remainder tree

Let us first focus on multi-modular reconstruction. In order to make this process as efficient as possible, the arity should be selected as a function of and so as to make as large as possible in (8), while remaining negligible with respect to . Let

For inborn FFT schemes, we notice that

For the root of the remainder tree, we take

Using the same formula recursively for the other levels of the remainder tree, it is natural to define the sequences and by , and for ; the construction stops as soon as . Notice that we always have . The precise choice of is motivated by the following lemmas.

Lemma 2. If , then .

Proof. We clearly cannot have , since otherwise . If , then

If , then

This proves the result in all cases.

Lemma 3. If , then we have for .

Proof. For , we have , so that and whenever is sufficiently large. By construction, we also have

By induction, it follows that

for . We also have .

Lemma 4. We have .

Proof. Let and let be smallest such that . For all , we have , whence

Let be a constant such that for all . We also have

so that

Since for all , it follows that

Adding up (16) and (17) while letting tend to zero, the result follows.

Lemma 5. We have .

Proof. By construction, . For large , this means that , since otherwise , which is impossible. Consequently, .

3.6.Complexity analysis of multi-modular reconstruction

Let be the complexity multi-modular reconstruction for fixed moduli with for .

Theorem 6. If integer multiplication is based on an inborn FFT scheme, then

This bound holds uniformly in for .

Proof. In the special case when , the bound (8) yields

and we are done. If , then (8) implies

for . By induction, and using the fact that , we get

The result now follows from Lemma 4 and (1).

Remark 7. For and , the bound (18) simplifies into

If , then the bound becomes

Remark 8. It is interesting to examine the cost of the precomputations as a function of the parameters , and . For a node of the remainder tree at level , we essentially need to compute cofactors and their transforms. This can be done in time . Since we have nodes at level , the total precomputation time at level is therefore bounded by . Now . Consequently, the total precomputation time is bounded by

3.7.Complexity analysis of multi-modular reduction

Let us now consider the complexity of multi-modular reduction for fixed moduli with for . In this case, it is most efficient to work with scaled remainders, so the algorithm contains three main steps:

1. The initial conversion of into (an approximation of) .

2. The computation of (approximations of) the scaled remainders .

3. The final conversions of (approximations of) into .

At a first stage, we will assume that , where is sufficiently small such that the final approximations of the scaled remainders allow us to recover the usual remainders .

Let denote the cost of step 2. The conversions in steps 1 and 3 boil down to multiplications with fixed arguments, so that

For step 2, we use the scaled remainder tree algorithm from subsection 3.4, while taking the arities as in subsection 3.5.

Our next task is to show that is small enough in order to recover all remainders .

Lemma 9. There exists a constant such that for all and , we have

Proof. For the result clearly holds, so assume that . In particular, if is sufficiently large, then it follows that . Now assume for contradiction that . Then we would get . This is impossible for .

Now starting with a circular -approximation of , the scaled reduction algorithm from subsection 3.3 yields circular -approximations for the scaled remainders at level . Lemma 9 now shows that is sufficiently small for a continued application of the same algorithm for the next level. By induction over , the same reasoning shows that the scaled remainders at the -th level are computed with an error below

At the very end, we obtain -approximations for the scaled remainders . Since , this allows us to reconstruct each remainder using a multiplication by . This shows that is indeed small enough for the algorithm to work.

Theorem 10. If integer multiplication is based on an inborn FFT scheme, then

This bound holds uniformly in for .

Proof. A similar cost analysis as in the proof of Theorem 6 yields

when and

when . In both cases, combination with (19) yields

Notice also that , by Lemma 5.

Given a number , we may construct a similar sequence when using in the role of . Taking minimal such that , we have

Moreover, , which implies . Plugging this into (21), the result follows, since and the assumption that is increasing implies .

Remark 11. For and , the bound (20) simplifies into

If , then the bound (20) becomes

For very large with , this yields

Remark 12. Using a similar analysis as in Remark 8, the cost of all precomputations as a function of , and is again bounded by

3.8.Base conversion

The approach of this section can also be used for the related problem such as base conversion. Let and be a fixed base and order. Given a number , the problem is to compute with

Inversely, one may wish to reconstruct from . It is well known that both problems can be solved using a similar remainder tree process as in the case of Chinese remainders. The analogues for the formulas (7) and (9) are (22) and

The analogue of the recursive process of subsection 3.4 reduces a problem of size to similar problems of size and one similar problem of size but for the base . A routine verification shows that the complexity bounds (18) and (20) also apply in this context.

Moreover, for nodes of the remainder tree at the same level, the analogues of the cofactors and the multiplicands in (9) do not vary as a function of the node. For this reason, the required precomputations as a function of and can actually be done much faster. This makes it possible to drop the hypothesis that and are fixed and consider these parameters as part of the input. Let us denote by and the complexities of computing as a function of and vice versa.

Theorem 13. If integer multiplication is based on an inborn FFT scheme, then

These bound holds uniformly in for .

Proof. Let us estimate the cost of the precomputations as a function of and . The analysis is similar as in Remark 8 except that we only have to do the precomputations for a single node of the tree at each level . Consequently, the precomputation time is now bounded by . Since the decrease with at least geometric speed, this cost is dominated by . This proves the result under the condition that .

If , then we need to construct in a slightly different way. Assuming that is sufficiently large, let be maximal such that

Notice that

We again set and . This time, we take for and proceed with the usual construction for . It follows that

and

Using the new bound for , a similar complexity analysis as in the proofs of Theorems 6 and 10 yields the bounds (24) and (25) when forgetting about the cost of the precomputations. Now the cost of the precomputations for the first levels is bounded by

and the cost for the remaining levels by

We conclude that the right-hand sides of (24) and (25) absorb the cost of the precomputations.

3.9.Polynomial analogues

It is quite straightforward to adapt the theory of this section to univariate polynomials instead of integers. An example of this kind of adaptations can be found in [2]. In particular, this yields efficient algorithms for multi-point evaluation and interpolation in the case when the evaluation points are fixed. The analogues of our algorithms for base conversion yield efficient methods for -adic expansions and reconstruction.

More precisely, let be an effective commutative ring and let be the cost of multiplying two polynomials of degree in . Assume that allows for inborn FFT multiplication. Then , where and satisfy similar properties as in the integer case. Let be monic polynomials of degree . Given a polynomial of degree in we may then compute the remainders for in time

The reconstruction of from these remainders can be done in time

The assumption that admits a suitable inborn FFT scheme is in particular satisfied if is a finite field [14]. When working in an algebraic complexity model, this is still the case if is any field of positive characteristic [14]. For general fields of characteristic zero, the best known FFT schemes rely on the adjunction of artificial roots of unity [6]. In that case, our techniques only result in an asymptotic speed-up by a factor instead of . Nevertheless, the field of complex numbers does admit roots of unity of any order, and our algorithms can be used for fixed point approximations of complex numbers at any precision.

Multi-point evaluation has several interesting applications, but it is important to keep in mind that our speed-ups only apply when the moduli are fixed. For instance, assume that we computed approximate zeros to a complex polynomial of degree , using a bit-precision . Then we may use multi-point evaluation in order to apply Newton's method simultaneously to all roots and find approximations of bit-precision . Unfortunately, our speed-up does not work in this case, since the approximate zeros are not fixed. On the other hand, if the polynomial has degree instead of and we are still given approximate zeros (say all zeros in some disk), then the same moduli are used times, and one may hope for some speed-up when using our methods.

At another extremity, it is instructive to consider the approximate evaluation of a fixed polynomial with fixed point coefficients at a single point . We may thus assume that suitable Fourier transforms of the have been precomputed. Now we rewrite with , and . In order to evaluate each of the at , it suffices to transform the numbers , to perform the evaluations in the Fourier representation and then transform the results back. This can be achieved in time . We may finally compute using Horner's method, in time . For large , the dominant term of the computation time is .

4.Gentle moduli

4.1.The base algorithms revisited for special moduli

Let us now reconsider the naive algorithms from section 3.2, but in the case when the moduli are all close to a specific power of two. More precisely, we assume that

where and a small number. As usual, we assume that the are pairwise coprime and we let .

For such moduli, the naive algorithm for the euclidean division of a number by becomes particularly simple and essentially boils down to the multiplication of with the quotient of this division. In other words, the remainder can be computed in time instead of . For small values of , and , this gives rise to a speedup by a factor at least. More generally, the computation of remainders can be done in time .

Multi-modular reconstruction can also be done faster, as follows, using a similar technique as in [5]. Let . Besides the usual binary representation of and the multi-modular representation , it is also possible to use the Newton representation

where . Let us now show how to obtain efficiently from . Since , we must take . Assume that have been computed. For we next compute

using and

Notice that can be computed in time . We have

Now the inverse of modulo can be precomputed. We finally compute

which can be done in time . For small values of , we notice that it may be faster to divide successively by modulo instead of multiplying with . In total, the computation of the Newton representation can be done in time . Having computed the Newton representation, we next compute

for , using the recurrence relation

Since , the computation of takes a time . Altogether, the computation of from can therefore be done in time .

4.2.The gentle modulus hunt

For practical applications, we usually wish to work with moduli that fit into one word or half a word. Now the algorithm from the previous subsection is particularly efficient if the numbers also fit into one word or half a word. This means that we need to impose the further requirement that each modulus can be factored

with . If this is possible, then the are called -gentle moduli. For given bitsizes and , the main questions are now: do such moduli indeed exist? If so, then how to find them?

If , then it is easy to construct -gentle moduli by taking , where is odd. Indeed,

and . Unfortunately, this trick does not generalize to higher values . Indeed, consider a product

where are small compared to . If the coefficient of vanishes, then the coefficient of becomes the opposite of a sum of squares. In particular, both coefficients cannot vanish simultaneously, unless .

If , then we are left with the option to search -gentle moduli by brute force. As long as is “reasonably small” (say ), the probability to hit an -gentle modulus for a randomly chosen often remains significantly larger than . We may then use sieving to find such moduli. By what precedes, it is also desirable to systematically take for . This has the additional benefit that we “only” have to consider possibilities for .

4.3.The sieving procedure

We implemented a sieving procedure in Mathemagix [18] that uses the Mpari package with an interface to Pari-GP [22]. Given parameters and , the goal of our procedure is to find -gentle moduli of the form

with the constraints that

for , and . The parameter is small and even. One should interpret and as the intended and maximal bitsize of the small moduli . The parameter stands for the minimal bitsize of a prime factor of . The parameter should be such that fits into a machine word.

In Table 1 below we have shown some experimental results for this sieving procedure in the case when , , and . For , the table provides us with , the moduli , as well as the smallest prime power factors of the product . Many hits admit small prime factors, which increases the risk that different hits are not coprime. For instance, the number divides both and , whence these -gentle moduli cannot be selected simultaneously (except if one is ready to sacrifice a few bits by working modulo instead of ).

In the case when we use multi-modular arithmetic for computations with rational numbers instead of integers (see [10, section 5 and, more particularly, section 5.10]), then small prime factors should completely be prohibited, since they increase the probability of divisions by zero. For such applications, it is therefore desirable that are all prime. In our table, this occurs for (we indicated this by highlighting the list of prime factors of ).

In order to make multi-modular reduction and reconstruction as efficient as possible, a desirable property of the moduli is that they either divide or . In our table, we highlighted the for which this happens. We notice that this is automatically the case if are all prime. If only a small number of (say a single one) do not divide either or , then we remark that it should still be possible to design reasonably efficient ad hoc algorithms for multi-modular reduction and reconstruction.

Another desirable property of the moduli is that is as small as possible: the spare bits can for instance be used to speed up matrix multiplication modulo . Notice however that one “occasional” large modulus only impacts on one out of modular matrix products; this alleviates the negative impact of such moduli. We refer to section 4.5 below for more details.

For actual applications, one should select gentle moduli that combine all desirable properties mentioned above. If not enough such moduli can be found, then it it depends on the application which criteria are most important and which ones can be released.

 27657 28867 4365919 6343559 13248371 20526577 25042063 57267 416459 1278617 2041469 6879443 25754563 28268089 77565 7759 8077463 8261833 18751793 19509473 28741799 95253 724567 965411 3993107 4382527 19140643 23236813 294537 190297 283729 8804561 19522819 19861189 29537129 311385 145991 4440391 4888427 6812881 7796203 32346631 348597 114299 643619 6190673 11389121 32355397 32442427 376563 175897 1785527 2715133 7047419 30030061 30168739 462165 39841 3746641 7550339 13195943 18119681 20203643 559713 353201 873023 2595031 11217163 18624077 32569529 649485 21727 1186571 14199517 15248119 31033397 31430173 656997 233341 1523807 5654437 8563679 17566069 18001723 735753 115151 923207 3040187 23655187 26289379 27088541 801687 873767 1136111 3245041 7357871 8826871 26023391 826863 187177 943099 6839467 11439319 12923753 30502721 862143 15373 3115219 11890829 18563267 19622017 26248351 877623 514649 654749 4034687 4276583 27931549 33525223 892455 91453 2660297 3448999 12237457 21065299 25169783

Table 1. List of -gentle moduli for , , and .

4.4.Influence of the parameters , and

Ideally speaking, we want to be as large as possible. Furthermore, in order to waste as few bits as possible, should be close to the word size (or half of it) and should be minimized. When using double precision floating point arithmetic, this means that we wish to take . Whenever we have efficient hardware support for integer arithmetic, then we might prefer .

Let us start by studying the influence of on the number of hits. In Table 2, we have increased by one with respect to Table 1. This results in an approximate increase of the “capacity” of the modulus . On the one hand, we observe that the hit rate of the sieve procedure roughly decreases by a factor of thirty. On the other hand, we notice that the rare gentle moduli that we do find are often of high quality (on four occasions the moduli are all prime in Table 2).

 936465 543889 4920329 12408421 15115957 24645539 28167253 2475879 867689 4051001 11023091 13219163 24046943 28290833 3205689 110161 12290741 16762897 22976783 25740731 25958183 3932205 4244431 5180213 5474789 8058377 14140817 25402873 5665359 241739 5084221 18693097 21474613 23893447 29558531 5998191 30971 21307063 21919111 22953967 31415123 33407281 6762459 3905819 5996041 7513223 7911173 8584189 29160587 9245919 2749717 4002833 8274689 9800633 15046937 25943587 9655335 119809 9512309 20179259 21664469 22954369 30468101 12356475 1842887 2720359 7216357 13607779 23538769 30069449 15257781 1012619 5408467 9547273 11431841 20472121 28474807

Table 2. List of -gentle moduli for , , and .

Without surprise, the hit rate also sharply decreases if we attempt to increase . The results for and are shown in Table 3. A further infortunate side effect is that the quality of the gentle moduli that we do find also decreases. Indeed, on the one hand, tends to systematically admit at least one small prime factor. On the other hand, it is rarely the case that each divides either or (this might nevertheless be the case for other recombinations of the prime factors of , but only modulo a further increase of ).

 889305 50551 1146547 4312709 5888899 14533283 16044143 16257529 17164793 2447427 53407 689303 3666613 4837253 7944481 21607589 25976179 32897273 2674557 109841 1843447 2624971 5653049 7030883 8334373 18557837 29313433 3964365 10501 2464403 6335801 9625841 10329269 13186219 17436197 25553771 4237383 10859 3248809 5940709 6557599 9566959 11249039 22707323 28518509 5312763 517877 616529 879169 4689089 9034687 11849077 24539909 27699229 6785367 22013 1408219 4466089 7867589 9176941 12150997 26724877 29507689 7929033 30781 730859 4756351 9404807 13807231 15433939 19766077 22596193 8168565 10667 3133103 3245621 6663029 15270019 18957559 20791819 22018021 8186205 41047 2122039 2410867 6611533 9515951 14582849 16507739 30115277

Table 3. List of -gentle moduli for , , and .

An increase of while maintaining and fixed also results in a decrease of the hit rate. Nevertheless, when going from (floating point arithmetic) to (integer arithmetic), this is counterbalanced by the fact that can also be taken larger (namely ); see Table 4 for a concrete example. When doubling and while keeping the same upper bound for , the hit rate remains more or less unchanged, but the rate of high quality hits tends to decrease somewhat: see Table 5.

It should be possible to analyze the hit rate as a function of the parameters , , and from a probabilistic point of view, using the idea that a random number is prime with probability . However, from a practical point of view, the priority is to focus on the case when . For the most significant choices of parameters and , it seems in principle to be possible to compile full tables of -gentle moduli. Unfortunately, our current implementation is still somewhat inefficient for . A helpful feature for upcoming versions of Pari would be a function to find all prime factors of an integer below a specified maximum (the current version only does this for prime factors that can be tabulated).

 303513 42947057 53568313 331496959 382981453 1089261409 1176003149 851463 10195123 213437143 470595299 522887483 692654273 1008798563 1001373 307261 611187931 936166801 1137875633 1196117147 1563634747 1422507 3950603 349507391 490215667 684876553 693342113 1164052193 1446963 7068563 94667021 313871791 877885639 1009764377 2009551553 1551267 303551 383417351 610444753 1178193077 2101890797 2126487631 1555365 16360997 65165071 369550981 507979403 1067200639 1751653069 4003545 20601941 144707873 203956547 624375041 655374931 1503716491 4325475 11677753 139113383 210843443 659463289 936654347 1768402001 4702665 8221903 131321017 296701997 496437899 1485084431 1584149417 5231445 25265791 49122743 433700843 474825677 907918279 1612324823 5425527 37197571 145692101 250849363 291039937 456174539 2072965393 6883797 97798097 124868683 180349291 234776683 842430863 858917923 7989543 4833137 50181011 604045619 638131951 1986024421 2015143349

Table 4. List of -gentle moduli for , , and . Followed by some of the next gentle moduli for which each divides either or .

 15123 380344780931 774267432193 463904018985637 591951338196847 34023 9053503517 13181369695139 680835893479031 723236090375863 40617 3500059133 510738813367 824394263006533 1039946916817703 87363 745270007 55797244348441 224580313861483 886387548974947 95007 40134716987 2565724842229 130760921456911 393701833767607 101307 72633113401 12070694419543 95036720090209 183377870340761 140313 13370367761 202513228811 397041457462499 897476961701171 193533 35210831 15416115621749 727365428298107 770048329509499 519747 34123521053 685883716741 705516472454581 836861326275781 637863 554285276371 1345202287357 344203886091451 463103013579761 775173 322131291353 379775454593 194236314135719 1026557288284007 913113 704777248393 1413212491811 217740328855369 261977228819083 1400583 21426322331 42328735049 411780268096919 626448556280293

Table 5. List of -gentle moduli for , , and . Followed by some of the next gentle moduli for which each divides either or .

4.5.Application to matrix multiplication

One of our favourite applications of multi-modular arithmetic is the multiplication of integer matrices . We proceed as follows:

1. Compute and for , using multi-modular reductions.

2. Multiply for .

3. Reconstruct using multi-modular reconstructions.

If is larger than for all and , then can be read off from .

From a practical point of view, the second step can be implemented very efficiently if fits into the size of a word. When using floating point arithmetic, this means that we should have for all . For large values of , this is unrealistic; in that case, we subdivide the matrices into smaller matrices with . The fact that may depend on is very significant. First of all, the larger we can take , the faster we can multiply matrices modulo . Secondly, the in the tables from the previous sections often vary in bitsize. It frequently happens that we may take all large except for the last modulus . The fact that matrix multiplications modulo the worst modulus are somewhat slower is compensated by the fact that they only account for one out of every modular matrix products.

Several of the tables in the previous subsections were made with the application to integer matrix multiplication in mind. Consider for instance the modulus from Table 1. When using floating point arithmetic, we obtain , , , , and . Clearly, there is a trade-off between the efficiency of the modular matrix multiplications (high values of are better) and the bitsize of (larger capacities are better).

If is large with respect to , then the modular matrix multiplication step is the main bottleneck, so it is important to take all approximately of the same size (i.e. should be small) and in such a way that the corresponding lead to the best complexity ( tends to work well). This can often only be achieved by lowering to or . For closer to , the Chinese remaindering steps become more and more expensive, which makes it interesting to consider larger values of and to increase the difference . For significantly below , we resort to FFT-based matrix multiplication. This corresponds to taking roughly twice as many moduli, but the transformations become approximately times less expensive.

4.6.Implementation issues

We have not yet implemented any of the algorithms in this paper. The implementation that we envision selects appropriate code as a function of the parameters , and . The parameters and highly depend on the application: see the above discussion in the case of integer matrix multiplication. Generally speaking, is bounded by the bitsize of a machine word or by .

So far, we have described four main strategies for solving Chinese remaindering problems for moduli :

G

The “gentle modulus strategy” requires to be an -gentle modulus with . Conversions between and can be done fast. The further conversions from and to are done in a naive manner.

P

The “gentle product strategy” strategy requires to be gentle moduli. We now perform the multi-modular reductions and recombinations using the algorithms from subsection 4.1.

N

The “naive strategy” for Chinese remaindering has been described in subsections 3.2 and 3.3.

F

The asymptotically fast “FFT-based strategy” from Theorems 6 and 10, which attempts to do as much work as possible using Fourier representations.

The idea is now to apply each of these strategies at the appropriate levels of the remainder tree. At the very bottom, we use G, followed by P and N. At the top levels, we use F.

As a function of , we need to decide how much work we wish to perform at each of these levels. For small , it suffices to combine the strategies G and P. As soon as exceeds , some of the modular reductions in the strategy G may become expensive. For this reason, it is generally better to let P do a bit more work, i.e. to take . It is also a good practice to use as the “soft word size” for multiple integer arithmetic (see subsection 2.4).

As soon as starts to get somewhat larger, say , then some intermediate levels may be necessary before that FFT multiplication becomes plainly efficient. For these levels, we use strategy N with arity two. It still has to be found out how large this intermediate region is, exactly. Indeed, the ability to do more work using the Fourier representation often lowers the threshold at which such methods become efficient. If one manages to design extremely efficient implementations for the strategies G and P (e.g. if multi-modular reduction and reconstruction can approximately be done as fast as multiplication itself), then one may also consider the use of Chinese remaindering instead of Fourier transforms in F.

For large , the FFT-based algorithms should become fastest and we expect that the theoretical speed-up should result in practical gains of a factor three at least. As emphasized before, it is crucial to rely on inborn FFT strategies. For very large , we may also run out of -gentle moduli. In that case, we may need to resort to lower values of , with the consequence that some of the lower levels may become somewhat more expensive.

4.7.Alternative moduli

A very intriguing question is whether it is possible to select moduli that allow for even faster Chinese remaindering. In the analogue case of polynomials, highly efficient algorithms exist for multi-point evaluation and interpolation at points that form a geometric sequence [5]. Geometric sequences of moduli do not work in our case, since they violate the mutual coprime requirement and they grow too fast anyway. Arithmetic progressions are more promising, although only algorithms with a better constant factor are known in the polynomial case [5], and the elements of such sequences generically admit small prime factors in common that have to be treated with care.

Another natural idea is to chose products of the form , where is highly composite. Each divisor of naturally induces a divisor of , where denotes the -th cyclotomic polynomial. For instance, the number is divisible by , , and . Now euclidean division by and is easy since these numbers admit an extremely sparse and regular binary representations. Furthermore, the numbers and can both be factored into products of two integers of less than bits. Taking and , it should therefore be possible to design a reasonably efficient Chinese remaindering scheme for on a 64-bit architecture.

There are several downsides to this approach. Most importantly, the largest prime divisor of grows quite quickly with , even if is very smooth. For instance, the four largest for which are and the four largest for which are . By construction, the number also admits many divisors if is smooth, so several moduli of this type are generally not coprime. Both obstructions together imply that we can only construct Chinese remaindering schemes with a small number of moduli in this way. For instance, on a 64 bit architecture, one of the best schemes would use and twelve moduli of an average size of bits. In comparison, there are many -gentle and -gentle moduli that are products of prime numbers of approximately bits (or more), and combining a few of them leads to efficient Chinese remaindering schemes for twelve (or more) prime moduli of approximately bits.

Bibliography

[1] D. Bernstein. Removing redundancy in high precision Newton iteration. Available from http://cr.yp.to/fastnewton.html, 2000.

[2] D. Bernstein. Scaled remainder trees. Available from https://cr.yp.to/arith/scaledmod-20040820.pdf, 2004.

[3] A. Borodin and R.T. Moenck. Fast modular transforms. Journal of Computer and System Sciences, 8:366–386, 1974.

[4] A. Bostan, G. Lecerf, and É. Schost. Tellegen's principle into practice. In Proceedings of ISSAC 2003, pages 37–44. ACM Press, 2003.

[5] A. Bostan and É. Schost. Polynomial evaluation and interpolation on special sets of points. Journal of Complexity, 21(4):420–446, August 2005. Festschrift for the 70th Birthday of Arnold Schönhage.

[6] D.G. Cantor and E. Kaltofen. On fast multiplication of polynomials over arbitrary algebras. Acta Informatica, 28:693–701, 1991.

[7] J.W. Cooley and J.W. Tukey. An algorithm for the machine calculation of complex Fourier series. Math. Computat., 19:297–301, 1965.

[8] C.M. Fiduccia. Polynomial evaluation via the division algorithm: the fast fourier transform revisited. In A.L. Rosenberg, editor, Fourth annual ACM symposium on theory of computing, pages 88–93, 1972.

[9] M. Fürer. Faster integer multiplication. SIAM J. Comp., 39(3):979–1005, 2009.

[10] J. von zur Gathen and J. Gerhard. Modern Computer Algebra. Cambridge University Press, New York, NY, USA, 3rd edition, 2013.

[11] T. Granlund et al. GMP, the GNU multiple precision arithmetic library. http://www.swox.com/gmp, 1991.

[12] G. Hanrot, M. Quercia, and P. Zimmermann. The middle product algorithm I. speeding up the division and square root of power series. Accepted for publication in AAECC, 2002.

[13] D. Harvey and J. van der Hoeven. On the complexity of integer matrix multiplication. Technical report, HAL, 2014. http://hal.archives-ouvertes.fr/hal-01071191, accepted for publication in JSC.

[14] D. Harvey, J. van der Hoeven, and G. Lecerf. Faster polynomial multiplication over finite fields. Technical report, ArXiv, 2014. http://arxiv.org/abs/1407.3361, accepted for publication in JACM.

[15] D. Harvey, J. van der Hoeven, and G. Lecerf. Even faster integer multiplication. Journal of Complexity, 36:1–30, 2016.

[16] J. van der Hoeven. Newton's method and FFT trading. JSC, 45(8):857–878, 2010.

[17] J. van der Hoeven and G. Lecerf. Faster FFTs in medium precision. In 22nd Symposium on Computer Arithmetic (ARITH), pages 75–82, June 2015.

[18] J. van der Hoeven, G. Lecerf, B. Mourrain, et al. Mathemagix, 2002. http://www.mathemagix.org.

[19] J. van der Hoeven, G. Lecerf, and G. Quintin. Modular SIMD arithmetic in Mathemagix. ACM Trans. Math. Softw., 43(1):5:1–5:37, 2016.

[20] R.T. Moenck and A. Borodin. Fast modular transforms via division. In Thirteenth annual IEEE symposium on switching and automata theory, pages 90–96, Univ. Maryland, College Park, Md., 1972.