Jean Claude Bajard
email: jean.bajard@inria.fr

Kazuhide Fukushima
email: ka-fukushima@kddi-research.jp

Thomas Plantard
email: thomas.plantard@gmail.com

Arnaud Sipasseuth
email: ar-sipasseuth@kddi-research.jp

Generating Very Large RNS Bases

Keywords: Residue Number Systems, Setwise Coprime, Modular Arithmetic, Cryptography. !

Residue Number Systems (RNS) are proven to be effective in speeding up computations involving additions and products. For these representations, there exists efficient modular reduction algorithms that can be used in the context of arithmetic over finite fields or modulo large numbers, especially when used in the context of cryptographic engineering. Their independence allows random draws of bases, which also makes it possible to protect against side-channel attacks, or even to detect them using redundancy. These systems are easily scalable, however the existence of large bases for some specific uses remains a difficult question. In this paper, we present four techniques to extract RNS bases from specific sets of integers, giving better performance and flexibility to previous works in the litterature. While our techniques do not allow to solve efficiently every possible case, we provide techniques to provably and efficiently find the largest possible available RNS bases in several cases, improving the state-of-the-art on various works of the recent literature.

INTRODUCTION

Scientifc Context: Unconventional arithmetic offers original and effective approaches in various domains of application [START_REF] Sousa | Nonconventional computer arithmetic circuits, systems and applications[END_REF]. The use of Residue Number System (RNS) appeared in the 1950s as arithmetic [START_REF] Svoboda | Operational circuits[END_REF], [START_REF] Garner | The residue number system[END_REF]. This representation is directly inspired by the Chinese remainder theorem [START_REF] Szabo | Residue Arithmetic and its Applications to Computer Technology[END_REF], [START_REF] Knuth | Seminumerical Algorithms[END_REF]. The main interest lies in the speed of the addition and multiplication operations which are distributed on small values that are the modular remainders on a set of pairwise coprime numbers named RNS base. Despite a conversion cost which in the worst case is quadratic in the size of the base, calculations involving additions and products become extremely profitable (i.e. the inner product, convolution products, operations on large numbers, etc.). However, division and comparison remain costly operations in RNS, recent works propose interesting approaches [START_REF] Kawamura | Efficient algorithms for sign detection in rns using approximate reciprocals[END_REF], [START_REF] Isupov | Interval estimation of relative values in residue number system[END_REF], [START_REF] Babenko | Improved modular division implementation with the akushsky core function[END_REF]. The three main areas of application are signal processing [START_REF] Jullien | Application of the residue number system to computer processing of digital signals[END_REF], [START_REF] Chaves | 2 n + 1, 2 n+k , 2 n -1: A new RNS moduli set extension[END_REF], [START_REF] Hiasat | On the design of RNS intermodulo processing units for the arithmetic-friendly moduli sets {2 n+k , 2 n -1, 2 n+1 -1}[END_REF], cryptography but also in theoretical computer science to reach complexity bounds [START_REF] Beame | Log depth circuits for division and related problems[END_REF], [START_REF] Van Der Hoeven | Fast chinese remaindering in practice[END_REF]. The work in this paper is relevant for all applications on large numbers including cryptography since the 90's [START_REF] Schoinianakis | Residue arithmetic systems in cryptography: a survey on modern security applications[END_REF] with RSA, DH, ECC [START_REF] Bigou | Single base modular multiplication for efficient hardware RNS implementations of ECC[END_REF], [START_REF] Asif | Highly parallel modular multiplier for elliptic curve cryptography in residue number system[END_REF], pairing [START_REF] Duquesne | Rns arithmetic in F p k and application to fast pairing computation[END_REF], Euclidean lattices, homomorphic protocols [START_REF] Halevi | An improved rns variant of the bfv homomorphic encryption scheme[END_REF], [START_REF] Kim | Revisiting homomorphic encryption schemes for finite fields[END_REF], etc.

In a domain such as cryptography, arithmetic operations are performed modulo large numbers that are often prime, the use of RNS becomes more complicate as modular reduction requires a conversion of RNS bases [START_REF] Posch | Modulo reduction in residue number systems[END_REF], [START_REF] Kawamura | Cox-rower architecture for fast parallel montgomery multiplication[END_REF], [START_REF] Djath | Hierarchical approach in rns base extension for asymmetric cryptography[END_REF], [START_REF] Ahsan | Efficient fpga implementation of rns montgomery multiplication using balanced rns bases[END_REF]. This last point has generated a rich literature, in particular around the choice of bases for efficient implementations [START_REF] Bigou | Single base modular multiplication for efficient hardware RNS implementations of ECC[END_REF], [START_REF] Cheung | FPGA implementation of pairings using residue number system and lazy reduction[END_REF], [START_REF] Lesavourey | Efficient leak resistant modular exponentiation in RNS[END_REF], [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF]. RNS are also particularly interesting for countering attacks by faults, as the addition of redundancy elements at the base level makes it possible to set up fault detection [START_REF] Bajard | Multi-fault attack detection for RNS cryptographic architecture[END_REF]. Finally, the random drawing of bases ensures that the same calculation produces different patterns at each evaluation, making learning possible leakage of information more difficult [START_REF] Courtois | Resilience of randomized rns arithmetic with respect to side-channel leaks of cryptographic computation[END_REF].

Implementations are constrained by the limitations of the architectures we deals with. The question we are faced with is what sizes can be reached. This is very dependent on the operators of the target architecture which limit the size of the elements of the RNS base. Finding the largest possible base with this constraint on its elements thus becomes a major issue [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF], [START_REF] Bigou | Hybrid position-residues number system[END_REF]. In order to answer this challenge, an approach proposes after a filtering to build a co-primality graph and to perform a brute force search of maximal cliques [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]. This problem is known to be NP-complete [START_REF] Karp | Reducibility among combinatorial problems[END_REF], [START_REF] Tomita | Efficient algorithms for finding maximum and maximal cliques and their applications[END_REF].

In order to achieve larger RNS base sizes, considering the cost of finding a maximal clique, filtering to reduce the graph size becomes crucial.

Main results: In this paper, we introduce four filtering methods that reduce the size of the graph. Our approach significantly improves the one of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]. The RNS bases we are able to build are larger than those proposed in the literature by maintaining the same constraints. One can find in [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] a large overview of RNS bases with their advantages. The first algorithm is targeted towards intervals. The second algorithm is targeted towards random sets. The third algorithm is an adaptation of the first algorithm for "non-interval semirandom" sets. The fourth one is a heuristic approach, aimed to accelerate computations.

In the literature: Different approaches found in the literature deal with obtaining an easy reduction for each modulo and an efficient RNS base extension. Most of them use classical pseudo-Mersenne or Solinas moduli. We mention here some particular approaches that condition the moduli. In the Cox-Rover algorithm [START_REF] Kawamura | Cox-rower architecture for fast parallel montgomery multiplication[END_REF] authors suggest to use pseudo-Mersenne coprimes with an interdependence between the number of moduli, their pseudo-Mersenne form and the truncation error of the computation. This could reduce significantly the possibilities depending on the target architecture. With the double Montgomery approach [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF], the authors use a Montgomery reduction at the modulo level, they relax the constraints which allows to significantly increase the number of moduli. Then quadratic RNS [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] offers an efficient bases exchange but adding a quadratic property of the moduli. Finally, a recent approach suggest friendly Montgomery numbers [START_REF] Bajard | Montgomery-friendly primes and applications to cryptography[END_REF] to increase the number of coprimes while keeping an efficient base extension. We will illustrate how our methods allow to build large RNS bases for each case.

Organization of the paper: After introducing some useful background for the understanding of this article, we introduce in Section 3 four filtering approaches. We present two basic alternatives to [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] for intervals and random sets, then specialize those basic alternatives to propose another two alternatives to deal with specialized sets for a total of four filtering techniques. Then, we propose new bases for the most relevant cases of families of moduli used in the state of the art of RNS optimization [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF].

BACKGROUND

Residue Number Systems and notations

Definition 1 (RNS). A Residue Number

System is defined by a set of pairwise coprime integers M = {m 1 , .., m d } ∈ (N * ) d named RNS base of size d. We denote M = m i the product size. For x ∈ Z, we denote x M = x 1 , ..., x d the RNS representation of x, where x i = x mod m i and call the set of x i the residues of x in M.

The Chinese remainder theorem (CRT) ensures that there is an isomorphism between d i=1 Z mi to Z M . Thus, for each x ∈ Z, x M is unique and corresponds to x M = x mod M . Sometimes the same notation is used for both M and M . Note that the search for a large RNS base, i.e. reaching large sizes for M , often means increasing d, as the elements of an RNS base are limited in size depending on the architectures used. If the machine words are of a given size n in binary, then each element of the database is smaller than 2 n . In practice, we would like to have the same size for all residues, in order to balance the load between all computing units.

We start by introducing some specific notations.

Notation 1. We note P the set of prime numbers. For any set I, we note I min its smallest element, I max its largest element. We note I P the set of primes that can divide at least one nonzero element of I (e.g. to simplify in some algorithms, I P can be replaced by the set of primes smaller than I max ).

