Reference / Paper · 2010
Factorization of a 768-bit RSA Modulus (IACR ePrint 2010/006)
The RSA-768 factorisation by the number field sieve.
Reference / Paper · 2010
The RSA-768 factorisation by the number field sieve.
Factorization of a 768-bit RSA modulus version 1.4, February 18, 2010 Thorsten Kleinjung1 , Kazumaro Aoki2 , Jens Franke3 , Arjen K. Lenstra1 , Emmanuel Thomé4 , Joppe W. Bos1 , Pierrick Gaudry4 , Alexander Kruppa4 , Peter L. Montgomery5,6 , Dag Arne Osvik1 , Herman te Riele6 , Andrey Timofeev6 , and Paul Zimmermann4 1 EPFL IC LACAL, Station 14, CH-1015 Lausanne, Switzerland NTT, 3-9-11 Midori-cho, Musashino-shi, Tokyo, 180-8585 Japan 3 University of Bonn, Department of Mathematics, Beringstraße 1, D-53115 Bonn, Germany 4 INRIA CNRS LORIA, Équipe CARAMEL - bâtiment A, 615 rue du jardin botanique, F-54602 Villers-lès-Nancy Cedex, France 5 Microsoft Research, One Microsoft Way, Redmond, WA 98052, USA 6 CWI, P.O. Box 94079, 1090 GB Amsterdam, The Netherlands 2 Abstract. This paper reports on the factorization of the 768-bit number RSA-768 by the number field sieve factoring method and discusses some implications for RSA. Keywords: RSA, number field sieve 1 Introduction On December 12, 2009, we factored the 768-bit, 232-digit number RSA-768 by the number field sieve (NFS, [20]). The number RSA-768 was taken from the now obsolete RSA Challenge list [38] as a representative 768-bit RSA modulus (cf. [37]). This result is a record for factoring general integers. Factoring a 1024-bit RSA modulus would be about a thousand times harder, and a 768-bit RSA modulus is several thousands times harder to factor than a 512-bit one. Because the first factorization of a 512-bit RSA modulus was reported only a decade ago (cf. [7]) it is not unreasonable to expect that 1024-bit RSA moduli can be factored well within the next decade by an academic effort such as ours or the one in [7]. Thus, it would be prudent to phase out usage of 1024-bit RSA within the next three to four years. The previous record NFS factorization was that of the 663-bit, 200-digit number RSA-200 (cf. [4]), established on May 9, 2005. That 663-bit and the current 768-bit record NFS factorizations should not be confused with record special NFS (SNFS) factorizations, the most recent of which was the complete factorization of the 1039-bit number 21039 − 1, obtained in the spring of 2007 (cf. [2]). Although 21039 − 1 is much bigger than RSA-768, its special form made it an order of magnitude easier to factor than RSA-768. The following effort was involved. We spent half a year on 80 processors on polynomial selection. This was about 3% of the main task, the sieving, which was done on many hundreds of machines and took almost two years. On a single core 2.2 GHz AMD Opteron processor with 2 GB RAM, sieving would have taken about fifteen hundred years. This included a generous amount of oversieving, to make the most cumbersome step, the matrix step, more manageable. Preparing the sieving data for the matrix step took a couple of weeks on a few processors, the final step after the matrix step took less than half a day of computing, but took about four days of intensive labor because a few bugs had to be fixed. It turned out that we had done about twice the sieving strictly necessary to obtain a usable matrix, and that the extra data allowed generation of a matrix that was quite a bit easier than anticipated at the outset of the project. Although we spent more computer time on the sieving than required, sieving is a rather laid back process that, once running, does not require much care beyond occasionally restarting a machine. The matrix step, on the other hand, is a more subtle affair where a slight disturbance easily causes major trouble, in particular if the problem is by its sheer size stretching the available resources. Thus, our approach to overspend on an easygoing part of the computation led to a matrix that could be handled relatively smoothly, thereby saving us considerable headaches. More importantly, and another reason behind the oversieving, the extra sieving data allow us to conduct various experiments aimed at getting a better understanding about the relation between sieving and matrix efforts and the effect on NFS feasibility and overall performance. This is ongoing research, the results of which will be reported elsewhere. All in all, the extra sieving cycles were well spent. For 21039 −1 the matrix step was, for the first time, performed on a number of different clusters: as reported in [2], the main part of the calculation was run in parallel at two locations on four clusters. This is possible due to our usage of the block Wiedemann algorithm for the matrix step (cf. [9]), where the main calculation consists of two consecutive sequences of matrix times vector multiplications. As discussed in [2], a greater degree of independent parallelization leads to a variety of challenges that need to be overcome before larger problems can be tackled (such as, for instance, 1024-bit moduli). Here we set a first step toward a scalable solution by showing how greater flexibility can be achieved in the number and contributions of the various independent clusters. As a result a nine times harder (than for 21039 −1) matrix step was solved in less than four months running up to eight independent jobs on clusters located in France, Japan, and Switzerland. During one of the substeps about a terabyte of memory was needed. These figures imply that much larger matrices are already within reach, leaving preciously little doubt about the feasibility by the year 2020 of a matrix required for a 1024-bit NFS factorization. As part of the experiments mentioned above we also intend to study if a single large cluster would be able to handle such matrices using the block Lanczos algorithm (cf. [8]). Compared to block Wiedemann this has advantages (a shorter, single sequence of iterations and no tedious and memory-hungry central Berlekamp-Massey step [41]) but disadvantages as well (it cannot be run on separate clusters and each iteration consists of a multiplication by both a matrix and its transpose). The factorization problem and the steps taken to solve it using NFS are described in Section 2. The factorization is presented in Section 2.5. Implications for moduli larger than RSA-768 are briefly discussed in the concluding Section 3. Appendix A presents some of the details of the sieving approach that we used, and Appendix B describes a new twist of the block Wiedemann algorithm that makes it easier to share large calculations among different parties. There are many aspects to an effort as described in this paper that are not interesting from a scientific or research point of view, but that are nevertheless crucial for its success. Although NFS can be run as a BOINC project (cf. [31]), we chose not to do so. Instead, we gathered a small but knowledgeable group of contributors who upfront committed themselves to finish the project, while dedicating at least a known minimum of computational resources for a substantial period of time. This allowed us to target a reasonable completion date, a date 2 that we managed to meet easily despite our ample oversieving. Although different sievingclients do not need to communicate, each client needs to communicate a fair amount of data to the central storage location. This required a bit more organizational efforts than expected, occasional recovery from mishaps such as unplugged network cables, switched off servers, or faulty raids, and a constantly growing farm of backup drives. We do not further comment on these more managerial issues in this article, but note that larger efforts of this sort would benefit from full-time professional supervision. 2 2.1 Factoring RSA-768 Factoring using the Morrison-Brillhart approach A composite integer n can be factored by finding integer solutions x, y of the congruence of squares x2 ≡ y 2 mod n, and by hoping that n is factored by writing it as a product gcd(x − y, n) · gcd(x + y, n). For a random such pair the probability is at least 12 that a non-trivial factor of n is found in this way. The Morrison-Brillhart approach [26] solves the equation x2 ≡ y 2 mod n by combining congruences of smooth squares. We explain this in the present section and, very rudimentarily, how NFS applies the same basic idea. A non-zero integer u is b-smooth if the prime factors of |u| are at most b. Each b-smooth integer u corresponds to the (π(b) + 1)-dimensional integer vector v(u) of exponents of the primes ≤ b in its factorization, where π(b) is the number of primes ≤ b and the “+1” accounts for inclusion of the exponent of −1. The factor base consists of the primes at most equal to the smoothness bound. Let n be a composite integer, b a smoothness bound, and t a positive integer. Suppose that π(b) + 1 + t different integers v have been found such that the least absolute remainders r(v) = v 2 mod n are b-smooth. Because the corresponding (π(b) + 1)-dimensional vectors v(r(v)) are linearly dependent, at least t independent Psubsets S of the v’s can be found using linear algebra such that, for each of these is an all-even vector. Each Q subsets, v∈Spv(r(v)) Q subset S then leads to a solution x = v∈S v and y = r(v) to x2 ≡ y 2 mod n and v∈S thus, overall, t chances of at least 21 to factor n. If v’s as above are generated by picking random integers and hoping that the least absolute remainder modulo n of their square is smooth, the search for the proper v’s depends on the smoothness probability of integers of the same order of magnitude as n. That is Dixon’s random squares method (cf. [12]), which has the advantage that its expected runtime can rigorously be proved. It was preceded by Morrison and Brillhart (cf. [26]) who had already obtained a much higher smoothness probability, and thus a much faster – but heuristic – factoring method, by using continued fractions to generate quadratic residues modulo n of roughly half the size of n, i.e., of order n1/2 . Richard Schroeppel with his linear sieve was the first, in about 1976, to combine essentially the same high smoothness probability with a fast sieving-based method to recognize the smooth numbers (cf. [34, Section 6]). He was also the first to analyze the resulting expected runtime (cf. [18, Section 4.2.6]). Later this culminated in Carl Pomerance’s quadratic sieve (cf. [34,35]) which, like Schroeppel’s method, relies on fast sieving based recognition of smooth quadratic residues of rough order n1/2 . 3 Factoring methods of this sort that rely on smoothness of residues of order nθ(1) can be shown to be stuck at expected runtimes of the form e(c+o(1))(ln n) 1/2 (ln ln n)1/2 , for positive constants c and asymptotically for n → ∞. The number field sieve [20] was the first, and so far the only, practical factoring method to break through the barrier of the ln n-exponent of 21 . This is achieved by looking at more contrived congruences that involve smoothness of numbers of order no(1) , for n → ∞, that can, as usual, be combined into a congruence of squares x2 ≡ y 2 mod n. As a result, the number field sieve factors a composite integer n in heuristic expected time 1/3 +o(1))(ln n)1/3 (ln ln n)2/3 e((64/9) , asymptotically for n → ∞. It is currently the algorithm of choice to factor numbers without special properties, such as those on the RSA Challenge List [38]. Of this list, we concentrate on the number RSA-768, a 768-bit RSA modulus with 232-digit decimal representation 123018668453011775513049495838496272077285356959533479219732245215172640050726 365751874520219978646938995647494277406384592519255732630345373154826850791702 6122142913461670429214311602221240479274737794080665351419597459856902143413. Similar to Schroeppel’s linear sieve, the most important steps of NFS are sieving and the matrix step. In the former relations are collected, congruences involving smooth values similar to the smooth r(v)-values above. In the latter linear dependencies are found among the exponent vectors of the smooth values. Unlike Schroeppel’s method, however, NFS requires two nontrivial additional steps: a pre-processing polynomial selection step before the sieving can start, and a post-processing square root step to convert the linear dependencies into congruences of squares. A rough operational description of these steps as applied to RSA-768 is given in the remainder of this section. For an explanation why these steps work, the reader may consult one of the many expositions on NFS in the literature (cf. [20,21,36]). 2.2 Polynomial selection Let n be a composite integer to be factored. In NFS, relations are given by coprime pairs of integers (a, b) with b > 0 such that two integer values v1 (a, b) and v2 (a, b) that depend on a, b are simultaneously smooth, v1 (a, b) with respect to some bound b1 and v2 (a, b) with respect to some bound b2 . The values v1 (a, b) and v2 (a, b) are defined as follows. Let f1 (X), f2 (X) ∈ Z[X] be two irreducible integer polynomials of degrees d1 and d2 , respectively, with a common root m modulo n, i.e., f1 (m) ≡ f2 (m) ≡ 0 mod n. For our simplified exposition we assume that f1 and f2 are monic, despite the fact that the actual f1 and f2 will not be monic. Then vk (a, b) = bdk fk (a/b) ∈ Z. Sufficiently many more than π(b1 ) + π(b2 ) + 2 relations lead to enough chances to factor n (cf. [20]), as sketched in the next two paragraphs. Let Q(αk ) = Q[X]/(fk (X)) for k = 1, 2 be two algebraic number fields. The elements a−bαk ∈ Z[αk ] have norm vk (a, b) and belong to the first degree prime ideals in Q(αk ) of (prime) norms equal to the prime factors of vk (a, b). These prime ideals in Q(αk ) correspond bijectively to the 4 pairs (p, r mod p) where p is prime and fk (r) ≡ 0 mod p: excluding factors of fk ’s discriminant, a first degree prime ideal corresponding to such a pair (p, r mod p) has norm p and is generated by p and r − αk . Because fk (m) ≡ 0Pmod n, there are two natural homomorphisms φk : Z[αk ] → Z/nZ for Pring k −1 k −1 k = 1, 2 that map di=0 ai αki (with ai ∈ Z) to di=0 ai mi mod n and for which φ1 (a−bα1 ) ≡ φ2 (a − bα2 ) mod n. Finding a linear dependency modulo 2 among the vectors consisting of the exponents of the primes in the b1 -smooth v1 (a, b), b2 -smooth v2 (a, b) pairs (and after, at some stage, slightly extending the vectors to make sure that the square norms lead to squares Q of algebraic numbers, cf. [1]), subsets S of the set of relations are constructed such that (a,b)∈S (a − bαk ) is a square σk in Q(αk ), both for k = 1 and for k = 2. With φ1 (σ1 ) ≡ 1/2 φ2 (σ2 ) mod n it then remains to compute square roots τk = σk a solution x = φ1 (τ1 ) and y = φ2 (τ2 ) to x2 ≡ y 2 mod n. ∈ Q(αk ) for k = 1, 2 to find It is easy to find polynomials f1 and f2 that lead to smoothness of numbers of order no(1) , for 3 ln n 1/3 , let d2 = 1, let m be an integer slightly smaller n → ∞. Let d1 ∈ N be of order ( ln ln n ) P d1 i 1/d 1 than n , and let n = i=0 ni m with 0 ≤ ni < m be the radix m representation of n. Then Pd1 f1 (X) = i=0 ni X i and f2 (X) = X − m have coefficients that are no(1) for n → ∞, they have m as common root modulo n, and for the particular choice of d1 the values a, b that suffice to generate enough relations are small enough to keep bd1 f1 (a/b) and bd2 f2 (a/b) of order no(1) as well. Finally, if f1 is not irreducible, it can be used to directly factor n or, if that fails, one of its factors can be used instead of f1 . If d1 > 1 and d2 = 1 we refer to “k = 1” as the algebraic side and “k = 2” as the rational side. With d2 = 1 the algebraic number field Q(α2 ) is simply Q, the first degree prime ideals in Q are the regular primes and, with f2 (X) = X − m, the element a − bα2 of Z[α2 ] is a − bm = v2 (a, b) ∈ Z with φ2 (a − bα2 ) = a − bm mod n. Although with the above polynomials NFS achieves the asymptotic runtime given in the previous section, the amount of freedom in the choice of m, f1 , and f2 means that from many different choices the best may be selected. Here we keep interpretation of the term “best” intuitive and make no attempt to define it precisely. When comparing different possibilities, however, it is usually obvious which one is better. Most often this is strongly correlated with the smoothness probability over the relevant range of a, b pairs, and thus with the size of the coefficients and the number of roots modulo small primes, smoothness properties of the leading coefficients, and the number of real roots. Sieving experiments may be carried out to break a tie. Given this somewhat inadequate explanation of what we are looking for, the way one should be looking for good candidates is the subject of active research. Current best methods involve extensive searches, are guided by experience, helped by luck, and profit from patience. Only one method is known that produces two good polynomials of degrees greater than one (namely, twice degree two, cf. [5]), but its results are not competitive with the current best d1 > 1, d2 = 1 methods which are all based on refinements described in [16] of the approach from [25] and [27] as summarized in [7, Section 3.1]. A search of three months on a cluster of 80 Opteron 3 cores (i.e., 12 · 80 = 20 core years), conducted at BSI in 2005 already and thus not including the idea from [17], produced three pairs of polynomials of comparable quality. We used the 5 pair f1 (X) = 265482057982680X 6 + 1276509360768321888X 5 − 5006815697800138351796828X 4 − 46477854471727854271772677450X 3 + 6525437261935989397109667371894785X 2 − 18185779352088594356726018862434803054X − 277565266791543881995216199713801103343120, f2 (X) = 34661003550492501851445829X − 1291187456580021223163547791574810881 from which the common root m follows as the negative of the ratio modulo RSA-768 of the constant and leading coefficients of f2 . The leading coefficients factor as 23 · 32 · 5 · 72 · 11 · 17 · 23 · 31 · 112 877 and 13 · 37 · 79 · 97 · 103 · 331 · 601 · 619 · 769 · 907 · 1 063, respectively. The discriminant of f1 equals 212 · 32 · 52 · 13 · 17 · 17 722 398 737 · c273, where c273 denotes a 273-digit composite number that probably does not have a factor of 30 digits or less and that is expected to be square-free. The discriminant of f2 equals one. A renewed search at EPFL in the spring of 2007 (also not using the idea from [17]), after we had decided to factor RSA-768, produced a couple of candidates of similar quality, again after spending about 20 core years. Following the algorithm from [16], polynomials with the following properties were considered: the leading coefficient of f2 allowed 11 (first search at BSI) or 10 (second search at EPFL) prime factors equal to 1 mod 6 with at most one other factor < 215.5 , and the leading coefficient of f1 was a multiple of 258 060 = 22 ·3·5·11·17·23. Overall, at least 2·1018 pairs of polynomials were considered. Given the polynomials, smoothness bounds have to be selected. Also, we need to specify a sieving region S of Z × Z>0 where the search for relations will be conducted. We conclude this section with a discussion of this final part of the parameter selection process. During any practical search for smooth numbers, one always encounters near misses: candidates that are smooth with the exception of a small number of prime factors that are somewhat larger than the smoothness bound. Not only can these large primes, as these factors are commonly referred to, be recognized at no or little extra cost, they also turn out to be useful in the factorization process. As a result the smoothness bounds can be chosen substantially lower than they would have to be chosen otherwise and relation collection goes faster while requiring less memory. On the negative side, it is somewhat harder to decide whether enough relations have been found – as the simple criterion that more than π(b1 ) + π(b2 ) + 2 are needed is no longer adequate – but that is an issue of relatively minor concern. In principle the decision requires duplicate removal, repeated singleton removal, and counting, all briefly touched upon at the end of Section 2.3. In practice it is a simple matter of experience. Let b` denote the upper bound on the large primes that we want to keep, which we assume to be the same for the algebraic and rational sides. Thus, b` ≥ max(b1 , b2 ). We say that a non-zero integer x is (bk , b` )-smooth if with the exception of, say, four prime factors between bk and b` , all remaining prime factors of |x| are at most bk . Consequently, we change the definition of a relation into a coprime pair of integers a, b with b > 0 such that bd1 f1 (a/b) is (b1 , b` )-smooth and bd2 f2 (a/b) is (b2 , b` )-smooth. 6 Selection of the smoothness bounds b1 , b2 , and b` , of the search area S for a, b, and how many primes between b1 or b2 and b` will be permitted and looked for is also mostly guided by experience, and supported by sieving experiments to check if the yield is adequate. For RSA-768 we used smoothness bounds b1 = 11 · 108 , b2 = 2 · 108 and b` = 240 on cores with at least 2 GB RAM (which turned out to be the majority). On cores with less memory (but at least a GB RAM, since otherwise it would not be worth our trouble), b1 = 4.5 · 108 and b2 = 108 were used instead (with the same b` ). For either choice it was expected that it would suffice to use as sieving region the bounded subset S of Z × Z>0 containing about 11 · 1018 coprime pairs a, b with |a| ≤ 3 · 109 · κ1/2 ≈ 6.3 · 1011 and 0 < b ≤ 3 · 109 /κ1/2 ≈ 1.4 · 107 . Here κ = 44 000 approximates the skewness of the polynomial f1 , and serves to approximately minimize the largest norm v1 (a, b) that can be encountered in the sieving region. Although prime ideal norms up to b` = 240 were accepted, the parameters were optimized for norms up to 237 . Most jobs attempted to split after the sieving algebraic and rational cofactors up to 2140 and 2110 , respectively, only considering the most promising candidates (cf. [15]). As far as we know, this was the first NFS factorization where more than three large primes were allowed on the algebraic side. 2.3 Sieving With the parameters as chosen above, we need to find relations: coprime pairs of integers (a, b) ∈ S such that bd1 f1 (a/b) is (b1 , b` )-smooth and bd2 f2 (a/b) is (b2 , b` )-smooth. In this section we describe this process. A prime p dividing fk (r) is equivalent to (r mod p) being a root of fk modulo p (as already noted above with the first degree prime ideals, and disregarding prime factors of fk ’s discriminant). Because d2 = 1, the polynomial f2 has one root modulo p for each prime p not dividing the leading coefficient of f2 , and each such prime p divides f2 (j) once every p consecutive values of j. For f1 the number of roots modulo p may vary from 0 (if f1 has no roots modulo p) to d1 (if f1 has d1 linear factors modulo p). Thus, some primes do not divide f1 (j) for all j, whereas other primes p may divide f1 (j) a total of d1 times for every p consecutive j-values. All (prime,root) pairs with prime at most b1 for f1 and at most b2 for f2 can be precomputed. Given a (prime,root) pair (p, r) for fk , the prime p divides bdk fk (a/b) if a/b ≡ r mod p, i.e., for a fixed b whenever a = rb mod p. This leads to the line sieving approach where per bvalue for each (p, r) pair for f1 (or f2 ) one marks with “p” all relevant a-values (of the form rb + ip for i ∈ Z). After this sieving step, the good locations are those that have been marked for many different p’s. These locations are remembered, after which the process is repeated for f2 (or f1 ), and the intersection of the new good locations with the old ones is further studied, since those are the locations where relations may be found. With 1.4 · 107 b-values to be processed (the lines) and about 2 · κ times that many a-values per b, the lines can be partitioned among almost any number of different processors a factoring effort of this limited scale can get access to. This is how the earliest implementations of distributed NFS sieving worked. A more efficient but also more complicated approach has gained popularity since the mid 1990s: the lattice sieve as described in [33]. For a (prime,root) pair q = (q, s) define Lq as the lattice of integer linear combinations of the 2-dimensional integer (row-)vectors (q, 0), (s, 1) ∈ Z2 , and let Sq be the intersection of S and Lq . Fix a (prime,root) pair q = (q, s) for, say, f1 . The 7 special prime q (as this lattice defining large prime was referred to in [33]) is typically chosen of the same order of magnitude as b1 , possibly bigger – but smaller than b` . It follows that q divides bd1 f1 (a/b) for each (a, b) ∈ Sq . Lattice sieving consists of marking, for each other (prime,root) pair p for f1 for which the prime is bounded by b1 , the points in the intersection Lp ∩ Sq . The good locations are remembered, after which the process is repeated for all the (prime,root) pairs for f2 with prime bounded by b2 , and the intersection of the new good locations with the old ones is further studied. In this way relations in Sq are found. Thus, for each of these relations it is the case that q divides v1 (a, b), and the lattice sieve is repeated for other (prime,root) pairs q until enough relations have been found. A relation thus found and involving some prime q 0 6= q dividing v1 (a, b) may also be found when lattice sieving with q0 = (q 0 , s0 ). Thus, unlike line sieving, when lattice sieving duplicates will be found: these have to be removed (while, obviously, keeping one copy). In practice one fixes bounds I and J independent of q and defines Sq = {iu + jv : i, j ∈ Z, −I/2 ≤ i < I/2, 0 < j < J}, where u, v form a basis for Lq that minimizes the norms v1 (a, b) for (a, b) ∈ Sq . Such a basis is found by partially reducing the initial basis (q, 0), (s, 1) for Lq such that, roughly speaking, the first coordinate is about κ times bigger than the second, thereby minimizing the norms v1 (a, b) for (a, b) ∈ Sq , cf. skewness of the sieving region S. Actual sieving is carried out over the set {(i, j) ∈ Z × Z>0 : −I/2 ≤ i < I/2, 0 < j < J}, interpreted as Sq in the above manner. For RSA-768 we used I = 216 and J = 215 which implies that we did lattice sieving over an area of size roughly 231 ≈ 2 · 109 (we did not use any line sieving). With b1 = 11 · 108 and b2 = 2 · 108 , the majority of the primes that are being sieved with can be expected to hit Sq only a few times. Thus, for any p that one will be sieving with, only a few of the j-values (the lines) will be hit, unlike the situation in plain line sieving where the lines are so long that each line will be hit several times by each prime. It implies that, when lattice sieving, one cannot afford to look at all lines 0 < j < J per p, and that a sieving method more sophisticated than line sieving must be used. This sieving by vectors, as it was referred to in [33], was first implemented in [14] and used for many factorizations in the 1990s, such as those reported in [11] and [7]. We used a particularly efficient implementation of sieving by vectors, the details of which have never appeared in print (cf. [13]). They are summarized in Appendix A. Lattice sieving with sieving by vectors as above was done for most of the about 480 million (prime,root) pairs (q, s) for special qs between 450 million and 11.1 billion (and some special qs below 450 million, with a smaller b1 -value) by eight contributing parties during the period August 2007 until April 2009. Table 1 lists the percentages contributed by the eight participants. Scaled to a 2.2 GHz Opteron core with 2 GB RAM, a single (prime,root) pair (q, s) was processed in less than a hundred seconds on average and produced about 134 relations, for an average of about four relations every three seconds. This average rate varies by a factor of about two between both ends of the special q range that we considered. In total 64 334 489 730 relations were collected, each requiring about 150 bytes. Compressed they occupied about 5 terabytes of disk space, backed up at various locations. Of these relations, 27.4% were duplicates. Duplicate removal can simply be done using bucket sorting (on b-values, for instance) followed by application of UNIX commands sort and uniq, but may be complicated by odd data formats, and may be done in any way one deems convenient. We 8 relations matrix contribution stage 1 stage 2 stage 3 total Bonn (University and BSI) 8.14% CWI 3.44% EPFL 29.78% 34.3% 100% 78.2% 51.9% INRIA LORIA (ALADDIN-G5K) 37.97% 46.8% 17.3% 35.0% NTT 15.01% 18.9% 4.5% 13.1% Scott Contini (Australia) 0.43% Paul Leyland (UK) 0.69% Heinz Stockinger (Enabling Grids for E-sciencE) 4.54% Table 1: Percentages contributed (with matrix stages 1, 2, and 3 contributing 35 th, 0, and 25 th to total). name of contributor used hashing. Most uniqueing was done during the sieving. Overall it took less than 10 days on a 2.66 GHz Core2 processor with ten 1TB hard disks. After inclusion of 57 223 462 free relations (cf. [20]) the uniqueing resulted in 47 762 243 404 relations involving 35 288 334 017 prime ideals. Given the set of unique relations, those that involve a prime ideal that does not occur in any other relation, the singletons, cannot be part of a dependency and can thus be removed. Recognizing singletons can be done using a simple hash table with a count (counting “0”, “1”, “many”). Doing this once reduced the set of relations to 28 984 986 047 elements, with 14 498 007 183 prime ideals. However, removal of singletons usually leads to new singletons, and singleton removal must be repeated until no more relations are removed. Note that the first (usually high volume) runs can be made to accept different prime ideals with the same hash to lower memory requirements for the hash table. Once singletons have been removed, we have to make sure that the number of remaining relations is larger than the total number of prime ideals occurring in them, so that there will indeed be dependencies among the exponent vectors. After a few more singleton removals 24 615 168 385 relations involving at most 9 976 671 468 prime ideals were left. At this point further singleton removal iterations were combined with clique removal, i.e., search of combinations with matching first degree prime ideals of norms larger than bk . Removal of almost all singletons together with clique removal (cf. [6]) led to 2 458 248 361 relations involving a total of 1 697 618 199 prime ideals, and still containing an almost insignificant number 604 423 of free relations. This not only clearly indicated that we had done enough sieving, but also gave us lots of flexibility while creating the matrix. In the next section it is described how dependencies were found. Overall singleton removal and clique finding took less than 10 days on the same platform as used for the uniqueing. 2.4 The matrix step Finding dependencies among the rows of a sparse matrix using block Wiedemann or block Lanzcos, currently the two favorite methods, takes time proportional to the product of the dimension and the weight (i.e., number of non-zero entries) of the matrix. Merging is a generic term for the set of strategies developed to build a matrix for which close to optimal performance can be expected for the dependency search of one’s choice. It is described in [6]. We ran about 10 merging jobs, aiming for various optimizations (low dimension, low weight, best-of-both, etc.), which each took a couple of days on a single core per node of a 37 node 2.66 GHz 9 Core2 cluster with 16 GB RAM per node, and a not particularly fast interconnection network. The matrix that we ended up using was generated by a 5-day process running on two to three cores on the same 37-node cluster, where the long duration was probably due to the larger communication overhead. It produced a 192 796 550 ×192 795 550-matrix of total weight 27 797 115 920 (thus, on average 144 non-zeros per row). Storage of the matrix required about 105 gigabytes. When we started the project, we expected dimension about a quarter billion and density about 150, which would have been about 74 times harder than what we managed to find (thanks to our large amount of oversieving). We obtained various other matrices that would all lead to worse performance for the next step, determining dependencies. For instance, for the full collection of relations we found a matrix with about 181 million rows and average weight about 171 per row, which is only slightly worse. Restriction to the set of 1 296 488 663 relations (involving fewer than 1 321 104 619 prime ideals) satisfying the much smaller b` -value of 234 still led, after singleton removal, to more relations (905 325 141) than prime ideals (894 248 046), as required, and ultimately to a 252 735 215×252 734 215-matrix of weight 37 268 770 998. For b` = 235 there were 3 169 194 001 relations on at most 2 620 011 949 prime ideals, of which 2 487 635 021 relations on 1 892 034 766 prime ideals remained after singleton removal, and which led to a 208 065 007 × 208 064 007matrix of weight 31 440 035 830. Neither matrix was competitive. We used the block Wiedemann algorithm [9] with block width 8 · 64 to find dependencies modulo 2 among the approximately 193 million rows resulting from the filtering. The details of this algorithm can be found in [9,41] and [2, Section 5.1]. Below we give a higher level explanation of the three basic steps of the Wiedemann algorithm (cf. [19, Section 2.19]). Given a non-singular d × d matrix M over the finite field F2 of two elements and b ∈ Fd2 , we wish to solve the system M x = b. The minimal polynomial F of M on the vector space spanned by b, M b, M 2 b, . . . has degree at most d, so that F (M )b = d X Fi M i b = 0. i=0 From F0 = 1 it follows that x = Pd i=1 Fi M i−1 b, so it suffices to find the F ’s. i Denoting by mi,j the jth coordinate of the vector M i b, it follows from the above equation that for each j with 1 ≤ j ≤ d the sequence (mi,j )∞ i=0 satisfies a linear recurrence relation of order at most d defined by the coefficients Fi : for any t ≥ 0 and 1 ≤ j ≤ d we have that d X Fi mi+t,j = 0. i=0 It is well known that 2d + 1 consecutive terms of an order d linear recurrence suffice to compute the coefficients generating the recurrence, and that this can be done efficiently using the Berlekamp-Massey method (cf. [23,41]). Each particular j may lead to a polynomial of smaller degree than F , but taking, if necessary, the least common multiple of the polynomials found for a small number of different indices j, the correct minimal polynomial will be found. Summarizing the above, there are three major stages: a first iteration consisting of 2d matrix times vector multiplication steps to generate 2d + 1 consecutive terms of the linear recurrence, the Berlekamp-Massey stage to calculate the Fi ’s, and finally the second iteration consisting 10 of d matrix times vector multiplication steps to calculate the solution using the Fi ’s. For large matrices the two iterations, i.e., the first and the final stage, are the most time consuming. In practice it is common to use blocking, taking advantage of the fact that on 64-bit machines 64 different vectors b over F2 can be processed simultaneously, at little or no extra cost compared to a single vector (cf. [9]), while using the same three main stages. If the vector b̄ is 64 bits wide and in the first stage, the first 64 coordinates of each of the generated 64 bits wide vectors M i b̄ are kept, the number of matrix (M ) times vector (b̄) multiplication steps in both the first and the final stage is reduced by a factor of 64 compared to the number of M times b multiplication steps, while making the central Berlekamp-Massey stage a bit more cumbersome. It is less common to take the blocking a step further and run both iteration stages spread over a small number n0 of different sequences, possibly run on different clusters at different locations; in [2] this was done with n0 = 4 sequences run on three clusters at two locations. If for each sequence one keeps the first 64 · n0 coordinates of each of the 64 bits wide vectors they generate during the first stage, the number of steps to be carried out (per sequence) is further reduced by a factor of n0 , while allowing independent and simultaneous execution on possibly n0 different clusters. After the first stage the data generated for the n0 sequences have to be gathered at a central location where the Berlekamp-Massey stage will be carried out. While keeping the first 64 · n0 coordinates per step for each sequence results in a reduction of the number of steps per sequence by a factor of 64 · n0 , keeping a different number of coordinates while using n0 sequences results in another reduction in the number of steps for the first stage. Following [2, Section 5.1], if the first 64 · m0 coordinates are kept of the 64 d d n0 d bits wide vectors for n0 sequences, the numbers of steps become 64·m 0 + 64·n0 = ( m0 + 1) 64·n0 d 0 and 64·n 0 for the first and third stage, respectively and for each of the n sequences. The choices of m0 and n0 should be weighed off against the cost of the Berlekamp-Massey step 0 0 )3 0 0 )2 with time and space complexities proportional to (m +n d1+o(1) and (m +n d, respectively n0 n0 and for d → ∞, and where the exponent “3” may be replaced by the matrix multiplication exponent (our implementation uses “3”). When actually running the first stage of block Wiedemann in this way using n0 different sequences, the effect of non-identical resources used for the different sequences quickly becomes apparent: some locations finish their allotted work faster than others (illustrated in Fig. 1). To keep the fast contributors busy and to reduce the work of the slower ones (thereby reducing the wall-clock time), it is possible to let a quickly processed first stage sequence continue n0 d for s steps beyond ( m 0 + 1) 64·n0 while reducing the number of steps in another first stage sequence by the same s. As described in Appendix B, this can be done in a very flexible way, as long as the overall number of steps over all first stage sequences adds up to n0 · d n0 (m 0 + 1) 64·n0 . The termination points of the sequences in the third stage need to be adapted accordingly. This is easily arranged for, since the third stage allows much easier and much wider parallelization anyhow (assuming checkpoints from the first stage are kept). Another way to keep all parties busy is swapping jobs, thus requiring data exchanges, synchronization, and more human interaction, making it a less attractive option altogether. For our matrix with d ≈ 193·106 we used, as in [2], m0 = 2n0 . But where n0 = 4 was used in [2], we used n0 = 8, thereby quadrupling the Berlekamp-Massey runtime and doubling its memory compared to the matrix from [2], on top of the increased runtime and memory demands caused by the larger dimension of the matrix. On the other hand, the compute intensive first and 11 (a) First stage contributions. (b) Final shot of third stage bookkeeping tool. Fig. 1: Contributions to sequences 0-7 with blue indicating INRIA, orange EPFL, and pink NTT. third stages could be split up into twice as many independent jobs as before. For the first stage 6 8 on average ( 16 + 1) 193·10 64·8 ≈ 565 000 steps needed to be taken per sequence (for 8 sequences), 6 for the third stage the average was about 193·10 64·8 ≈ 380 000 steps. The actual numbers of steps varied, approximately, between 490 000 and 650 000 for the first stage and between 300 000 and 430 000 for the third stage. The calculation of these stages was carried out on a wide variety of clusters accessed from three locations: a 56-node cluster of 2.2GHz dual hex-core AMD processors with Infiniband at EPFL (installed while the first stage was in progress), a variety of ALADDIN-G5K clusters in France accessed from INRIA LORIA, and a cluster of 110 3.0GHz Pentium-D processors on a Gb Ethernet at NTT. A comprehensive overview of clusters and timings is given in Appendix E. It follows f