Notation 2 (Sets of maximum possible pairwise coprimes).

We denote Φ(S, I) the set of all sets of pairwise coprime numbers of maximal size containing elements of S ∪ I. I is the initial set of numbers considered, and S is a set of filtered values. Almost all the values of S are pairwise coprimes.

Example 1. Let S = {} and I = {2, 3, 4, 11, 17, 121}:

• M 1 = {121, 17, 2, 3} ∈ Φ(S, I) • M 2 = {11, 17, 4, 3, 2} / ∈ Φ(S, I) as gcd(2, 4) = 1 • M 3 = {11, 17, 3} / ∈ Φ(S, I): |M 3 | < |M 1 | • M 4 = {11, 17, 2, 3} ∈ Φ(S, I) • M 5 = {11, 17, 4, 3} ∈ Φ(S, I) • M 6 = {121, 17, 4, 3} ∈ Φ(S, I)
By bruteforce, we can easily see that {M 1 , M 4 , M 5 , M 6 } = Φ(S, I).

The goal of this paper is to find at least one element of Φ({}, I) for specific I. Thus, we try to reduce the problem by finding smaller sets S ⊆ I and I ⊆ I such that Φ(S , I ) ⊆ Φ({}, I). In the previous example, If S = {121} and I = {2, 3, 4, 17} then Φ(S , I ) = {M 1 , M 6 } ⊆ Φ({}, I). At the end of the process of filtering we could obtain S = {2, 3, 4, 17, 121} and I = {}. The work of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] proposed a solution when I is an interval.

Maximum RNS bases are maximum cliques

In [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF], methods were proposed to find a maximum size RNS base within an interval of consecutive numbers. The authors of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] saw the problem of finding an element of Φ({}, I) as a graph problem. The approach is the following: 1) ∀a ∈ I, create a node tagged with a.

2) ∀a 1 , a 2 ∈ I s.t gcd(a 1 , a 2 ) = 1, add an edge (a 1 , a 2 ).

3) A maximum size RNS base from I is a clique from the resulting graph.

Example 2. Let I = [2, 13] be a set of 12 consecutive numbers. On the left is the graph constructed as per the above described process, and on the right is a maximum clique of the left graph, i.e., a subgraph of the largest possible size where every node is connected: it is also a graph that can be constructed from an element of Φ({}, I). Note that the solution is not unique: for example, node 2 can be replaced by node 4 or 8 and 3 by 9, thus Φ({}, I) has at least six elements. Φ({2, 3}, I) however, is a singleton. A clique, maximal or not, will always represent a set of pairwise coprime numbers in this model.

As also pointed by [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF], one issue with this initial approach is that the maximum clique problem is known to be NP-complete, thus before tackling the graph problem, [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] proposes to simplify the search space by "filtering" Φ(∅, I) into a smaller set. It is very clear from the above examples, that given a set I, we have Φ(I ∩P, I) ⊆ Φ(∅, I), thus part of the problem in finding an element of Φ(∅, I) can be reduced to find an element of Φ(I ∩ P, I), the difficulty comes from filtering Φ(I ∩ P, I) even further. Thus, they propose some solutions to ease the computations. Their main observation is the following lemma: In other words, [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] can construct a subset where there is at least one maximum RNS base containing this set S, reducing the search space. After constructing S, the search space can be even further simplified: one only needs to replace I with the set I S = {x ∈ I|∀m ∈ S, gcd(m, x) = 1} i.e the elements of I coprime to S. we thus have Φ(S, I S ) = Φ(S, I) ⊆ Φ({}, I), and finding an element of Φ(S, I S ) is computationally much easier.

As pointed by [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] the set S is easily computable. In fact, for some small parameters, [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] reduce I S to an empty set, outputting S ∈ Φ({}, I) as a solution. The approach however can have some limitations after computing S:

• For large I, I S can be too large to be reasonably processed and a clique algorithm might still be required.

• For large I, I S can be itself expensive to compute.

• If I is not an interval, the approach cannot be used.

To alleviate the first issue, [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] propose to heuristically prune the graph: constructing S was essentially pruning the graph from elements containing only one or two distinct primes, thus their suggestion was to continue further by pruning combinations of 3 distinct prime factors, then 4, etc... However, as they also pointed out, this heuristic pruning method beyond S does no longer guarantee a maximal size RNS base.

In our work, we propose upgraded pruning methods to help alleviate all three issues. Those methods, while not proven to be able to avoid the resolution of an NPcomplete problem, are sure to provide an optimal solution. We also provide a highly efficient greedy version of our methods which do not guarantee optionality, even if often achieved. They do simplify the computations by essentially giving a graph of a much lower size as an input to the maximum clique algorithm. Interestingly, in many of our applications, our methods transform our original input into a graph already corresponding to a clique, thus avoiding the necessity to call a maximum clique algorithm in the first place.

TWO NOVEL BASIC FILTERING TECHNIQUES

Here we propose two alternatives filtering techniques, initially proposed at [START_REF] Bajard | Generating residue number system bases[END_REF], that directly improve over the theoretical work of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]. Similarly to [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF], the aim is to minimize the size of the set we give as an input to a graph algorithm, while relying on a same principle of picking a provably optimal choice at each iteration to guarantee the maximality of the output.

The first method modifies [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] in a way such that computing gcd is never necessary, and that the initial set I is never stored, only prime numbers are (i.e. I P ). This number of primes can be significantly smaller than the size of the initial set I. The second method is not a direct improvement of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] but rather an alternative method to deal with sets that were not intervals: Lemma 1 can only applies if the initial set given as an entry is an interval. It is also more efficient than the other methods for sets of small size, even intervals, as it does not require any knowledge of the prime numbers.

Important note: These filtering techniques give at least one RNS base with the maximum number of moduli extracted from an initial set, they do not select the best base according to particular criteria (for example the size of the product of the moduli), among the, potentialy numerous, bases with the maximum number of moduli.

First method: filtering by factors recomposition

Constructing coprime moduli sets within an interval by analyzing their factorization was the core idea of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]. Their method call a clique algorithm over all remaining possible numbers, constructing S and I S beforehand. Our first improvement over [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] is the following: when dealing with an interval, there is no need to store the interval at all. We just need the list of any prime number that can be found in a decomposition of an integer within the interval. Prime lists and code to compute them can be found in various places on the internet ( [START_REF] Silva | Tables of values of pi(x) and of pi2(x)[END_REF], [START_REF] Caldwell | The first fifty million primes[END_REF], [37]) if primes required are not overly large (i.e beyond 2 64 ), and checking whether an integer is within an interval can be done with the help of only two values: I min and I max .

Here, the first method we present is essentially the following: we thin the leftover primes, thus lowering the possible combinations, accelerating the process. This algorithm can ensure a maximum RNS base after a call to a clique algorithm to leftover conflicting combinations, similarly to [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]. In [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF], their core observation is that any integer prime p > δ (with δ = I max -I min ) can only appear once in the decomposition of integers of I. We add a generalized observation, which severely limits the amount of prime numbers we have to consider.

Knowing that any prime p > √ I max cannot be powered in a said decomposition, we can further limit the amount of primes to consider. Let us denote β = max(δ, I max ). We explain below our method to compute E 1 the set of prime powers within I with the algorithm FirstStep, keeping information about the primes not being used for E 1 as those are necessary later.

Algorithm 1 FirstStep

Input: P all the primes smaller than β Here we reduced the problem to Φ(E 1 , I ) ⊆ Φ({}, I), where I is I minus all the primes or power of primes. Note that we do not store I nor I , this is one of the main advantage compared to the approach of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF].

E 1 ← P ∩ I for p ∈ P do if there is a k ≥ 1, such that p k ∈ I then Put the largest p k into E 1 Eliminate p from P Output: E 1 , P . Example 3. If I = [2 16 -2 8 , 2 16 ] then 1 P = {2, 3, 5, ..., 251} and E 1 = {65287, ..., 65521} 2
After FirstStep, if E 1 is not a solution by itself (most often not), then we can continue our process using the leftover primes within P to find candidates to expand E 1 .

In a SecondStep, we then construct E 2 the set of products of two distinct primes taking one prime factor as large as possible, and fill in a set RES product of two integers that lie outside [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]'s proof (lemma 1), but we wish to keep information on the primes we pick for E 2 and those we will not pick as those represent the necessary information to complete the missing elements of a RNS base of maximal size. When a product of two primes does not satisfy lemma 1, if this product belongs to I then we store it in a set RES, and if one of this two primes could appear in a product of three primes, which belongs to I, we store them in a set R. RES and R are used in the ThirdStep.

Construction of E 2 , R, RES is done with Algorithm Sec-ondStep. Note that in the case δ = √ I max , RES is empty at the end of SecondStep. We have now simplified the problem to Φ(E 1 ∪ E 2 , I) ⊆ Φ({}, I), and the only elements we need to check are now combinations of elements of R and elements of RES. ,I). However, we still need to further process the remaining candidates in the case RES and/or R are not empty. We basically reconstruct I S by bruteforce for S = E 1 ∪ E 2 in order to reduce the search to Φ(E 1 ∪ E 2 , I) = Φ(S, I S ) where I S is essentially composed of the integers within I coprime to S which should be exactly RES and the missing combinations involving primes in R. Constructing I S , if not proven to be empty, is done by algorithm ThirdStep, which stores in RES any possible product combination of the powers of the remaining primes within R fitting I.

After SecondStep we have constructed E 2 such that Φ(E 1 ∪ E 2 , I) ⊆ Φ(E 1 , I) ⊆ Φ({}

Algorithm 3 ThirdStep Input

: I min , I max , k, R,RES if k ≤ I max then if k ≥ I min then RES = RES ∪ {k} else if R = ∅ then p ← min(R) R ← R\{p} while k ≤ I max do RES ←ThirdStep(I min , I max , k, R,RES) k ← k * p Output: RES
ThirdStep can have a very large complexity for large R. However, the recursion falls short when the elements of R are large: the level of recursion excluding calls to K = 1 is bounded by log r (I max ) for r = min(R). Thus it is very important that previous steps eliminate as many small primes as possible before ThirdStep.

Note that the output of ThirdStep does not theoretically guarantee a set of pairwise coprime integers since given primes p, q, x, w, integers a = p × q × z and b = p × x × w could be valid integers within I. While this scenario has never occurred in our tests, we do not have any proof showing its impossibility. Thus, in theory, a maximum clique algorithm could still be needed. The entire process is thus the concatenation of the three above algorithms, described in Algorithm 4 FactorFilter. Do note that at the end of Algorithm SecondStep, if both R and RES are empty, then C is a maximum clique and ThirdStep would actually do nothing as there is nothing to search within, as E 1 and E 2 have been constructed using only "proven picks" that do not reduce the size of an optimal solution.

Algorithm 4 FactorFilter

Input: I min , I max the interval range. Output: C a set of pairwise coprime numbers, RES potential additions to grow C

1: E 1 , P ← FirstStep(I min , I max ) 2: E 2 , R, RES ← SecondStep(I min , I max , P ) 3: C ← E 1 ∪ E 2 4: RES← ThirdStep(I min , I max , R, 1, RES) 5: return C, RES If RES is empty, then C is maximal
During the whole process, not a single gcd has been computed, nor was any value of I beyond its extremities needed. The remainder set I S can also be emptied early i.e not constructed at all. Also note that the function Previous-Prime, while expensive in theory for arbitrary numbers in arbitrary applications, can be very efficiently called in our case if we already have access to the list of all necessary primes: we will never need a prime larger than I min beyond the construction of E 1 .

Second method: generalizing filtering for any set

This second proposed method works for any set, and is better than the first method on intervals anytime the set of integers I we need to manage is lower than the set of primes I P we have to manage. To the best of our knowledge, there is no other comparable preexisting technique besides the direct clique method.

This second method relies on the following lemma:

Lemma 2 (Divisor pick for general sets).

If ∃y = ab ∈ I and {gcd(y, x), x ∈ S} = {1}, s.t ∃a > 1, ∀x ∈ I, a| gcd(x, y) or gcd(x, y) = 1.

Then Φ(S ∪ {y}, I) ⊆ Φ(S, I)

This lemma can be understood by seeing that if such y exists, then for any element M ∈ Φ(S, A), either there is an element divisible by a which can then be replaced by y without changing pairwise coprimality relations, either there is not and thus M is not maximal as you could include y (a contradiction). The method to find an element of Φ(S, I) is then the following: The full process is described by Algorithm 5. Note that this filter can also be applied after the previous filter since the output of ThirdStep is also a set, thus our two algorithms are not mutually exclusive. After this filtering algorithm over I, to further increase the set size of S, we would need to call a maximum clique algorithm over the leftover I S : however, in our experiments, I S is often empty, leaving S as the final base. Unlike algorithm f ← gcd(d, f )

10: if f = 1 then Break #divisors ≥ 2 11: if f = 1 then #divisors < 2 12:
S ← S ∪ {m} 13:

I S ← {m ∈ I S |gcd(m, m ) = 1} 14:
Break 15: return S, I S

FITERING NON-INTERVAL SPECIFIC SETS

We propose here another two filtering methods to complement those presented above, in case the target set is neither an interval, nor is it a completely random set. First, one called in this paper the third method, which introduces a modification of the first method. This third method adapted to work on any set of integers, which can replace the above two methods for very large sets (but does not outperform the previous methods otherwise). We will show that this method also guarantees maximality of the result. The second one, called fourth method, is a simpler, but distinct, method for extracting a RNS base from any random set. This method does not guarantee maximality of the result and therefore cannot replace any of the algorithms we propose nor the original work of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]. However, the interesting part of this technique is its efficiency and simplicity and the possibility to prove ad-hoc maximality under certain conditions.

General Idea

Those new algorithms were designed to address some of the shortcomings of the previous two algorithms: one is effective for any given small set, the other starts to be effective for any given large interval but when the bounds have a similar binary size. Both algorithms do not give satisfactory results when the bounds of the interval have a different size, or when the set is a large combination of disjoint intervals. Note, that solving the problem efficiently under any set is difficult: any "maximum clique problem" can be transformed into an equivalent instance of a "largest pairwise coprime subset problem" in polynomial time, thus an efficient algorithm for any set is an efficient algorithm for a general NP-hard problem. However, we present here adapted algorithms when the entry sets are neither intervals nor completely random. Thus those supplementary algorithms are designed following two observations:

•

The efficiency of the first method is limited when the small factors have many possibilities, failing to pick a choice and resorts to exhaustive search in the final steps of the algorithm.

•

The efficiency of the second algorithm is limited by the fact that processing the factorization of each element of the entry set has to be done, thus limiting the size of the sets that can be reasonably be processed.

Thus we alleviate those issues with the following:

• Instead of just iterating on small factors, it might be better to also start with the largest factors, which are guaranteed to have very limited choices no matter how wide are the bounds. Choosing large factors first will also eliminate some small factors that could have had many solutions.

•

If we have a low-cost function (in memory and speed) which determines whether an integer is part of the initial set or not, we could recompose by factors and use that function without ever storing the set. Whenever this is possible, we can achieve an efficiency gain on any large structured set over the generic method.

Third method: generalizing the recomposition technique

The third method reuses the philosophy of the intervalbased algorithm. However as some properties on intervals might no longer hold on generic sets we have to make some generalizations.

Properties to justify maximality

As we are managing more diverse sets, we also need more general properties to justify our algorithms will still output a maximum size RNS.

Property 1 (Prime powers can always included).

Let a ∈ P, ∀x ∈ S gcd(a, x) = 1, and ∃k ≥ 1, a k ∈ I. Then ∀k ≥ 1 such that a k ∈ I, Φ(S ∪ {a k }, I\{x, gcd(x, a) > 1}) ⊆ Φ(S, I). Consequently, ∀M ∈ Φ(S, I), ∃m ∈ M such that a|m.

The above property while obvious is necessary when filtering. The property below is a generalization of this and the lemma 1 of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]:

Property 2 (Primes with limited representations). Let a, b ∈ P, (∀x ∈ S, gcd(ab, x) = 1), k, k ≥ 1 s.t a k b k ∈ I and (∀m ∈ I, a|m =⇒ b|m). Then Φ(S ∪ {a k b k }, I\{x s.t gcd(x, ab) > 1}) ⊆ Φ(S, I)
Property 2 does not make any size arguments nor does it require unicity of a representation: in a sense it is also a specialization of lemma 2 for prime recompositions.

The two above properties does not make use of the strongest feature of the first method which is being able to exclude primes per size argument. However, we can have a close version of it since any finite set has a lower and upper bound. Let P = {p 1 , ..., p s } where p 1 < ... < p s be the remaining primes factors we can use to construct new elements of I that are coprime to our previously selected integers. Then Property 3 (Excluding large elements). p s can be excluded whenever:

• p 1 p s > I max • ∀i ∈ [1, s], p i p s / ∈ I and p 2 1 p s > I max • Let X 3 = {m = p i p j p s s.t ∀i, j ∈ [1, s], m < p 4 1 }. If ∃m ∈ X 3 s.t m > I max and ∀m ∈ X 3 , m < m =⇒ m / ∈ I
Intuitively, if we note X k the set of combinations of k primes of R that are lower that p k+1 1 , we can generalize the property even further as: if for all k < k elements of X k are not in I and X k is as above then p s can be excluded. However, computing those sets X k is not trivial: for X 3 , we only know that p 2 1 < p 1 p 2 ≤ p i p j for i = j but we cannot compare p 1 p 3 and p 2 2 without actually computing their values. To avoid bruteforcing, we do not rely on the full property or its generalization.

Crafting the algorithm parts

The algorithm we craft here will roughly follow the same principles as the three-step method of recomposing by factors, with some slight modifications:

• First step gather all possible primes.

• Second step eliminate both largest and smallest primes simultaneously, iterating pincer-wise.

•

Third step reconstruct all remaining candidates in I.

Let f I be a function that outputs TRUE for an entry x if and only if x ∈ I. Note that I min , I max the integral range can be used to construct f I , even when I is not an interval: any finite set of numbers have a minimal and maximal value. We do not specify any f I to keep a general algorithm, but keep in mind that this whole section relies on such a function to test whether an element is part of any given set I.

Instead of managing all the primes from [2, I {max} ] ∩ P, we should only consider the primes in I P . For example, if I requires every factor to be 1 mod 8 or to be quadratic residues modulo some big prime, we can eliminate a lot of prime factors as preprocessing. However, for random generic sets, I P is not always easily available and we can choose to take P ∩ [2, I max ] instead: the algorithm will be less efficient, but the end result will still allow to find a M ∈ Φ(∅, I).

The reasoning used in previous algorithms are still valid and useful for optimization. Instead of storing all the primes into one set, we will split the storage into two sets for efficiency consideration: P L will contain all the primes within [2, √ I max [, and P H will contain all the primes within

[ √ I max , I max ].
In practice, if we had a complete description of I P , we can possibly reduce those sets. For a, b ∈ P H , we have ab / ∈ I: to create an element of I from a ∈ P H , either a ∈ I or one has to multiply a by elements of P L . The first step is then as described in algorithm 6. Before we explain the second step, we briefly explain the quite simple exhaustive search for the undecided combinations in algorithm 7.

R 3 ← R 3 ∪ {p 3 } 8: p 3 ← p 3 p 9: R 2 ← R 2 ∪ R 3 10: RES = {} 11: for p ∈ R 2 do 12: if f I (p) then RES ← RES ∪ {p}
13: return RES Before we describe the middle procedure we will name NewSecondStep between NewFirstStep and NewThird-Step, we need to explain how we plan to select the valid combinations and when we can terminate NewSecondStep and move on to NewThirdStep. Basically, we also proceed by a step-by-step construct-or/and-reject using the above properties:

1) If p 1 p s > I max , use 3, discard p s and repeat.

2) If there is one unique i such that p i p s ∈ I and p 2 1 p s > I max , then we use property 2 and update S ← S ∪ {p i p s }, discard p 1 , p s and repeat. 3) If ∀i, p i p s / ∈ I and p 2 1 p s > I max , then discard p s and repeat (property 3). If p 2 1 p s ∈ I and p 1 p 2 p s > I max , then discard p 1 , p s and update S ← S ∪ {p 2 1 p s } (property 2). 4) If there is at least 2 distinct values i, j with i < j < s such that p i p s , p j p s ∈ I, check if p 2 1 p s > I max . If it is not the case, we move to NewThirdStep. If it is the case, we need to check if for all valid i at least one of them can be easily computed to verify property 2 1 : we can then update with S ← S ∪{p i p s } and R ∪ R\{p i , p s } and repeat. If not (or there is no easier way than bruteforce), then we move to NewThirdStep.

i.e, among the remaining combinations

m, if p i |m =⇒ ps|m.
Note that easy depends on the properties on the entry: for any generic set there is no easy way.

Hence, NewSecondStep presented in Algorithm 8. Note that this version does not make use of property 2, as we are currently unaware of a way to use it generically. Now that we have introduced in details all three core building blocks, we can build by concatenation the general algorithm GenericFactorFilter in algorithm 9. s ← min(S)

Algorithm 8 NewSecondStep

17:

E 2 ← E 2 ∪ {sb}
1: E 1 , P L , P H ← NewFirstStep(I P , f I ) 2: E 2 , P L , P H ← NewSecondStep(I max , P L , P H , f I ) 3: C ← E 1 ∪ E 2 , R ← P L ∪ P H , RES ← {} 4: RES ← NewThirdStep(I max , R, 1, RES, f I ) 5: return C, RES
If RES is empty, then C is maximal

Fourth method: pick first, prove later

The three previous methods relied on proving the maximality of the result on a step-by-step iteration. We present here a new approach: we cannot prove that our intermediate choice will leave an optimal result, but we can easily verify in some cases that the final result is indeed optimal. In the last method, we chose to separate the remaining primes into two sets P L , P H , after E 1 was constructed. Then we have the simple yet very effective property Property 4 (Ad-hoc proof of maximality).

Let E 1 , P L , P H be the result of NewFirstStep on some initial set of numbers I. Let M ⊂ I be a RNS base. Then

|M| = |E 1 | + |P L | =⇒ M ∈ Φ({}, I)
We can trivially see that after NewFirstStep, only products of at least two distinct primes can form a valid element of I. By construction, elements of P H cannot be powered or combined within each other (as we would have an element larger than I max ). Thus, the largest possible outcome is for every other choice to be combining one element of P L and one element of P H . We assume P L and P H have their elements sorted in increasing order: the largest element is at the last position and the smallest element at the first position. We then provide some heuristic completion algorithms that do not guarantee optimality (but where some results can be proven to be maximal). The philosophy here is simple: construct couples ab ∈ I, regardless of whether or not a or b have multiple possibilities, and giving priority to the largest possible combinations. The algorithms can be constructed by swapping NewSecondStep by GreedyMerge (algorithm 10) when needed.

Algorithm 10 GreedyMerge

CONSTRUCTION OF ELEMENTARY RNS BASES

Base of Pseudo-Mersenne

In this section, we present the use of our basic alternative filtering techniques on Pseudo-Mersenne Numbers. Pseudo-Mersenne numbers are of the form 2 n -c, where c is a relatively small constant, and have been proposed for practical use in [START_REF] Crandall | Method and apparatus for public key exchange in a cryptographic system[END_REF]. We are continuing the work of [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF] for pseudo-Mersenne integers within an interval we denote solely for this section I n = [2 n -c, 2 n ], for simplicity, we keep n even to have a simple computation of c = 2 n/2 for now and only attempt to find the size of a RNS base M ∈ Φ(∅, I n ). We conduct experiments for n ∈ [START_REF] Asif | Highly parallel modular multiplier for elliptic curve cryptography in residue number system[END_REF]64]. Let R n be the set R outputted by SecondStep when the entry set is I n and the primes given were already filtered out by FirstStep. We list the number of primes within R n in Table 1. We can observe from Table 1 that the size of R n is low enough to make ThirdStep inexpensive. Note that R n = ∅ is the result of SecondStep and directly implies that The only cases where ThirdStep did not output an empty set are the following: In the above cases, calling a graph algorithm is clearly unnecessary.

{E 1 ∪ E 2 } ∈ Φ(E 1 ∪ E 2 , I n ) ⊆ Φ(∅, I n ).
We could conjecture that for n ≥ 40 then R n = ∅ in such cases but as the density of prime numbers vanishes as n grows we prefer not to.

We present some results for E 1 , E 2 and d = |M| for M ∈ Φ(∅, I n ) For c small, as c = 256, the amount of elements is not larger than the number of primes to consider, thus the Generic Filtering is a better choice. We present our experimental results on c = 256 in Table 3, listing the maximal size d = |M| for M ∈ Φ(∅, [2 n -256, 2 n ]). [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] in

ITALIC P c i = 2 n -m i c i = 2 n -m i k k' 192 5

6) At least d different values

Once the process is finished, we have obtained valid sQ-RNS bases M, M . We cannot prove our result is optimal, but experimentally we still improve the results given by section 4.4 of [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF]: table 6 shows larger sizes d can be found given (n, k), and table 7 shows smaller k for fixed (n, d).

dQ-RNS

As far as we understood, dQ-RNS was conceived to make base searching easier than sQ-RNS. The condition for dQ-RNS was to obtain a single base M of size 2d such that QR(p, M ) = 1 with moduli m i being square numbers. In that regard, our previous approach still works except we ignore the construction of M and focus on a larger M.

This simplification can lead to different approaches to find dQ-RNS bases. For example, to research two RNS base size of d for dQ-RNS within a set of pseudo-mersenne numbers I (defined by (n, k)) given p, we could:

1) Remove all m ∈ I such that QR(p, m) = 1 from I.

2) Remove all m ∈ I such that m is not square from I.

3) Find, if possible M ∈ Φ(∅, I) with |M | ≥ 2d. 4) Split M into M, M such that |M| = |M | = d
Using our filtering methods on both techniques, we should be able to improve the known scalability of [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] through larger RNS bases. The work of [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] is non-trivial, and obtaining provably maximal bases without an exhaustive search seems to be a hard task at hand, and we leave this as an interesting open question.

SUPPLEMENTARY APPLICATIONS

In this section we provide some examples applications where the third and fourth method can extract RNS bases more efficiently.

Large Interval of many possible bitsizes

It is known that for [2, p] the largest subset of pairwise coprime numbers are basically all the primes until p. We previously treated the cases [2 n -2 n/2 , 2 n ], we are now considering the cases [2 n-2 , 2 n ]: experimentally, the final graph given by the first method on those sets are not empty, simply because unicity of representation cannot be guaranteed in intermediate steps. In particular, we can show in those specific examples that NewFirstStep+GreedyMerge outperforms NewFirstStep+NewSecondStep, while proving maximality of the result using property 4 (i.e no need for further processing).

Large Intersection of many intervals

Let I = I 1 ∩ ... ∩ I k . If I is large, using the generic filtering is not ideal, and the first method does not apply whenever I is not an interval: we can imagine I being a mix of Montgomery-friendly integers and pseudo-mersenne integers of different sizes, or more exotic large sets. In this case, we can take a simple arbitrary example where I is a mix between pseudo-Mersenne numbers of various sizes:

I = {a ∈ [2 k -c, 2 k ] | k ∈ [19, 21], c < 2 16 }
It is very easy to verify whether or not an integer belongs to I, and it also seems that on this example the heuristic attempt using GreedyMerge also directly outputs an optimal RNS basis giving E 2 of size 210 for a total RNS base of 229 moduli (E 1 has size 19). However, there is a limit to this approach over the generic approach: the amount of primes to manage can be much larger than |I|. If that is the case, the Generic Filtering method (algorithm 5) could be preferred.

RNS base coprime with a very large smooth integer

This application is more interesting outside the field of cryptography, where we attempt to compute ring operations over a non-prime moduli that might be smooth. In cryptography we usually have to deal with large primes, thus we usually focus on constructing on RNS bases on small moduli that are guaranteed to be coprime to the large prime (to keep operations valid). Suppose for example we want to simplify operations on random rings Z N using a RNS approach where N 's full factorization is unknown: we then have to create novel RNS bases on the fly for every new N given. We can reuse our algorithms by preselecting primes, excluding small prime numbers that are not coprime with N . If a largest base is not needed but rather a base with a total bitsize of 2 log 2 (N ), finding good RNS basis on the fly could be done more efficiently than with previous methods. is indicating when computational time was too long to be relevant () is indicating when algorithm did not return optimal basis bold is indicating the best algorithm for corresponding set

Comparison

We present in Table 9 computation time of state of the art methods as well as our new proposals. Our test were perfromed on a i7 -8565 CPU based laptop and implemented in MAGMA 2 . We perform our test on 3 types of set of moduli.

• First, we simply grow the size of our set of moduli to show the impact of the size of the initial set on each algorithm.

• Secondly, we keep a fix size, but we increase the size of the moduli themselves to show how it impacts those methods.

• Finally, we change the density of moduli in a set, once again to show its impact.

The first four algorithms (from left to right in the table) guarantee an optimal basis, the last two do not. Note that we do not test Alg 4 as some of the samples are not intervals, but test instead Alg 9 as it is simply a slightly slower generic version.

Table 9 clearly indicates that the best choice:

• when the size of the set is large, is to recompose by factors when available

• when the size of elements is large, is to apply the generic technique for random sets

We note also that our new greedy algorithm allows to often obtain optimal bases in an efficient time if the set is reasonably dense. Previous greedy solutions rarely give an optimal basis, whatever is the type of sets studied.

CONCLUSION

This work expands the scope of possibilities over research on the usage of RNS bases, their efficiency and security applications in two ways: firstly, by presenting more effective methods to obtain larger bases over any set; and secondly, by improving the bases exhibited in the recent literature using the aforementioned methods.

In this extended version, we also presented two novel algorithms, offering more adapted tools to situations the previous works could not handle. One of those algorithms 2. http://magma.maths.usyd.edu.au/magma is heuristic and can prove maximality only under certain conditions, which differs from the previous approaches. Those can also still be improved and provide applications that can apply not just to cryptography but to practical integer computer arithmetic in the broad sense, for example computations under integral rings whenever the moduli are random and extremely large.

Lemma 1 (

 1 [START_REF] Gérard | Contributions to the design of residue number system architectures[END_REF]'s maximality preprocessing).Let I = [I min , I max ] and δ = I max -I min > 0. Given S = {a k b ∈ I s.t a ≤ δ, b > δ or b = 1, {a, b} ⊂ P ∪ {1}}we have Φ(S, I) ⊆ Φ({}, I) and Φ(S, I) = ∅.

  2 16 ∈ I, so eliminate 2 from P and add 2 16 to E 1 3 P = {3, 5, ..., 251} and E 1 = {65287, ..., 65521, 2 16 } 4 No other p ∈ P can be transferred and eliminated End |P | = 53 and |E 1 | = 22

Algorithm 2

 2 SecondStep Input: the remaining primes P , output of FirstStep while P = ∅ do a ← min(P ) and b ← PreviousPrime( Imax a ) if ab ∈ I and a ≤ δ and b > δ then (lemma 1) puts ab into E 2 Remove a, b from P (if b ∈ P ) else s ← a 2 and break. s the smallest product of two primes, is used for selecting elements for R while P = ∅ do a ← min(P ) and b ← PreviousPrime( Imax a ) Remove a, b from P (if b ∈ P ) if ab ∈ I and a ≤ δ and b > δ then (lemma 1) puts ab into E 2 else if ab ∈ I and a < b then store ab in RES if as ≤ I max then store a in R if bs ≤ I max then store b in R Output: E 2 , R, RES Example 4. For I = [2 10 -46, 2 10 ], P = {3, 5, . . . , 41, 43} and E 1 = {983, 991, 997, 1009, 1013, 1019, 1021, 1024}. We find E 2 = {(979 = 11 * 89), (995 = 5 * 199), (1003 = 17 * 59), (1011 = 3 * 337)} with R = {7, 13, 19} and RES = {(989 = 23 * 43)}. The break has occurred for a = 7 where we passed in the second while loop. When P is finally empty, |E 2 | has 4 elements, R has 3 elements and RES one.

Example 5 .

 5 With the example for I = [2 16 -2 8 , 2 16 ], Algorithm SecondStep gives R = {29, 37, 47, 53, 61, 71, 73} and RES = ∅. Luckily here, only one combination within elements of R can fit within I, giving 29 × 37 × 61 = 65453. The maximum clique from a singleton is itself, thus we have actually our maximum clique and E 1 ∪ E 2 ∪ {65453} ∈ Φ({}, I).

1 )

 1 Scan elements of I until such y is found. If no such y exists quit, else go to step 2 2) S ← S ∪ {y} and I ← {m ∈ I s.t gcd(y, m) = 1}, return to step 1.

Algorithm 5

 5 Generic Filtering Input: I = {m i } a set of integers Output: S ⊂ I an RNS base, and leftovers I S 1: S ← {} 2: I S ← I 3: while I S decreases in size do

Example 7 .

 7 Let us take the example of:I = {a ∈ [2 16 -2 8 , 216 ], with c = 216 -a and hw(c) < 4} where hw(c) is the number of 1 in the binary representation of c. We obtain E 1 = {65407, 65519, 65536} and the sets of primes |P L | = 42 and |P H | = 56 whose numbers are smaller than if we do not check the belonging to I P ( where |P L | = 53 and |P H | = 169).

18 :PExample 8 .

 188 H ← P H \{b}, P L ← P L \{s} return E 2 , P L , P H Let us take the example of:I = {a ∈ [2 16 -2 8 , 216 ], with c = 216 -a and hw(c) < 4} where hw(c) is the number of 1 in the binary representation of c. We obtain with NewSecondStep E 2 = {65501 = 17 * 3853, 65503 = 31 * 2113, 65515 = 5 * 13103, 65523 = 3 * 21841, 65531 = 19 * 3449} and the sets of primes are now such |P L | = 37 and |P H | = 13. Then we apply NewThirdStep with R = P L ∪ P H , and K = 1, we obtain RES = {65471 = 7 * 47 * 199, 65527 = 7 * 11 * 23 * 37, 65533 = 13 * 71 2 } Algorithm 9 GenericFactorFilter Input: I max , f I , I P Output: C a set of pairwise coprime numbers, RES potential additions to grow C

  -a and hw(c) < 4} where hw(c) is the number of 1 in the binary representation of c. We obtain with GreedyMerge E 2 = {65501 = 17 * 3853, 65503 = 31 * 2113, 65515 = 5 * 13103, 65523 = 3 * 21841, 65531 = 19 * 3449} and the sets of primes are now such |P L | = 37 and |P H | = 51. Then we apply NewThirdStep with R = P L ∪ P H , and K = 1, we obtain RES = {65471 = 7 * 47 * 199, 65527 = 7 * 11 * 23 * 37, 65533 = 13 * 71 2 }

• I 16 :

 16 65453 = 29 × 37 × 61 • I 20 : 1048207 = 73 × 83 × 173 • I 22 : 4193923 = 73 2 × 787 and 4193993 = 109 2 × 353

  Let us select integers step-by-step: first let us set I = [968, 972, 3328, 1701, 875, 1445, 2873, 539, 493, 1573]. 1) The first integer, 968 has gcd 4 with 3328 and 11 with 539. Since gcd(4, 11) = 1, we cannot select 968. 2) The same goes on for every integer up to 539: their possible gcd cannot be reduced to powers of a same prime. 3) We can however pick 493: it shares 17 as gcd with 1445 and 2873, and 1 with other integers. With this information, we update to S = {493}, I = [968, 972, 3328, 1701, 875, 539, 1573]. 4) We retry from 968 to 1701 with no success. 5) 875 has gcd 7 with both 1701 and 539, but 1 with others.

4, algorithm 5 has a straightforward update of Φ(S, I) ← Φ(S , I ) after each loop iteration. Example 6. With this information we can set S = {493, 875} and I = [968, 972, 3328, 1573]. 6) We ignore 968 again, but pick 972, its gcd is 4 with 968 and 3328 but 1 with others. We update to S = {493, 875, 972} and I = [1573]. 7) We finish with S = {493, 875, 972, 1573} and I = ∅, which completes the search.

  Algorithm 6 NewFirstStep, build initial primes Input: I P , f I s.t f I (a) = TRUE ⇐⇒ a ∈ I Output: Initial base E 1 and unused primes P L , P H . 1: E 1 ← {}, P L ← {}, P H ← {} 2: for a ∈ I P do

	4: 5:	E 1 ← E 1 ∪ {a k } else if a < √ I max then
	6:	P L ← P L ∪ {a}
	7:	else
	8:	P H ← P H ∪ {a}
	9: return E 1 ,P L ,P H
	Algorithm 7 NewThirdStep, a recursive procedure
	Input: I max the set range, R a list of leftover primes and f I
		a function to verify inclusion in I
	Output: RES contains the remaining clique candidates.
	1: R 2 = {1}
	2: for p ∈ R do
	3:	R 3 = {}
	4:	for p 2 ∈ R 2 do
	5:	p 3 ← p 2
	6:

3:

if ∃k ≥ 1 s.t f I (a k ) then while p 3 ≤ I max do 7:

  , build E 2 and filter primes Input: I max , P L , P H two sets of primes, f I to verify inclusion within I Output: E 2 , P L , P H the leftover primes 1: E 2 ← {} 2: while P L = ∅ and P H = ∅ do if |P L | > 1 and ab × min(P L \{a}) ≤ I max then

	3:	a ← min(P L ), b ← max(P H )
	4:
	5:	Break	Too many choices to explore
	6:	else
	7:	S ← {s ∈ P L , f I (sb)}
	8:	if f I (a 2 b) and S ⊆ {a} then
	9:	E 2 ← E 2 ∪ {a 2 b}
	10:
	11:	else if |S| = 0 then
	12:
	14:	Break
	15:	else
	16:

P H ← P H \{b}, P L ← P L \{a} P H ← P H \{b} 13: else if f I (a 2 b) or |S| > 1 then

  Input: f I , P L , P H Output: E 2 pairwise coprimes, P L and P H leftover primes 1: E2 = {} 2: aP L = P L 3: while aP L = ∅ do P L ← P L \{a}, P H ← P H \{b} Let us take the example of:I = {a ∈ [2 16 -2 8 , 216 ], with c = 216

	4:	a ← min(aP L )
	5:	aP H = P H
	6:	while aP H = ∅ do
	7:	b ← max(aP H )
	8:	if f I (ab) then
	9:	E 2 ← E 2 ∪ {ab},
	10:
	11:	break
	12:	aP H ← aP H \{b}
	13:	aP L ← aP L \{a}
	14: return E 2 , P L , P H
	Example 9.

Table 1

 1 Size of Rn, output of SecondStep

	n	16 18 20 22 24 26 28 30 32	34	36 38 40 to 64
	R n	7	0	10 60	6	21 19 13	1	283	0	1	0

Table 2

 2 Set size of E 1 , E 2 , and max RNS set size d for even n ∈[START_REF] Asif | Highly parallel modular multiplier for elliptic curve cryptography in residue number system[END_REF] 64]

	n E 1 \{2 n }	E 2	max size
	16	21	25	48
	24	251	198	450
	32	2931	1851	4783
	40	37798	19856	57655
	48	504634	226507	731142
	56 6920100	2724323	9644424
	64 96798093 34267158 131065252

Table 3

 3 Maximum set size d for m i ∈ [2 n -2 8 , 2 n ] for even n ∈[START_REF] Asif | Highly parallel modular multiplier for elliptic curve cryptography in residue number system[END_REF] 64]

	n 16 18 20 22 24 26 28 30 32 34 36 38 40
	d 48 52 45 46 50 50 46 48 49 50 47 52 47
	n 42 44 46 48 50 52 54 56 58 60 62 64
	d 48 50 50 50 48 48 50 49 48 46 49 46

Table 7

 7 Values c i , c i and their size k for sQ-RNS with d = 4 and same n. Our work is in bold,

Table 8 E

 8 2 's size within I = [2 n-2 , 2 n ] after NewFirstStep constructs E 1

	n	16 17 18 19 20	21	22	23	24
	GreedyMerge	24 32 41 59 81 110 152 204 282
	NewSecondStep	23 31 30 50 64 109 123 183 216

Table 9

 9 Time Comparison Set Factors (Alg 9) Generic (Alg 5) Filtering [29] Clique [29] New Greedy (Alg 10) Greedy [26], [33] [2 9 , 2 10 ]

		<0.01s	0.06s	< 0.01s	> 1h	< 0.01s	(< 0.01s)
	[2 10 , 2 11 ]	<0.01s	0.11s	86s	-	< 0.01s	(< 0.01s)
	[2 11 , 2 12 ]	<0.01s	0.34s	2612	-	0.12s	(0.01s)
	[2 13 , 2 14 ]	0.01s	30s	-	-	0.12s	(0.01s)
	[2 14 , 2 15 ]	8.76s	141s	-	-	0.02s	(0.44s)
	[2 24 -2 12 , 2 24 ]	13.4s	0.26s	0.03s	-	(16.7s)	(0.05s)
	[2 28 -2 12 , 2 26 ]	1257s	0.36s	0.05s	-	(811s)	(0.08s)
	[2 52 -2 12 , 2 64 ]	-	0.24s	156s	-	-	(0.07s)
	[2 512 -2 12 , 2 512 ]	-	0.34s	> 1h	-	-	(0.07s)
	10% of [2 18 , 2 19 ]	0.5s	48s	n/a	-	(0.21)	(1.22s)
	50% of [2 16 , 2 17 ]	0.52s	154s	n/a	-	(0.05s)	(1.7s)
	80% of [2 15 , 2 16 ]	4.9s	248s	n/a	-	0.02s	(1.14s)

Base of Solinas numbers

In this part we exhibit our results for the research of Solinas numbers: while every number can be written as x = x i 2 i where x i ∈ {-1, 0, 1}, Solinas numbers are such that the amount of x i = 0 is low, where w called the weight of x is the number of non-zero x i . [START_REF] Bajard | Selected RNS bases for modular multiplication[END_REF] exhibited how to use those numbers efficiently: in particular, they showed how one can avoid using multiplications using those numbers and replace them by more efficient additions, shifts and substractions. However, [START_REF] Bajard | Selected RNS bases for modular multiplication[END_REF] only exhibited RNS bases with 6 integers.

We use the Generic Algorithm to determine d w the size of the elements of Φ(∅, I n,w ) where I n,w are the elements of [2 n -2 n/2 , 2 n ] with weight w. Our experimental results are presented in 4.

CONSTRUCTION OF SPECIFIC RNS BASES

Base of Montgomery-friendly numbers

This section considers "Montgomery-friendly" numbers: introduced in [START_REF] Hamburg | Fast and compact elliptic-curve cryptography[END_REF], [START_REF] Bos | Fast cryptography in genus 2[END_REF], [START_REF] Bos | Selecting elliptic curves for cryptography: an efficiency and security analysis[END_REF], they exhibit efficient properties whenever Montgomery reduction is involved within the computations. [START_REF] Bajard | Montgomery-friendly primes and applications to cryptography[END_REF] showcased that those numbers can have comparable or better performance and scalability compared to pseudo-Mersenne numbers, focusing on numbers of n-bits of the form

We use the Generic Filtering to further improve their results, showcasing that their approach has better scalability than they previously showed, exhibiting the size d of the elements of Φ(∅, I) while [START_REF] Bajard | Montgomery-friendly primes and applications to cryptography[END_REF]'s found sizes |M| are sometimes two times smaller: see the comparison table in Table 5. Parameters n, k, t are taken straight of section 5.3 of [START_REF] Bajard | Montgomery-friendly primes and applications to cryptography[END_REF].

Base of quadratic residues

Previous examples focused on finding exact elements of Φ(∅, I). In some cases, the problematic is different for finding large RNS bases: coprimality is nothing more than one condition rather than the sole condition for a base to be valid. Such is the case of the RNS bases required by [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF]: they exhibited powerful algorithms, tested on FPGA, which outperforms state-of-the-art algorithms at that time. However, specific bases with further requirements were necessary. While they managed to exhibit basis that were large enough for some practical applications, the scalability of their approach was relatively unknown: researching better bases was left as an open question. [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] proposed two algorithms: sQ-RNS and dQ-RNS, with two specific sets of conditions. While we cannot prove to find optimal bases, we show we can still improve the currently known results.

Q-RNS Definitions

The work of [START_REF] Kawamura | RNS montgomery reduction algorithms using quadratic residuosity[END_REF] relies on a bivariate function (close to the Jacobi symbol), the QR function, which creates asymetric relations between the moduli: it is the main reason simple graph algorithms will not work.

Definition 2 (QR function).

Let a, m ∈ N * and gcd(a, m) = 1.

Then, QR(a, m) = 1 ⇐⇒ ∃x s.t. x 2 = a mod m.

Otherwise, QR(a, m) = 0.

Elements of a basis M have to be of the form m i = 2 n -c. Previously, we either aimed to maximize the size d of M.

Here, we also aim to bound c to c < 2 k if |M| = d is fixed.

sQ-RNS sQ-RNS requires two distinct bases

is equal to 1. The prime p is usually fixed by the (cryptographic) system we plan to use the base for. Since QR(M , M )QR(p, M ) = 1, we believe M should be generated first: p, which is fixed, influences M the same way M influences M . Our strategy is then the following:

First, we aim to construct M minimizing k given values d, p fixed by the (crypto)system verifying:

In order to do so, we increment k from k = 1 until the process can successfully complete: first filter with p, k to verify 1), then apply generic filtering to find moduli verifying 2) and repeat until we have d moduli m i to verify 3).

Once a suitable M is found, M is constructed similarly by incrementing from k = k given M until we can verify the following conditions:


Vasiliki Stergiopoulou

Luca Calatroni

Sébastien Schaub

Laure Blanc-Féraud

3D IMAGE SUPER-RESOLUTION BY FLUOROPHORE FLUCTUATIONS AND MA-TIRF MICROSCOPY RECONSTRUCTION (3D-COL0RME)

Keywords: 3D Super-Resolution, Fluorescence microscopy, Blinking fluorophores, Sparse optimisation, MA-TIRF microscopy

We propose a 3D super-resolution approach to improve both lateral and axial spatial resolution of a thin layer adjacent to the coverslip in Total Internal Reflection Fluorescence (TIRF) imaging applications. Our approach, called 3D-COL0RME (3D -Covariance-based ℓ 0 super-Resolution Microscopy with intensity Estimation) improves both lateral and axial resolution by combining sparsity-based modelling for precise molecule localisation and intensity estimation in the lateral plane with a 3D reconstruction procedure in the axial one using Multi-Angle TIRF (MA-TIRF). Differently from state-of-the-art approaches, 3D-COL0RME does not require specific fluorophores or specific illumination. We validate 3D-COL0RME on simulated MA-TIRF blinking-type data and on challenging real MA-TIRF acquisitions, showing significant resolution improvements.

INTRODUCTION

Fluorescence microscopy is an imaging technique which allows the real-time observation of sub-cellular entities in live samples. Due to light diffraction, however, the spatial resolution normally achieved by means of light microscopy techniques is limited. In standard acquisition settings, the diffraction limit, i.e. the size of the biological entities we are able to distinguish, is approximately equal to 200 in the lateral plane and to 500 in the axial one. Several biological entities of interest, however, have size smaller than these values. The implementation of appropriate super-resolution techniques can thus significantly improve the visualisations of invisible sub-cellular entities, by allowing the understanding of the biological functions happening at the molecular level.

There are many super-resolution techniques appeared in the literature. Some of them are able to gain resolution in 3D. The widely known are 3D-STORM [START_REF] Huang | Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy[END_REF] and 3D-StED [START_REF] Harke | Three-dimensional nanoscopy of colloidal crystals[END_REF] that can achieve very good levels of spatial resolution but at the expense of long acquisition times (especially for StORM), fast photobleaching and they require special fluorophores (es-pecially for StED). The most resolutive method in 3D, to our knowledge, is iPALM [START_REF] Shtengel | Interferometric fluorescent super-resolution microscopy resolves 3d cellular ultrastructure[END_REF] which combines StORM approach and interferometry at the price of a highly complex prototype. On the other side 3D SIM [START_REF] Gustafsson | Three-dimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF] has short acquisition times, so improved temporal resolution, but limited spatial resolution.

In this work, we relax the requirement of special equipment and propose a super-resolution method for live-cell imaging with fine spatio-temporal resolution, which exploits the stochastic temporal fluctuations of standard fluorophores. The diffraction-limited images we aim to process are acquired by a Multi-Angle Total Internal Reflection Fluorescent (MA-TIRF) microscope that benefits from low photobleaching, low photodamage and can excite a broad spectrum of common fluorophores. However, due to the fast decay of the evanescent field produced by total internal reflection, TIRF imaging is limited to a thin layer (500-800 nm) adjacent to the coverslip, as only the fluorophores in this layer are excited.

Compared to standard TIRF microscopy acquisitions [START_REF] Axelrod | Cell-substrate contacts illuminated by total internal reflection fluorescence[END_REF][START_REF] Axelrod | Chapter 7 total internal reflection fluorescence microscopy[END_REF], MA-TIRF further provides depth information as it allows the sample's optical sectioning at different axial levels according to the chosen incident angle of the illumination. Given a stack of MA-TIRF acquisitions, a 3D super-resolved image can be estimated by means of suitable super-resolution methods. For improving the spatial resolution in the lateral plane, we consider in this work the Covariance-based ℓ 0 super-Resolution Microscopy approach with intensity Estimation (COL0RME) [START_REF] Stergiopoulou | COL0RME: Covariance-based ℓ 0 super-resolution microscopy with intensity estimation[END_REF][START_REF] Stergiopoulou | COL0RME: Super-resolution microscopy based on sparse blinking fluorophore localization and intensity estimation[END_REF], while for the axial plane, we consider the MA-TIRF reconstruction algorithm proposed in [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF][START_REF] Soubies | Improving 3D MA-TIRF Reconstruction with Deconvolution and Background Estimation[END_REF]. Our combined method is referred to as 3D-COL0RME.

Our paper is structured as follows: in Section 2, we describe in mathematical terms the 3D MA-TIRF superresolution problem. In Section 3 we discuss the COL0RME and MA-TIRF reconstruction models allowing for improved resolution in the lateral and axial plane, respectively. Finally, in Section 4 we report results computed both on simulated and real data.

PROBLEM FORMULATION

To start with, TIRF microscopy exploits the properties of the evanescent field appearing in a small region adjacent to the interface between two mediums with different refractive indices. When the light propagates from a medium with a high refractive index to a medium with a lower refractive index, there is a critical angle ( ) beyond which the total amount of light is reflected. For angles > the light beam is totally reflected, so the electromagnetic field (evanescent wave) penetrates into the medium with an intensity ( , ) that exponentially decays in the axial direction and is given by:

( , ) = 0 ( ) -( ) , ( ) = (sin 2 ( ) -sin 2 ( )) (1
) where ( ) is the inverse of the penetration depth and is a parameter that depends on the index of the incident medium and the excitation wavelength. See [START_REF] Axelrod | Cell-substrate contacts illuminated by total internal reflection fluorescence[END_REF][START_REF] Axelrod | Chapter 7 total internal reflection fluorescence microscopy[END_REF] for more details.

Given a set of images {x ∈ R 2 } =1 obtained with different incident angles of the illumination beam { > } =1 , we are able to retrieve the elevation in the axial direction ( ) in each pixel f ∈ R 2 × , by considering the following problem:

find f ∈ R 2 × s.t. x = W f, (2)
where W : R

2 × → R
2 is an operator representing the weighted summation of the slices of f, with weights related to the angle and equal to the intensity factors { ( , )} =1 of the evanescent field.

Due to light diffraction, the resolution of the acquisitions is still limited in the lateral plane. Therefore, we propose the computation of highly-resolved images in the lateral plane for each one of the incident angles . By acquiring a sequence of images {g , ∈ R 2 } =1 , for each angle , and by exploiting the independence of the random fluctuations of individual fluorescent molecules, we aim to find x ∈ R 2 , the laterally super-resolved image defined on a -times finer grid, if = , > 1. Mathematically, we can write:

x := 1 =1 u , with g , = SHu , + b + n , , (3)
f ∈ R 2 ×
is thus estimated given a stack of blinking acquisitions {{g , } =1 } =1 acquired at different angles > , ∀ ∈ 1, . . . , .

3D-COL0RME

We describe in this section how the two ill-posed inverse problems ( 2) and ( 3) can be solved in a sequential way by means of appropriate sparse regularisation models. Note that by solving (3) a super-resolved image x can be estimated for each incident angle of the illumination beam. Those images serve as input data for problem [START_REF] Harke | Three-dimensional nanoscopy of colloidal crystals[END_REF] where the objective consists in finding the desired 3D image f with improved spatial and axial resolution.

For the resolution of (3) we use the COL0RME approach [START_REF] Stergiopoulou | COL0RME: Covariance-based ℓ 0 super-resolution microscopy with intensity estimation[END_REF][START_REF] Stergiopoulou | COL0RME: Super-resolution microscopy based on sparse blinking fluorophore localization and intensity estimation[END_REF] which is based on the formulation of a two-step procedure which, by solving suitable sparse optimisation problems, computes an accurate estimation of both sample support and intensity. For solving [START_REF] Harke | Three-dimensional nanoscopy of colloidal crystals[END_REF] we use the 3D MA-TIRF reconstruction algorithm proposed in previous works [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF][START_REF] Soubies | Improving 3D MA-TIRF Reconstruction with Deconvolution and Background Estimation[END_REF].

Support and intensity estimation via COL0RME

The method COL0RME [START_REF] Stergiopoulou | COL0RME: Covariance-based ℓ 0 super-resolution microscopy with intensity estimation[END_REF][START_REF] Stergiopoulou | COL0RME: Super-resolution microscopy based on sparse blinking fluorophore localization and intensity estimation[END_REF] is composed of two steps: a support estimation step and an intensity estimation one. The two steps are performed sequentially. The main idea of COL0RME consists in exploiting the temporal and spatial independence of the fluorescent emitters by a sparse approximation of their second-order statistics. A formulation of problem (3) in the covariance domain is the following (see [START_REF] Stergiopoulou | COL0RME: Covariance-based ℓ 0 super-resolution microscopy with intensity estimation[END_REF][START_REF] Stergiopoulou | COL0RME: Super-resolution microscopy based on sparse blinking fluorophore localization and intensity estimation[END_REF] for more details):

r q g = ( )r q u + v I , (4)
where, for all angles , r q g ∈ R 4 is the vectorised form of the covariance matrix of the raw data {g , } =1 , r q u ∈ R Note that due to the MA-TIRF setup, the support of the entire sample can be found by solving (4) only in correspondence with the angle ˜ closest to the critical one , as it will contain information of molecules located in the whole depth of investigation. This corresponds to a significant computational gain as (4) needs thus to be solved only once. By denoting with rg the covariance matrix associated to {g ˜ , } =1 , then the following problem can be considered for solving (4):

(r u , ˆ ) ∈ arg min ru ∈R 2 + , ˜ ∈R+ 1 2 rg -( )r u -˜ v I 2 2 + R (r u ; ) , (5) where
> 0 is a regularization parameter, R (•; ) is a sparsity-promoting penalty, that could be either the ℓ 1 -norm or the continuous exact relaxations of the ℓ 0 (CEL0) [START_REF] Soubies | A continuous exact ℓ 0 penalty (CEL0) for least squares regularized problem[END_REF] penalty. Numerically, problem (5) can be solved by alternate minimisation and by means of standard sparse optimisation solvers (such as FISTA, iterative reweighted ℓ 1 . . . ), see, e.g., [START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF] for a review.

Given ru , the support of interest Ω is simply Ω = { : (r u ) ≠ 0} = : ( x ) ≠ 0 . By an accurate estimation of the noise variance ˆ and by means of the discrepancy principle [START_REF] Hansen | Discrete Inverse Problems: Insight and Algorithms[END_REF], an automatic intensity estimation procedure can be now designed as a second step of COL0RME. In our scenario, for each angle , ∈ 1, ..., , the mean intensity image x restricted to the estimated support Ω and the smoothly varying background b can be estimated from the empirical temporal mean of the acquired stack {g } =1 by solving:

( x , b ) ∈ arg min x∈R |Ω| + , b∈R 2 + 1 2 x -(g -b) 2 2 (6) + 2 ∇ Ω x 2 2 + 2 ∇b 2 2 , ∀ ∈ 1, . . . , .
Here, the parameter can be automatically estimated via discrepancy principle while > 0 does not require very fine tuning. By ∇ ∈ R 2 2 × 2 we denote the discrete gradient operator, while ∈ R 2 × |Ω| is a matrix whose -th column is extracted from for all indexes ∈ Ω and ∇ Ω is the discrete gradient operator restricted to Ω. Problem ( 6) can be solved efficiently via (proximal) gradient-type algorithms.

MA-TIRF reconstruction

Having all the estimated COL0RME images x , ∀ = {1, . . . , } at hand, we can use them to solve the problem (2). To estimate the 3D super-resolved image f ∈ R 2 × , we thus follow [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF][START_REF] Soubies | Improving 3D MA-TIRF Reconstruction with Deconvolution and Background Estimation[END_REF] and look for solutions of

f ∈ arg min f ∈R 2 × + =1 1 2 W f -x 2 2 + R (∇f)) , (7)
where > 0, { x } =1 are the super-resolved COL0RME images and W : R

2 × → R
2 is the discrete TIRF operator related to the angle , while R (∇•) can be the Hessian Shatten-norm or the Total-Variation (TV) regularization, depending on whether Hessian-or gradient-sparsity is seeked, respectively.

RESULTS

Simulated MA-TIRF SOFI data. We start by applying 3D-COL0RME to process simulated 3D tubulin images. To simulate the data, we first set the 3D spatial pattern (see Figure 1a) using the SMLM 2016 MT0 microtubules dataset 1 ; temporal fluctuations are simulated by using the SOFI simulation tool [START_REF] Girsault | SOFI simulation tool: A software package for simulating and testing super-resolution optical fluctuation imaging[END_REF] upon a specific choice of parameters (see below). For five different angles { } 5 =1 of the illumination beam, with < 1 < • • • < 5 , a stack of 500 frames is simulated. The fluctuations' parameters are chosen as: 20ms for on-state average lifetime, 40ms for off-state average lifetime, 35s for average time until bleaching (so that little bleaching, around 18% is practically observed ) and frame rate of 100 frames per second (fps). The PSF used has a full-width-half-maximum (FWHM) of approximately 229nm while the pixel size is chosen to be equal to 100nm. Spatially varying background was added to the acquisition as well as additive Gaussian noise of signal-to-noise ration (SNR) equal to 14.75 dB.

In Figure 2, a single frame of the acquired stack as well as its temporal mean (diffraction-limited image) is shown for each incident angle of the illumination beam in the first and second line, respectively. In the third line, the super-resolved images { x } 5 =1 computed by solving COL0RME models (5)-( 6) at each angle with a super-resolution factor of 4. Finally, in Figure 3b we show the 3D MA-TIRF reconstruction result obtained by solving (7) using { x } 5 =1 as input. The regularization parameters that do not allow automatic estimation have been chosen empirically, while ℓ 1 -and TVregularization has been used for the problems ( 5) and ( 7) respectively.

For comparison, we further plot in Figure 3a the result of the MA-TIRF approach without deconvolution previously considered [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF][START_REF] Soubies | Improving 3D MA-TIRF Reconstruction with Deconvolution and Background Estimation[END_REF]: we can clearly observe that, compared to 3D-COL0RME, shows fine axial, but poor spatial resolution. Real MA-TIRF data. We further apply 3D-COL0RME on a dataset of images of tubulins stained with Alexa Fluor 488 that are acquired by a real MA-TIRF microscope. A sequence of = 500 frames is processed for each , = 0, . . . , 4. As 0 < , no propagation of the evanescent wave is observed for it. However, such angle is used only for a more precise support estimation. Namely, as a final support we consider the superposition of the supports estimated for angles 0 and 1 . The total acquisition time of the whole dataset is approximately 2min, the pixel size of the CCD camera used is 106nm and the FWHM of the PSF has been measured experimentally and is equal to 292nm.

In Figure 4 (first row), the diffraction limited images are shown, together with the spatially super-resolved COL0RME images (second row). Finally, in Figure 5, a comparison between the 3D reconstruction computed using 3D-COL0RME (left upper part) and the standard MA-TIRF approach with only background removal [START_REF] Soubies | Improving 3D MA-TIRF Reconstruction with Deconvolution and Background Estimation[END_REF] (right lower part) is shown.

2

  is the vector of the auto-covariances of the high-resolution images {u , } =1 , ∈ R + is the (unknown) Gaussian noise variance, := SH ∈ R 2 × 2 , denotes the column-wise Kronecker product and v I ∈ R 4 is the vectorised form of the identity matrix I M 2 .

Fig. 1 :

 1 Fig. 1: (a) Ground truth tubulin image with depth information, (b) The exponential decay of the global intensity of the diffraction limited and the super-resolved COL0RME images with respect to the inverse of the penetration depth { ( )} 5 =1 , see (1) .

Fig. 2 :

 2 Fig. 2: For each { }, = 1, . . . , 5: (first row) one frame of the acquired stack g , , (second row) the temporal mean ḡ , (third line) 2D COL0RME results x .

Fig. 3 :

 3 Fig. 3: (a) MA-TIRF reconstruction result (b) Super-resolved 3D-COL0RME reconstruction. Colour quantifies sample depth.

Fig. 4 :

 4 Fig. 4: For each , = 1, . . . , 4 of the illumination beam: (first row) the temporal mean ḡ of the stack, (second row) the 2D COL0RME results x .

Fig. 5 :

 5 Fig. 5: (a) Super-resolved 3D-COL0RME image, (b) Standard MA-TIRF reconstruction. Colour quantifies sample depth.

5 .

 5 COMPLIANCE WITH ETHICAL STANDARDSThis work was conducted using biological data available in open access by EPFL SMLM datasets and data from cultured cells where no ethical questions are concerned.

ACKNOWLEDGEMENTS

The authors are thankful to E. Soubies for providing code for MA-TIRF reconstruction. VS and LBF are supported by the 3IA Côte d'Azur Investments (with reference number ANR-19-P3IA-0002). LC acknowledges the EU H2020 RISE NoMADS grant 777826 and the SPLIN grant by GdR ISIS. Support for development of the microscope was received from IBiSA (Infrastructures en Biologie Santé et Agronomie) to the MICA microscopy platform.