A Montgomery-like Square Root for the Number Field Sieve - CiteSeerX
Algorithmic Number Theory { Proceedings of ANTS-III (June 21{25, 1998, Portland, Oregon, USA) J. Buhler (Ed.), vol. 1423 of Lecture Notes in Computer Science, pages 151{168
c Springer-Verlag (http://www.springer.de/comp/lncs/index.html)
A Montgomery-like Square Root for the Number Field Sieve Phong Nguyen [email protected] http://www.dmi.ens.fr/~pnguyen/
Ecole Normale Superieure Laboratoire d'Informatique 45, rue d'Ulm F - 75230 Paris Cedex 05
Abstract. The Number Field Sieve (NFS) is the asymptotically fastest
factoring algorithm known. It had spectacular successes in factoring numbers of a special form. Then the method was adapted for general numbers, and recently applied to the RSA-130 number [6], setting a new world record in factorization. The NFS has undergone several modi cations since its appearance. One of these modi cations concerns the last stage: the computation of the square root of a huge algebraic number given as a product of hundreds of thousands of small ones. This problem was not satisfactorily solved until the appearance of an algorithm by Peter Montgomery. Unfortunately, Montgomery only published a preliminary version of his algorithm [15], while a description of his own implementation can be found in [7]. In this paper, we present a variant of the algorithm, compare it with the original algorithm, and discuss its complexity.
1 Introduction The number eld sieve [8] is the most powerful known factoring method. It was rst introduced in 1988 by John Pollard [17] to factor numbers of form x3 + k. Then it was modi ed to handle numbers of the form re ? s for small positive r and jsj: this was successfully applied to the Fermat number F9 = 2512 + 1 (see [11]). This version of the algorithm is now called the special number eld sieve (SNFS) [10], in contrast with the general number eld sieve (GNFS) [3] which integers. n in heuristic time can handle 1arbitrary GNFS factors1=integers = 3 2=3 3 exp (cg + o(1)) ln n ln ln n with cg = (64=9) 1:9. Let n be the composite integer we wish to factor. We assume that n is not a prime power. Let Zn denote the ring Z=nZ. Like many factoring algorithms, the number eld sieve attempts to nd pairs (x; y) 2 Z2n such that x2 y2 (mod n). For such a pair, gcd(x ? y; n) is a nontrivial factor of n with a probability of P d 1 j at least 2 . The NFS rst selects a primitive polynomial f (X ) = j=0 cj X 2 Z[X ] irreducible over Z, and an integer m with f (m) 0 (mod n). Denote by F (X; Y ) = Y df (X=Y ) in Z[X; Y ] the homogenous form of f . Let 2 C be a
2 root of f , and K = Q ( ) be the corresponding number eld. There is a natural ring homomorphism from Z[ ] to Zn induced by ( ) m (mod n). We will do as if mapped the whole K . If ever ( ) is not de ned for some 2 K , then we have found an integer not invertible in Zn, and thus, a factor N of n which should not be trivial. If n0 = n=N is prime, the factorization is over, and if not, we replace n by n0 , and by 0 induced by 0 ( ) m (mod n0 ). By means of sieving, the several integer pairs (ai ; bi ) and a nite Q NFS nds ? bi ) and Qi2S (aiQ? bi m) are squares in K nonempty set S such that i2S (ai?Q and in Z, respectively. We have i2S (ai ? bi ) i2S (ai ? bi m) (mod n); therefore
13 2s 3 2 0s Y Y 4 @ (ai ? bi )A5 4 (ai ? bi m)5 2
i2S
2
i2S
(mod n)
after extracting the square roots, which gives rise to a suitable pair (x; y). The NFS does root of the Q not specify how to evaluate these square roots. The square prime factorizainteger i2S (ai ? bi m) mod n can be found using the known Q tions of each ai ? bi m. But extracting the square root of i2S (aiQ? bi ) is much more complicated and is the subject of this paper. We note = i2S (ai ? bi ): The following facts should be stressed: { the cardinality jS j is large, roughly equal to the square root of the run time of the number eld sieve. It is over 106 for n larger than 100 digits. { the integers ai; bi are coprime, and t in a computer word. { the prime factorization of each F (ai ; bi) is known. { for every prime number p dividing cd or some F (ai; bi ), we know the set R(p) consisting of roots of f modulo p, together with 1 if p divides cd . The remainder of the paper is organized as follows. In Section 2, we review former methods to solve the square root problem, one of these is used in the last stage of the algorithm. Section 3 presents a few de nitions and results. In Section 4, we describe the square root algorithm, which is a variant of Montgomery's original algorithm, and point out their di erences and similarities. We discuss its complexity in Section 5. Finally, we make some remarks about the implementation in Section 6, and the appendix includes the missing proofs.
2 Former methods UFD method. If is an algebraic integer and the ring Z[ ] is a unique factorization domain (UFD), then each ai ? bi can be factored into primes and units,
and so can be , which allows us to extract a square root of . Unfortunately, the ring Z[ ] is not necessarily a UFD for the arbitrary number elds GNFS encounters. And even though Z[ ] is a UFD, computing a system of fundamental units is not an obvious task (see [4]). The method was nevertheless applied with success to the factorization of F9 [11].
3
Brute-force method. One factorizes the polynomial P (X ) = X ? over
K [X ].
2
To do so, one has to explicitly write the algebraic number , for instance by expanding the product: one thus gets the (rational) coecients of as a polynomial of degree at most d ? 1 in . But there are two serious obstructions: the coecients that one keeps track of during the development of the product have O(jS j) digits. Hence, the single computation of the coecients of can dominate the cost of the whole NFS. And even if we are able to compute , it remains to factorize P (X ). One can overcome the rst obstruction by working with integers instead of rationals: let f^(X ) be the monic polynomial F (X; cd ), and ^ be the algebraic integer cd which is a root of f^. If is a square in K then 0 = c2ddjSj=2e f^0 (^ )2 is a square in Z[^ ], where f^0 denotes the formal derivative of f^. It has integral coecients as a polynomial of degree at most d?1 in ^, and these can be obtained with the Chinese Remainder Theorem, using several inert primes (that is, f is irreducible modulo this prime) if there exist inert primes (which is generally true). This avoids computations with very large numbers. However, one still has to factorize the polynomial Q(X ) = X 2 ? 0 ; whose coecients remain huge, so the second obstruction holds. Furthermore, a large number of primes is required for the Chinese Remainder Theorem, due to the size of the coecients.
Couveignes's method. This method overcomes the second obstruction. If f has odd degree d, Couveignes [5] remarks that one is able to distinguish the two p square roots of any square in K , by specifying its norm. Let 0 be the square root with positive norm. Since the prime factorization of N ( 0 ) is known, the p 0 integer any prime q. If q is inert p N ( ) can be eciently computed modulo (mod q). From the Chinese then 0 (mod q) can be computed after expanding 0 p Remainder Theorem, one recovers the coecients of 0 2 Z[^ ]. One can show that the complexity of the algorithm is at best O(M (jS j) ln jS j), where M (jS j) is the time required to multiply two jS j-bit integers. The algorithm appears to
be impractical for the sets S now in use, and it requires an odd degree. Montgomery's strategy [15, 14, 7] can be viewed as a mix of UFD and bruteforce methods. It bears some resemblance to the square root algorithm sketched in [3] (pages 75-76). It works for all values of d, and does not make any particular assumption (apart from the existence of inert primes) about the number eld.
3 Algebraic preliminaries Our number eld is K = Q ( ) = Q (^ ), where is an algebraic number and ^ = cd is an algebraic integer. Let O be its ring of integers, and I be the abelian group of fractional ideals of O. For x1 ; : : : ; xm 2 K , we note < x1 ; : : : ; xm > the element of I generated by x1 ; : : : ; xm . For every prime ideal p, we denote by vp the p-adic valuation that maps I to Z. We de ne the Q numerator and denominator of I 2 I to be the integral ideals numer(I ) = vp (I )>0 pvp (I ) and Q denom(I ) = vp (I ) in I . We follow the notations of [3] and recall some results. Let R be an order in O. By a \prime of R" we mean a non-zero prime ideal of R. We denote by flp;R : K ! Zgp the unique collection (where p ranges over the set of all primes of R) of group homomorphisms such that: { lp;R(x) 0 for all x 2 R; x 6= 0; { if x is a non-zero element of R, then lp;R(x) > 0 if and only if x 2 p; { for Y eachlpx;R(2x) K one has lp;R(x) = 0 for all but nitely many p, and N (p) = jNK (x)j; where p ranges over the set of all primes of R. p
lp;O (x) coincide with vpP ( < x > ). Let i = cd d?1?i + cd?1 d?2?i + + ci+1 . ?2 i Z is an order of O, which is in fact Z[ ] \ Z[ ?1]. We know that A = Z + id=0 Its discriminant (A) is equal to (f ) and we have: (d?1)(d?2) 2 [
(Z[^ ]) = c(dd?1)(d?2)(A);
[O : Z[^ ]] = cd
O : A]:
Recall that for any prime number p, R(p) is de ned as the set consisting of roots of f modulo p, together with 1 if p divides cd. Note that this R(p) is denoted R0 (p) in [3]. The pairs consisting of a prime number p and an element r 2 R(p) are in bijective correspondence with the rst degree primes p of A: { if r 6= 1 then p is the intersection of A and the kernel of the ring homomorphism p;r : Z[ ] ! Fp that sends to r. { if r = 1 then p is the intersection of A and the kernel of the ring homomorphism p;1 : Z[ ?1] ! Fp that sends ?1 to 0. Let p be a prime number, r an element of R(p) and a, b be coprime integers. If a br (mod p) and r 6= 1, or if b 0 (mod p) and r = 1, we de ne ep;r (a; b) = vp (F (a; b)) where vp denotes the ordinary pQ-adic valuation. Otherwise, we set ep;r (a; b) = 0. We have NK (a ? b ) = c1d p;r pep;r (a;b) ; the product ranging over all pairs p, r with p prime and r 2 R(p). Furthermore, for any coprime integers a, b and any rst degree prime p of A corresponding to a pair p, r 2 R(p), we have: e (a; b) 1 lp;A (a ? b ) = ep;r (a; b) ? v (c ) ifif rr 6= =1 p;r
p d
Theorem 1. Let a and b be coprime integers, and p be a prime number. Let p be a prime ideal of O above p such that vp ( ) 6= 0. If p does not divide [O : A] then:
5 1. For every r 2 R(p), there is a unique prime ideal pr of O that lies over the rst degree prime ideal qr of A corresponding to the pair p, r. pr is a rst degree prime ideal, given by pr = : Furthermore, we have vpr ( ) = lqr ;A (a ? b ). 2. There is at most one nite r 2 R(p) such that ep;r (a; b) 6= 0: 3. If p does not divide cd, such a nite r exists and p = pr . 4. If p divides cd, then either p is p1 , or pr for r nite. 5. p divides F (a; b) or cd . Proof. Let r 2 R(p) and qr be the rst degree prime ideal of A corresponding to the pair p,r. Since P p does not divide [O : A], we have from [3] (Proposition 7.3, pages 65-66): pr jqr f (pr =qr ) = 1; where pr ranges over all primes of O lying over qr and f denotes the residual degree. This proves that pr is unique and is a rst degree prime ideal. From [3] (Proposition 7.2, page 65), we also have:
lqr ;A (a ? b ) =
X
p0 jqr
f (p0 =qr )lp0 ;O (a ? b ) = lpr ;O (a ? b ):
Hence, vpr (a ? b ) = lqr ;A (a ? b ). Moreover, we know a Z-basis for any ideal qr of A, namely (p; 0 ? p;r ( 0 ); : : : ; d?2 ? p;r ( d?2 )): Since pr lies over qr , this Z-basis is a system of O-generators for pr . We therefore proved 1. From the de nition of i , one sees that i = ci ?1 + ci?1 ?2 + + c0 ?i?1 ; which proves that p;1 ( i ) = 0. This simpli es the formula when r = 1. One obtains 2 from the de nition of ep;r . Denote by q the intersection of p and A. q is a prime of A and p lies over q. We have lq;A (a ? b ) 6= 0 since vp (a ? b ) 6= 0. From [3] (page 89), this proves that q is a rst degree prime ideal of A. Hence, there exists r 2 R(p) such that q = qr . From 1, this proves that p = pr . This r is nite or in nite, and if r is nite, it is the r of 2. This proves 3 and 4. From the formula expressing lq;A (a ? b ) in terms of ep;r (a; b), we obtain 5. ut
4 The square root algorithm We that we want to compute a square root of the algebraic number = Q recall ( a ? i2S i bi ): The algorithm is split as follows:
1. Transform in order to make < > simpler. The running time of the rest of the algorithm heuristically depends on C ( < > ). p 2. Compute < > from the prime ideal factorization of < > given by the prime factorization of each F (ai ; bi ). p p 3. Approximate from < > : using lattice reductions, construct a sequence of algebraic Q integers 1; : : : ; L in O and signs s1; : : : ; sL in f1g such that = L`=1 `?2s` is a \small" algebraic integer. Qcan be thought as the square of the \guessing-error" inpthe approximation L`=1 `s` of p . 4. Since is a square, so is . Compute using brute-force method. One is able to explicitly write because is a \small" algebraic integer.
6
We thus obtain p as a product of algebraic integers with exponents 1: L p = p Y s` : `=1
`
This enables to compute (p ) without explicitly calculating p , and hopefully some factors of n. Although formalized di erently, Montgomery's algorithm uses the same strategy. Only the steps change. We use another heuristic approach in Step 1, which seems to be more e ective in practice. We use a new process in Step 2, derived from Section 3. Montgomery used a process which was as ecient, but only heuristic. Step 3 is the core of the algorithm. We modi ed this step by using the integral basis in a systematic manner, instead of the power basis. This simpli es the algorithm and the proofs. Heuristically, this should also improve the performances. We postpone the computation of the error in Step 4, while Montgomery included it in Step 3, by updating the computations during the approximation. This decreases the running-time because it is easier to estimate the necessary computations when Step 3 is over, and sometimes, Step 4 can be avoided (when the approximation is already perfect, which can be checked without additional computations). The new algorithm might be more suited to analysis, but like Montgomery's algorithm, its complexity has yet to be determined, even though they both display signi cantly better performances than former methods.
4.1 Computing in the number eld The ring of integers. During the whole algorithm, we need to work with ideals and algebraic integers. We rst have to compute an integral basis of O.
In general, this is a hopeless task (see [13, 2] for a survey), but for the number elds NFS encounters (small degree and large discriminant), this can be done by the so-called round algorithms [16, 4]. Given an order R and several primes pi , any round algorithm will enlarge this order for all these primes so that the new order Rb is pi -maximal for every pi . If we take for the pi all the primes p such that p2 divides (R), then Rb = O. To determine all these primes, a partial factorization of (R) suces, that is a factorization of the form df 2 where d is squarefree and f is factorized. Theoretically, a partial factorization is as hard to nd as a complete factorization and unfortunately, the discriminant is sometimes much larger than the number n we wish to factor. However, if one takes a \random" large number, and one removes all \small" prime factors from it (by trial division or by elliptic curves [12]), then in practice the result is quite likely to be squarefree. Furthermore, even in the case Rb 6= O, it will be true that Rb has almost all of the good properties of O for all ideals that we are likely to encounter in practice, like the fact that every ideal is a product of prime ideals. This is because every order satis es these properties for all ideals that are coprime to the index of the order in O. Hence, we can now assume that an integral basis (!1 ; : : : ; !d) of O has been computed.
7
Algebraic numbers and ideals. From this integral basis we can represent any algebraic number of K as a vector of Q d : thisPis the integral representation. If x 2 K we de ne x = [x1 ; : : : ; xd ]t where x = di=1 xi !i and xi 2 Q . We can also represent any algebraic number as a polynomial of degree at most d ? 1 in : this is the power representation. When dealing with algebraic integers, the integral representation is preferable. We will represent any integral ideal I by an integral matrix (with respect to (!1 ; : : : ; !d)) from a Z-basis or a system of O-generators. In the case of Z-basis, we use the Hermite normal form (HNF) of the square matrix for eciency reasons. We refer to [4] for algorithms concerning algebraic numbers and ideals. 4.2 Simplifying the principal ideal Q If is p a square in K , then so is any 0 = i2pS (ai ? bi p)ei ; when ei = 1. Since Q p (ai ? bi ); we can recover from 0 but actually, we only
= 0 ei =?1
look for a square identity. Fortunately:
2 s 3 3 2s Y Y 4 ( (ai ? bi )ei )5 4 (ai ? bim)ei 5 2
i2S
2
(mod n)
i2S
p
This replaces the computation of p by the computation of 0 . By cleverly selecting the ei , C ( < 0 > ) will be much smaller than C ( < > ): this is because many share the same prime ideals, since many NK (ai ? bi ) share the same primes (as a consequence of sieving). We now address the optimization problem of selecting the ei so that C ( < 0 > ) is small. Given a distribution of ei , the complexity of < 0 > can be computed by the following formula (which comes from the known \factorization" of each ai ? bi into primes of A):
Y jPi2S eiep;r ai;bi j Y jPi2S ei ep;1 ai;bi ?vp cd j p p : (
p;r6=1
)
[
pjcd
(
)
(
)]
The simplest method is a random strategy which selects randomly ei = 1. Another method is a greedy strategy (used in [7]): at every step, select ei = 1 according to the best complexity (whether we put ai ? bi in the numerator or in the denominator). This behaves better than the random strategy. But the best method so far in practice is based on simulated annealing [18], a well-known probabilistic solution method in the eld of combinatorial optimization. Here, the con guration space is E = f?1; +1gjSj, and the energy function U maps any e = (e1 ; : : : ; ejS j ) 2 E to ln C ( < > ) where corresponds to e. For any e 2 E , we de ne its neighbourhood V (e) = f(e1; : : : ; ei?1 ; ?ei ; ei+1 ; : : : ; ejSj) j i = 1; : : : ; jS jg. We try to minimize U by the following algorithm, which performances depend on three parameters i ; f (initial and nal temperatures) and : { select randomly e 2 E and set ? i . { choose randomly f 2 V (e) and set ? U (f ) ? U (e). If > 0, set p ? exp(?=), otherwise set p ? 1. Then set e ? f with probability p, and ? .
8
{ repeat previous step if > f . Although this method behaves better in practice than previous methods, theoretical estimates can hardly be given.
4.3 Ideal square root
Q
From now on, we forget about the initial and set = i2S (ai ? bi )ei : p We wish to obtain as a product of ideals with exponents lying in Z (this ideal is too large to be represented as a single matrix). This Q can be done by factoring into prime ideals the fractional ideal < > = < i2S (ai ? bi )ei > . We simplify the problem to the factorization of any linear expression < ai ? bi > with coprime ai ; bi . Such a factorization could be obtained by general ideal factorization algorithms (see [4]) but this would be too slow if we had to use these algorithms jS j times. Fortunately, we can do much of the work by ourself using the known factorization of each F (ai ; bi ) = f (ai =bi )bdi , as shown in the previous section. We say that a prime number p is exceptional if p divides the index = [O : A]. Otherwise, we say that p is normal. Naturally, a prime ideal of O is said to be exceptional (resp. normal) if it lies above an exceptional (resp. normal) prime. If m is the number of prime factors of , there are at most md exceptional prime ideals. We compute all the exceptional prime ideals (for example, by decomposing all the exceptional primes in O using the BuchmannLenstra algorithm described in [4]), along with some constants allowing us to compute eciently any valuation at these primes. From Theorem 1, we get the prime ideal factorization of < a ? b > as follows: for every prime number p dividing cd or such that there exists a nite r 2 R(p) satisfying ep;r (a; b) 6= 0, { if p is exceptional, compute the valuation of at all the exceptional ideals lying above p. { otherwise, p is normal. If there is a nite r 2 R(p) such that ep;r (a; b) 6= 0 (r is then unique), pick the prime ideal pr with exponent ep;r (a; b) where = : If 1 2 R(p), also pick the prime ideal p1 with exponent ep;1 (a; b) ? vp (cd ) where p1 = : We thus decompose < > pas a product of ideals where every exponent is necessarily even, which gives < > . Montgomery used a di erent ideal factorization process (see [7, 14]) by introducing a special ideal, but its correctness is not proved. pr
4.4 Square root approximation We now use the ideal square root p < > to approximate p . Since p < >
is a huge ideal, we will get an approximation through an iterative process, by selecting a small part of the ideal at each step: this small part will be alternatively
9 taken in the numerator and denominator. To lift an integral ideal to an algebraic integer, we use lattice reduction techniques. We associate several variables at each step `: { an algebraic number `. It canp be considered as the square of the error in the current approximation of . { a sign s` in f?1; +1g, indicating whether we take something in the denominator or in the numerator of the huge original ideal. { a fractional ideal G`, which is an approximation to p < ` > . p { an integral ideal H` of bounded norm. It di erentiates G` from < ` > . { an algebraic integer `. { an integral ideal I` of bounded norm. Q We initialize these variables by: 1 = = i2S (ai ? bi )ei , G1 = p < > , H1 = < 1 > , s1 = 1 if NK ( ) 1 and ?1 otherwise. Each step of the approximation makes `+1 in some sense smaller than ` , and G`+1 simpler than G` . After enough steps, G` is reduced to the unit ideal < 1 > , and ` becomes an algebraic integer suciently small that its integral representation can be determined explicitly (using Chinese Remainders) and a square root constructed using brute-force method. At the start of step `, we need to know the following: { approximations to the jj ( `)j for 1 j d, giving an approximation to jNK ( ` )j. { prime ideal factorization of G`. { Hermite normal form of H`. { value of s`. For ` = 1, these information are obtained from the initial values of the variables. Each step ` consists of: 1. Select an integral ideal I` of almost xed norm, by multiplying H` with another integral ideal dividing the numerator (resp. the denominator) of G` if s` = 1 (resp. s` = ?1). Compute its Hermite normal form. 2. Pick some \nice" ` in I` using lattice reductions. 3. De ne:
`+1 = ` `?2s` ; G`+1 =
I ?s` ` H`
G` ; H`+1 = ; s`+1 = ?s` : `
This allows to easily update necessary information: { compute the jj (` )j's to approximate the jj ( `+1 )j's. { the selection of I` is actually made in order to obtain the prime ideal factorization of G`+1 simply by updating the exponents of the prime ideal factorization of G` . { H`+1 and s`+1 are directly computed. 4. Store s` and the integral representation of ` . We now explain the meaning of the di erent then we detail the rst hQ`?1 svariables, i2 Q L two parts. By induction on `, = ` L=1 L : In other words, `L?=11 LsL is the approximation of p at step `. Each ` is a square and G` = H`s` p < ` > : Notice that C (G`+1 ) = N (I`1=H` ) C (G` ):
10
Ideal selection. We try to select an I` with norm as close as possible to a
constant LLLmax, set at the beginning of the iterative process, to be explained later on. To do so, we adopt a greedy strategy. Since we know the prime ideal factorization of G` , we can sort all the prime ideals (according to their norm) appearing in this factorization. We start with I` = H` , and we keep multiplying I` by the possibly largest prime ideal power in such manner that N (I` ) is less than LLLmax. In practice, this strategy behaves well because most of our prime ideals lie over small primes. At the same time, when we pick a prime ideal power to multiply with I` , we update its exponent in the prime ideal factorization of G` so that we obtain the prime ideal factorization of G`+1 . At the end of the approximation, when C (G` ) is small, we nd an I` of small norm (not close to LLLmax) such that HI`` equals the whole numerator or the whole denominator of G` .
Integer selection. We look for a nice element ` in the integral ideal I` , that is to say, an algebraic integer that looks like the ideal. For us, \looking like" will mainly mean \with norm almost alike". This really means something since the norm of any element is a multiple of the norm of the integral ideal. So we select ` in order to make N ( < ` > =I` ) as small as possible, which is the same as nding a short element in a given ideal. Fortunately an ideal is also a lattice, and there exists a famous polynomial-time algorithm for lattice reduction: LLL [9, 4]. We will use two features of the LLL-algorithm: computation of an LLL-reduced basis, and computation of a short vector (with respect to the Euclidean norm, not to the norm in a number eld). First, we reduce the basis of I` given by its HNF. In other words, we reduce the matrix of the integral representations (with respect to (!1 ; : : : ; !d )) of the elements of the basis. We do so because the HNF matrix is triangular, therefore not well-balanced: by applying an LLL reduction, coecients are smaller and better spread. Assume the obtained reduced basis is (v(j) )dj=1 . We specify a constant c > 0 by s LLL K ( ` )js` max cd = N (I ) jNj (K )j : ` Let j = jj ( `c)js`=2 for 1 j d. We de ne a linear transformation that maps P any v = di=1 vi !i 2 I` to v = [v1 ; : : : ; vd ; 1 1 (v); : : : ; d d (v)]t : This is when K is totally real. If f has complex roots: for any complex conjugate pairs i and p i , we replacepi (v) and i (v) in the de nition of by respectively, is huge, we use a hash-table and represent any normal prime ideal by its corresponding (p; r) pair. Exceptional prime ideals require more place, but there are very few exceptional primes. 2. It is only during the approximation process (namely, to obtain the Hermite normal form of I` ) that one needs to compute a system of O-generators for normal prime ideals. Such a computation is however very fast.
14 3. To avoid over ows, we do not compute jj ( ` )j, c and j but their logarithms. P One checks that dj=1 ln jj ( ` )j = ln jNK ( ` )j if one is in doubt about the precision. 4. To choose the constant LLLmax, one can compute the C constant from the formulas given in the proof of Theorem 3, but one can also perform some LLL reductions to obtain the practical value of C . Notice that when one knows C and LLLmax, one can estimate the number of iterations. 5. To know how many good primes are sucient to compute the last algebraic integer, one can compute the C 0 constant as shown in the proof of Proposition 4, which gives a bound for the coecients of the integral representation. 6. The last algebraic integer is often a small root of unity. This is because the last ideal I` is principal, and we know an approximation to the embeddings of one of its generators. This generator has unusual short norm in the corresponding lattice, therefore it is no surprise that the LLL algorithm nds this generator, making H`+1 equal to < 1 > . In the latter case, the last algebraic integer is often equal to 1: one should try to bypass the computation of the error and apply directly to nd some factors of n. The algorithm has been implemented using version 1.39 of the PARI library [1] developed by Henri Cohen et al. In December, 1996, it completed the factorization of the 100-digit cofactor of 17186 + 1, using the quadratic polynomials 5633687910X 2 ? 4024812630168572920172347X +482977515620225815833203056197828591062 and ?77869128383X 2 ? 2888634446047190834964717X + 346636133525639208946167278118238554489. Each dependency had about 1.5 million relations. It took the square root code about 10 hours to do both square roots on a 75Mhz Sparc 20.
7 Conclusion We presented an algorithm suitable for implementation to solve the square root problem of the number eld sieve. This algorithm is a variant of Montgomery's square root. We modi ed the square root approximation process by using an integral basis instead of the power basis: this allows to work with integers instead of rationals, and to search the algebraic integer ` in the whole ideal I` , not in some of its submodules. We introduced the simulated annealing method in the ideal simpli cation process. From results of [3], we proposed an ecient ideal square root process and proved its validity. We postponed the computation of the error to avoid useless computations. The present running time of the algorithm is negligible compared to other stages of the number eld sieve. In practice, the algorithm behaves as if it had linear complexity, but one should note that this is only heuristic as few things are proved about the complexity. It is an open problem to determine precisely the complexity of the algorithm. Acknowledgements. I am particularly grateful to both Arjen and Hendrik Lenstra for many explanations about the number eld sieve. I wish to thank Jean-Marc Couveignes and Peter Montgomery for enlightening discussions. I also thank Philippe Hoogvorst for his helpful comments, and for carrying out experiments.
15
A Proof of Theorem 3 This theorem is related to the classical result of the geometry of numbers which states that for any integral ideal I , there exists an algebraic integer 2 I such that jNK ()j M(K )N (I ) where M(K ) denotes the Minkowski constant of K . It relies on Minkowski's convex body theorem which can be viewed as a generalization of the pigeon-hole principle. Following an idea of Montgomery [14], we use the pigeon-hole principle to estimate precisely each component of ` . The only thing we need to know about LLL-reduced bases is that if (b1 ; : : : ; bd ) is an LLL-reduced basis of a lattice , then det()
Yd
kbi k 2d d? (
= det()
(1)
1) 4
i=1 kb1 k 2(d?1)=2kxk if x 2 ; x 6= 0
(2)
where det denotes the lattice determinant and k:k denotes the Euclidean norm. In the following, we will use the notation k:k even for vectors with di erent P d numberqof coordinates. Here, if x = i=1 xi !i is an algebraic number of K , then kxk = Pdi=1 x2i . We will use the notation (x)i to denote the i-th coordinate of x. From now on (all along the proof), we assume that K is totally real to simplify the de nition of , but a similar reasoning applies to other cases with a di erent choice of constants.
Lemma 5. There exists a computable constant C depending only on K such that for every x 2 K , and for any integer j = 1; : : : ; d: 1
jj (x)j C kxk j( x)d j j j C kxk
(3) (4)
1
+
1
P P Proof. We have x = di=1 xi !i where xi 2 Q . Therefore j (x) = di=1 xi j (!i ). Using triangle inequality and Cauchy-Schwarz, we obtain: jj (x)j
d X i=1
v u d uX jxi jjj (!i )j t jxi j
where C1 = max1jd de nition of .
i=1
2
v u d uX t jj (!i )j kxkC ; 2
1
i=1
qPd i jj (!i )j . This proves (3), which implies (4) by =1
2
ut
Lemma 6. There exists two computable constants C and C depending only 2
3
on K such that for any integral ideal I` , there exists a real M and an algebraic
16 integer z 2 I` ; z 6= 0 satisfying:
M d C2
Y j 2J
j
(5)
kzk M N (I` ) =d 8j 2 J j kzk M N (I` ) =d k z k C M N (I` ) =d
(6) (7) (8)
1
1
1
3
where J = fj = 1; : : : ; d = j > 1g. Proof. Let C2 = 2d(d?1)=4dd 2d+1. Since 2d(d?1)=4dd
Y Y dj e < C j by de j 2J j 2J Y dY d 2
nition of J , there exists M > 0 such that 2d(d?1)=4d
dj e < M C2 j : j 2J j 2J = (n1 ; : : : ; nd ) 2 N d such that each ni
This M satis es (5). The number of n satis es ni kv(i) k Md N (I` )1=d is at least
Yd M N (I`) =d Yd M N (I` ) =d d e 1
1
i=1
dkv i
( )
k
dkv i
( )
i=1
k
M d by (1) > Y d e: j d d d 2 (d?1)=4 j 2J
(i)
k c is a positive integer less than j . By the pigeonFor such an n, bj MnNi dk(Iv` )1=d hole principle, there therefore exists two distinct n = (n1 ; : : : ; nd) and n0 = (n01 ; : : : ; n0d) both in N d such that for all i = 1; : : : ; d:
P
ni kv(i) k Md N (I` )1=d n0i kv(i) k Md N (I` )1=d (i) 0 (i) 8j 2 J bj MniNdk(vI )1k=d c = bj MniNdk(vI )1k=d c `
`
(9) (10) (11)
De ne z = di=1 (ni ? n0i )v(i) . Then z 2 I` ; z 6= 0 and by (9) and (10), we have for all i = 1; : : : ; d:
jni ? n0i j:kv i k Md N (I` ) =d : ( )
1
This proves (6) by triangle inequality . Furthermore, for all j 2 J and for all i = 1; : : : ; d, the quantity j jni ? n0i j:kv(i) k is equal to
M N (I )1=d ni dkv(i) k ? n0i dkv(i) k ; d ` j M N (I` )1=d j M N (I` )1=d
17 which is, by (11), less than Md N (I` )1=d . This proves (7) by triangle inequality. Finally:
k z k = 2
d X j =1
j( z )j j + 2
kzk + 2
X
d X j =1
j( z )d j j +
j C1 kzk2 +
X
2
j C1 kzk2 by (4)
j 2J 0 j62J 1 h i X X @1 + C 1 + C 1A M N (I` ) =d 1
1
j 62J
1
j 2J
2
p
by (6), (7) and the de nition of J . This proves (8) with C3 = 1 + dC1 .
ut
Now, if is the algebraic integer output by the second LLL reduction, (2) implies that k k2 2d?1k z k2. Since kk k k, (8) implies that
kk 2 d? = C M N (I` ) =d : (
1) 2
1
3
Q Q Q Moreover, jNK ()j = dj=1 jj ()j = j2J jj ()j) j62J jj ()j : On the one hand, by (3): Y
j 62J
h
jj ()j (C kk)d?jJ j 2 d? = C C M N (I` ) =d 1
(
1) 2
1
1
3
id?jJ j
:
Qj2J j d+j j Qj2J j , where by the arithmetic-
Q
On the other hand, j2J jj ()j = geometric mean inequality:
(
)
0 1jJ j X j( )d j j @ j( )d j j A (k k )jJ j (2d? k z k )jJ j j 2J j 2J i h d? = =d jJ j Y
+
2
+
2
(
1) 2
2
2
C3 M N (I` )1
We collect these two inequalities: d?jJ j
h
1
by (8).
id?jJ j
jNK ()j QC 2 d? = C M N (I` ) =d j 2J j d Q ; C ) 2d d? = C dM dN (I` ) max(1 1
1
(
1) 2
(
3
1
jJ j
+
1) 2
3 j 2J j max(1; C1d )2d(d?1)=2C3d C2 N (I` ) by (5).
This completes the proof with C = 2d(d?1)=2 max(1; C1d )C2 C3d :
2
18
References 1. Batut, C., Bernardi, D., Cohen, H., and Olivier, M. Pari-gp computer package. Can be obtained by ftp at megrez.math.u-bordeaux.fr. 2. Buchmann, J. A., and Lenstra, Jr., H. W. Approximating rings of integers in number elds. J. Theor. Nombres Bordeaux 6, 2 (1994), 221{260. 3. Buhler, J. P., Lenstra, H. W., and Pomerance, C. Factoring integers with the number eld sieve. pages 50-94 in [8]. 4. Cohen, H. A course in computational algebraic number theory. Springer, 1993. 5. Couveignes, J.-M. Computing a square root for the number eld sieve. pages 95-102 in [8]. 6. Cowie, J., Dodson, B., Elkenbracht-Huizing, R. M., Lenstra, A. K., Montgomery, P. L., and Zayer, J. A world wide number eld sieve factoring record: On to 512 bits. In Proceedings of ASIACRYPT'96 (1996), vol. 1163 of Lecture Notes in Computer Science, Springer-Verlag, pp. 382{394. 7. Elkenbracht-Huizing, M. An implementation of the number eld sieve. Experimental Mathematics 5, 3 (1996), 231{253. 8. Lenstra, A. K., and Lenstra, Jr., H. W. The development of the Number Field Sieve, vol. 1554 of Lecture Notes in Mathematics. Springer-Verlag, 1993. 9. Lenstra, A. K., Lenstra, Jr., H. W., and Lovasz, L. Factoring polynomials with rational coecients. Math. Ann. 261 (1982), 515{534. 10. Lenstra, A. K., Lenstra, Jr., H. W., Manasse, M. S., and Pollard, J. M. The number eld sieve. pages 11-42 in [8]. 11. Lenstra, A. K., Lenstra, Jr., H. W., Manasse, M. S., and Pollard, J. M. The factorization of the ninth fermat number. Math. Comp. 61 (1993), 319{349. 12. Lenstra, Jr., H. W. Factoring integers with elliptic curves. Ann. of Math. 126 (1987), 649{673. 13. Lenstra, Jr., H. W. Algorithms in algebraic number theory. Bull. Amer. Math. Soc. 26 (1992), 211{244. 14. Montgomery, P. L. Square roots of products of algebraic numbers. Draft of June, 1995. Available at ftp://ftp.cwi.nl/pub/pmontgom/sqrt.ps.gz. 15. Montgomery, P. L. Square roots of products of algebraic numbers. In Mathematics of Computation 1943-1993: a Half-Century of Computational Mathematics (1994), W. Gautschi, Ed., Proceedings of Symposia in Applied Mathematics, American Mathematical Society, pp. 567{571. 16. Pohst, M., and Zassenhaus, H. Algorithmic algebraic number theory. Cambridge University Press, 1989. 17. Pollard, J. M. Factoring with cubic integers. pages 4-11 in [8]. 18. Reeves, C. R. Modern Heuristic Techniques for Combinatorial Problems. Blackwell Scienti c Publications, 1993.
c Springer-Verlag (http://www.springer.de/comp/lncs/index.html)
A Montgomery-like Square Root for the Number Field Sieve Phong Nguyen [email protected] http://www.dmi.ens.fr/~pnguyen/
Ecole Normale Superieure Laboratoire d'Informatique 45, rue d'Ulm F - 75230 Paris Cedex 05
Abstract. The Number Field Sieve (NFS) is the asymptotically fastest
factoring algorithm known. It had spectacular successes in factoring numbers of a special form. Then the method was adapted for general numbers, and recently applied to the RSA-130 number [6], setting a new world record in factorization. The NFS has undergone several modi cations since its appearance. One of these modi cations concerns the last stage: the computation of the square root of a huge algebraic number given as a product of hundreds of thousands of small ones. This problem was not satisfactorily solved until the appearance of an algorithm by Peter Montgomery. Unfortunately, Montgomery only published a preliminary version of his algorithm [15], while a description of his own implementation can be found in [7]. In this paper, we present a variant of the algorithm, compare it with the original algorithm, and discuss its complexity.
1 Introduction The number eld sieve [8] is the most powerful known factoring method. It was rst introduced in 1988 by John Pollard [17] to factor numbers of form x3 + k. Then it was modi ed to handle numbers of the form re ? s for small positive r and jsj: this was successfully applied to the Fermat number F9 = 2512 + 1 (see [11]). This version of the algorithm is now called the special number eld sieve (SNFS) [10], in contrast with the general number eld sieve (GNFS) [3] which integers. n in heuristic time can handle 1arbitrary GNFS factors1=integers = 3 2=3 3 exp (cg + o(1)) ln n ln ln n with cg = (64=9) 1:9. Let n be the composite integer we wish to factor. We assume that n is not a prime power. Let Zn denote the ring Z=nZ. Like many factoring algorithms, the number eld sieve attempts to nd pairs (x; y) 2 Z2n such that x2 y2 (mod n). For such a pair, gcd(x ? y; n) is a nontrivial factor of n with a probability of P d 1 j at least 2 . The NFS rst selects a primitive polynomial f (X ) = j=0 cj X 2 Z[X ] irreducible over Z, and an integer m with f (m) 0 (mod n). Denote by F (X; Y ) = Y df (X=Y ) in Z[X; Y ] the homogenous form of f . Let 2 C be a
2 root of f , and K = Q ( ) be the corresponding number eld. There is a natural ring homomorphism from Z[ ] to Zn induced by ( ) m (mod n). We will do as if mapped the whole K . If ever ( ) is not de ned for some 2 K , then we have found an integer not invertible in Zn, and thus, a factor N of n which should not be trivial. If n0 = n=N is prime, the factorization is over, and if not, we replace n by n0 , and by 0 induced by 0 ( ) m (mod n0 ). By means of sieving, the several integer pairs (ai ; bi ) and a nite Q NFS nds ? bi ) and Qi2S (aiQ? bi m) are squares in K nonempty set S such that i2S (ai?Q and in Z, respectively. We have i2S (ai ? bi ) i2S (ai ? bi m) (mod n); therefore
13 2s 3 2 0s Y Y 4 @ (ai ? bi )A5 4 (ai ? bi m)5 2
i2S
2
i2S
(mod n)
after extracting the square roots, which gives rise to a suitable pair (x; y). The NFS does root of the Q not specify how to evaluate these square roots. The square prime factorizainteger i2S (ai ? bi m) mod n can be found using the known Q tions of each ai ? bi m. But extracting the square root of i2S (aiQ? bi ) is much more complicated and is the subject of this paper. We note = i2S (ai ? bi ): The following facts should be stressed: { the cardinality jS j is large, roughly equal to the square root of the run time of the number eld sieve. It is over 106 for n larger than 100 digits. { the integers ai; bi are coprime, and t in a computer word. { the prime factorization of each F (ai ; bi) is known. { for every prime number p dividing cd or some F (ai; bi ), we know the set R(p) consisting of roots of f modulo p, together with 1 if p divides cd . The remainder of the paper is organized as follows. In Section 2, we review former methods to solve the square root problem, one of these is used in the last stage of the algorithm. Section 3 presents a few de nitions and results. In Section 4, we describe the square root algorithm, which is a variant of Montgomery's original algorithm, and point out their di erences and similarities. We discuss its complexity in Section 5. Finally, we make some remarks about the implementation in Section 6, and the appendix includes the missing proofs.
2 Former methods UFD method. If is an algebraic integer and the ring Z[ ] is a unique factorization domain (UFD), then each ai ? bi can be factored into primes and units,
and so can be , which allows us to extract a square root of . Unfortunately, the ring Z[ ] is not necessarily a UFD for the arbitrary number elds GNFS encounters. And even though Z[ ] is a UFD, computing a system of fundamental units is not an obvious task (see [4]). The method was nevertheless applied with success to the factorization of F9 [11].
3
Brute-force method. One factorizes the polynomial P (X ) = X ? over
K [X ].
2
To do so, one has to explicitly write the algebraic number , for instance by expanding the product: one thus gets the (rational) coecients of as a polynomial of degree at most d ? 1 in . But there are two serious obstructions: the coecients that one keeps track of during the development of the product have O(jS j) digits. Hence, the single computation of the coecients of can dominate the cost of the whole NFS. And even if we are able to compute , it remains to factorize P (X ). One can overcome the rst obstruction by working with integers instead of rationals: let f^(X ) be the monic polynomial F (X; cd ), and ^ be the algebraic integer cd which is a root of f^. If is a square in K then 0 = c2ddjSj=2e f^0 (^ )2 is a square in Z[^ ], where f^0 denotes the formal derivative of f^. It has integral coecients as a polynomial of degree at most d?1 in ^, and these can be obtained with the Chinese Remainder Theorem, using several inert primes (that is, f is irreducible modulo this prime) if there exist inert primes (which is generally true). This avoids computations with very large numbers. However, one still has to factorize the polynomial Q(X ) = X 2 ? 0 ; whose coecients remain huge, so the second obstruction holds. Furthermore, a large number of primes is required for the Chinese Remainder Theorem, due to the size of the coecients.
Couveignes's method. This method overcomes the second obstruction. If f has odd degree d, Couveignes [5] remarks that one is able to distinguish the two p square roots of any square in K , by specifying its norm. Let 0 be the square root with positive norm. Since the prime factorization of N ( 0 ) is known, the p 0 integer any prime q. If q is inert p N ( ) can be eciently computed modulo (mod q). From the Chinese then 0 (mod q) can be computed after expanding 0 p Remainder Theorem, one recovers the coecients of 0 2 Z[^ ]. One can show that the complexity of the algorithm is at best O(M (jS j) ln jS j), where M (jS j) is the time required to multiply two jS j-bit integers. The algorithm appears to
be impractical for the sets S now in use, and it requires an odd degree. Montgomery's strategy [15, 14, 7] can be viewed as a mix of UFD and bruteforce methods. It bears some resemblance to the square root algorithm sketched in [3] (pages 75-76). It works for all values of d, and does not make any particular assumption (apart from the existence of inert primes) about the number eld.
3 Algebraic preliminaries Our number eld is K = Q ( ) = Q (^ ), where is an algebraic number and ^ = cd is an algebraic integer. Let O be its ring of integers, and I be the abelian group of fractional ideals of O. For x1 ; : : : ; xm 2 K , we note < x1 ; : : : ; xm > the element of I generated by x1 ; : : : ; xm . For every prime ideal p, we denote by vp the p-adic valuation that maps I to Z. We de ne the Q numerator and denominator of I 2 I to be the integral ideals numer(I ) = vp (I )>0 pvp (I ) and Q denom(I ) = vp (I ) in I . We follow the notations of [3] and recall some results. Let R be an order in O. By a \prime of R" we mean a non-zero prime ideal of R. We denote by flp;R : K ! Zgp the unique collection (where p ranges over the set of all primes of R) of group homomorphisms such that: { lp;R(x) 0 for all x 2 R; x 6= 0; { if x is a non-zero element of R, then lp;R(x) > 0 if and only if x 2 p; { for Y eachlpx;R(2x) K one has lp;R(x) = 0 for all but nitely many p, and N (p) = jNK (x)j; where p ranges over the set of all primes of R. p
lp;O (x) coincide with vpP ( < x > ). Let i = cd d?1?i + cd?1 d?2?i + + ci+1 . ?2 i Z is an order of O, which is in fact Z[ ] \ Z[ ?1]. We know that A = Z + id=0 Its discriminant (A) is equal to (f ) and we have: (d?1)(d?2) 2 [
(Z[^ ]) = c(dd?1)(d?2)(A);
[O : Z[^ ]] = cd
O : A]:
Recall that for any prime number p, R(p) is de ned as the set consisting of roots of f modulo p, together with 1 if p divides cd. Note that this R(p) is denoted R0 (p) in [3]. The pairs consisting of a prime number p and an element r 2 R(p) are in bijective correspondence with the rst degree primes p of A: { if r 6= 1 then p is the intersection of A and the kernel of the ring homomorphism p;r : Z[ ] ! Fp that sends to r. { if r = 1 then p is the intersection of A and the kernel of the ring homomorphism p;1 : Z[ ?1] ! Fp that sends ?1 to 0. Let p be a prime number, r an element of R(p) and a, b be coprime integers. If a br (mod p) and r 6= 1, or if b 0 (mod p) and r = 1, we de ne ep;r (a; b) = vp (F (a; b)) where vp denotes the ordinary pQ-adic valuation. Otherwise, we set ep;r (a; b) = 0. We have NK (a ? b ) = c1d p;r pep;r (a;b) ; the product ranging over all pairs p, r with p prime and r 2 R(p). Furthermore, for any coprime integers a, b and any rst degree prime p of A corresponding to a pair p, r 2 R(p), we have: e (a; b) 1 lp;A (a ? b ) = ep;r (a; b) ? v (c ) ifif rr 6= =1 p;r
p d
Theorem 1. Let a and b be coprime integers, and p be a prime number. Let p be a prime ideal of O above p such that vp ( ) 6= 0. If p does not divide [O : A] then:
5 1. For every r 2 R(p), there is a unique prime ideal pr of O that lies over the rst degree prime ideal qr of A corresponding to the pair p, r. pr is a rst degree prime ideal, given by pr = : Furthermore, we have vpr ( ) = lqr ;A (a ? b ). 2. There is at most one nite r 2 R(p) such that ep;r (a; b) 6= 0: 3. If p does not divide cd, such a nite r exists and p = pr . 4. If p divides cd, then either p is p1 , or pr for r nite. 5. p divides F (a; b) or cd . Proof. Let r 2 R(p) and qr be the rst degree prime ideal of A corresponding to the pair p,r. Since P p does not divide [O : A], we have from [3] (Proposition 7.3, pages 65-66): pr jqr f (pr =qr ) = 1; where pr ranges over all primes of O lying over qr and f denotes the residual degree. This proves that pr is unique and is a rst degree prime ideal. From [3] (Proposition 7.2, page 65), we also have:
lqr ;A (a ? b ) =
X
p0 jqr
f (p0 =qr )lp0 ;O (a ? b ) = lpr ;O (a ? b ):
Hence, vpr (a ? b ) = lqr ;A (a ? b ). Moreover, we know a Z-basis for any ideal qr of A, namely (p; 0 ? p;r ( 0 ); : : : ; d?2 ? p;r ( d?2 )): Since pr lies over qr , this Z-basis is a system of O-generators for pr . We therefore proved 1. From the de nition of i , one sees that i = ci ?1 + ci?1 ?2 + + c0 ?i?1 ; which proves that p;1 ( i ) = 0. This simpli es the formula when r = 1. One obtains 2 from the de nition of ep;r . Denote by q the intersection of p and A. q is a prime of A and p lies over q. We have lq;A (a ? b ) 6= 0 since vp (a ? b ) 6= 0. From [3] (page 89), this proves that q is a rst degree prime ideal of A. Hence, there exists r 2 R(p) such that q = qr . From 1, this proves that p = pr . This r is nite or in nite, and if r is nite, it is the r of 2. This proves 3 and 4. From the formula expressing lq;A (a ? b ) in terms of ep;r (a; b), we obtain 5. ut
4 The square root algorithm We that we want to compute a square root of the algebraic number = Q recall ( a ? i2S i bi ): The algorithm is split as follows:
1. Transform in order to make < > simpler. The running time of the rest of the algorithm heuristically depends on C ( < > ). p 2. Compute < > from the prime ideal factorization of < > given by the prime factorization of each F (ai ; bi ). p p 3. Approximate from < > : using lattice reductions, construct a sequence of algebraic Q integers 1; : : : ; L in O and signs s1; : : : ; sL in f1g such that = L`=1 `?2s` is a \small" algebraic integer. Qcan be thought as the square of the \guessing-error" inpthe approximation L`=1 `s` of p . 4. Since is a square, so is . Compute using brute-force method. One is able to explicitly write because is a \small" algebraic integer.
6
We thus obtain p as a product of algebraic integers with exponents 1: L p = p Y s` : `=1
`
This enables to compute (p ) without explicitly calculating p , and hopefully some factors of n. Although formalized di erently, Montgomery's algorithm uses the same strategy. Only the steps change. We use another heuristic approach in Step 1, which seems to be more e ective in practice. We use a new process in Step 2, derived from Section 3. Montgomery used a process which was as ecient, but only heuristic. Step 3 is the core of the algorithm. We modi ed this step by using the integral basis in a systematic manner, instead of the power basis. This simpli es the algorithm and the proofs. Heuristically, this should also improve the performances. We postpone the computation of the error in Step 4, while Montgomery included it in Step 3, by updating the computations during the approximation. This decreases the running-time because it is easier to estimate the necessary computations when Step 3 is over, and sometimes, Step 4 can be avoided (when the approximation is already perfect, which can be checked without additional computations). The new algorithm might be more suited to analysis, but like Montgomery's algorithm, its complexity has yet to be determined, even though they both display signi cantly better performances than former methods.
4.1 Computing in the number eld The ring of integers. During the whole algorithm, we need to work with ideals and algebraic integers. We rst have to compute an integral basis of O.
In general, this is a hopeless task (see [13, 2] for a survey), but for the number elds NFS encounters (small degree and large discriminant), this can be done by the so-called round algorithms [16, 4]. Given an order R and several primes pi , any round algorithm will enlarge this order for all these primes so that the new order Rb is pi -maximal for every pi . If we take for the pi all the primes p such that p2 divides (R), then Rb = O. To determine all these primes, a partial factorization of (R) suces, that is a factorization of the form df 2 where d is squarefree and f is factorized. Theoretically, a partial factorization is as hard to nd as a complete factorization and unfortunately, the discriminant is sometimes much larger than the number n we wish to factor. However, if one takes a \random" large number, and one removes all \small" prime factors from it (by trial division or by elliptic curves [12]), then in practice the result is quite likely to be squarefree. Furthermore, even in the case Rb 6= O, it will be true that Rb has almost all of the good properties of O for all ideals that we are likely to encounter in practice, like the fact that every ideal is a product of prime ideals. This is because every order satis es these properties for all ideals that are coprime to the index of the order in O. Hence, we can now assume that an integral basis (!1 ; : : : ; !d) of O has been computed.
7
Algebraic numbers and ideals. From this integral basis we can represent any algebraic number of K as a vector of Q d : thisPis the integral representation. If x 2 K we de ne x = [x1 ; : : : ; xd ]t where x = di=1 xi !i and xi 2 Q . We can also represent any algebraic number as a polynomial of degree at most d ? 1 in : this is the power representation. When dealing with algebraic integers, the integral representation is preferable. We will represent any integral ideal I by an integral matrix (with respect to (!1 ; : : : ; !d)) from a Z-basis or a system of O-generators. In the case of Z-basis, we use the Hermite normal form (HNF) of the square matrix for eciency reasons. We refer to [4] for algorithms concerning algebraic numbers and ideals. 4.2 Simplifying the principal ideal Q If is p a square in K , then so is any 0 = i2pS (ai ? bi p)ei ; when ei = 1. Since Q p (ai ? bi ); we can recover from 0 but actually, we only
= 0 ei =?1
look for a square identity. Fortunately:
2 s 3 3 2s Y Y 4 ( (ai ? bi )ei )5 4 (ai ? bim)ei 5 2
i2S
2
(mod n)
i2S
p
This replaces the computation of p by the computation of 0 . By cleverly selecting the ei , C ( < 0 > ) will be much smaller than C ( < > ): this is because many share the same prime ideals, since many NK (ai ? bi ) share the same primes (as a consequence of sieving). We now address the optimization problem of selecting the ei so that C ( < 0 > ) is small. Given a distribution of ei , the complexity of < 0 > can be computed by the following formula (which comes from the known \factorization" of each ai ? bi into primes of A):
Y jPi2S eiep;r ai;bi j Y jPi2S ei ep;1 ai;bi ?vp cd j p p : (
p;r6=1
)
[
pjcd
(
)
(
)]
The simplest method is a random strategy which selects randomly ei = 1. Another method is a greedy strategy (used in [7]): at every step, select ei = 1 according to the best complexity (whether we put ai ? bi in the numerator or in the denominator). This behaves better than the random strategy. But the best method so far in practice is based on simulated annealing [18], a well-known probabilistic solution method in the eld of combinatorial optimization. Here, the con guration space is E = f?1; +1gjSj, and the energy function U maps any e = (e1 ; : : : ; ejS j ) 2 E to ln C ( < > ) where corresponds to e. For any e 2 E , we de ne its neighbourhood V (e) = f(e1; : : : ; ei?1 ; ?ei ; ei+1 ; : : : ; ejSj) j i = 1; : : : ; jS jg. We try to minimize U by the following algorithm, which performances depend on three parameters i ; f (initial and nal temperatures) and : { select randomly e 2 E and set ? i . { choose randomly f 2 V (e) and set ? U (f ) ? U (e). If > 0, set p ? exp(?=), otherwise set p ? 1. Then set e ? f with probability p, and ? .
8
{ repeat previous step if > f . Although this method behaves better in practice than previous methods, theoretical estimates can hardly be given.
4.3 Ideal square root
Q
From now on, we forget about the initial and set = i2S (ai ? bi )ei : p We wish to obtain as a product of ideals with exponents lying in Z (this ideal is too large to be represented as a single matrix). This Q can be done by factoring into prime ideals the fractional ideal < > = < i2S (ai ? bi )ei > . We simplify the problem to the factorization of any linear expression < ai ? bi > with coprime ai ; bi . Such a factorization could be obtained by general ideal factorization algorithms (see [4]) but this would be too slow if we had to use these algorithms jS j times. Fortunately, we can do much of the work by ourself using the known factorization of each F (ai ; bi ) = f (ai =bi )bdi , as shown in the previous section. We say that a prime number p is exceptional if p divides the index = [O : A]. Otherwise, we say that p is normal. Naturally, a prime ideal of O is said to be exceptional (resp. normal) if it lies above an exceptional (resp. normal) prime. If m is the number of prime factors of , there are at most md exceptional prime ideals. We compute all the exceptional prime ideals (for example, by decomposing all the exceptional primes in O using the BuchmannLenstra algorithm described in [4]), along with some constants allowing us to compute eciently any valuation at these primes. From Theorem 1, we get the prime ideal factorization of < a ? b > as follows: for every prime number p dividing cd or such that there exists a nite r 2 R(p) satisfying ep;r (a; b) 6= 0, { if p is exceptional, compute the valuation of at all the exceptional ideals lying above p. { otherwise, p is normal. If there is a nite r 2 R(p) such that ep;r (a; b) 6= 0 (r is then unique), pick the prime ideal pr with exponent ep;r (a; b) where = : If 1 2 R(p), also pick the prime ideal p1 with exponent ep;1 (a; b) ? vp (cd ) where p1 = : We thus decompose < > pas a product of ideals where every exponent is necessarily even, which gives < > . Montgomery used a di erent ideal factorization process (see [7, 14]) by introducing a special ideal, but its correctness is not proved. pr
4.4 Square root approximation We now use the ideal square root p < > to approximate p . Since p < >
is a huge ideal, we will get an approximation through an iterative process, by selecting a small part of the ideal at each step: this small part will be alternatively
9 taken in the numerator and denominator. To lift an integral ideal to an algebraic integer, we use lattice reduction techniques. We associate several variables at each step `: { an algebraic number `. It canp be considered as the square of the error in the current approximation of . { a sign s` in f?1; +1g, indicating whether we take something in the denominator or in the numerator of the huge original ideal. { a fractional ideal G`, which is an approximation to p < ` > . p { an integral ideal H` of bounded norm. It di erentiates G` from < ` > . { an algebraic integer `. { an integral ideal I` of bounded norm. Q We initialize these variables by: 1 = = i2S (ai ? bi )ei , G1 = p < > , H1 = < 1 > , s1 = 1 if NK ( ) 1 and ?1 otherwise. Each step of the approximation makes `+1 in some sense smaller than ` , and G`+1 simpler than G` . After enough steps, G` is reduced to the unit ideal < 1 > , and ` becomes an algebraic integer suciently small that its integral representation can be determined explicitly (using Chinese Remainders) and a square root constructed using brute-force method. At the start of step `, we need to know the following: { approximations to the jj ( `)j for 1 j d, giving an approximation to jNK ( ` )j. { prime ideal factorization of G`. { Hermite normal form of H`. { value of s`. For ` = 1, these information are obtained from the initial values of the variables. Each step ` consists of: 1. Select an integral ideal I` of almost xed norm, by multiplying H` with another integral ideal dividing the numerator (resp. the denominator) of G` if s` = 1 (resp. s` = ?1). Compute its Hermite normal form. 2. Pick some \nice" ` in I` using lattice reductions. 3. De ne:
`+1 = ` `?2s` ; G`+1 =
I ?s` ` H`
G` ; H`+1 = ; s`+1 = ?s` : `
This allows to easily update necessary information: { compute the jj (` )j's to approximate the jj ( `+1 )j's. { the selection of I` is actually made in order to obtain the prime ideal factorization of G`+1 simply by updating the exponents of the prime ideal factorization of G` . { H`+1 and s`+1 are directly computed. 4. Store s` and the integral representation of ` . We now explain the meaning of the di erent then we detail the rst hQ`?1 svariables, i2 Q L two parts. By induction on `, = ` L=1 L : In other words, `L?=11 LsL is the approximation of p at step `. Each ` is a square and G` = H`s` p < ` > : Notice that C (G`+1 ) = N (I`1=H` ) C (G` ):
10
Ideal selection. We try to select an I` with norm as close as possible to a
constant LLLmax, set at the beginning of the iterative process, to be explained later on. To do so, we adopt a greedy strategy. Since we know the prime ideal factorization of G` , we can sort all the prime ideals (according to their norm) appearing in this factorization. We start with I` = H` , and we keep multiplying I` by the possibly largest prime ideal power in such manner that N (I` ) is less than LLLmax. In practice, this strategy behaves well because most of our prime ideals lie over small primes. At the same time, when we pick a prime ideal power to multiply with I` , we update its exponent in the prime ideal factorization of G` so that we obtain the prime ideal factorization of G`+1 . At the end of the approximation, when C (G` ) is small, we nd an I` of small norm (not close to LLLmax) such that HI`` equals the whole numerator or the whole denominator of G` .
Integer selection. We look for a nice element ` in the integral ideal I` , that is to say, an algebraic integer that looks like the ideal. For us, \looking like" will mainly mean \with norm almost alike". This really means something since the norm of any element is a multiple of the norm of the integral ideal. So we select ` in order to make N ( < ` > =I` ) as small as possible, which is the same as nding a short element in a given ideal. Fortunately an ideal is also a lattice, and there exists a famous polynomial-time algorithm for lattice reduction: LLL [9, 4]. We will use two features of the LLL-algorithm: computation of an LLL-reduced basis, and computation of a short vector (with respect to the Euclidean norm, not to the norm in a number eld). First, we reduce the basis of I` given by its HNF. In other words, we reduce the matrix of the integral representations (with respect to (!1 ; : : : ; !d )) of the elements of the basis. We do so because the HNF matrix is triangular, therefore not well-balanced: by applying an LLL reduction, coecients are smaller and better spread. Assume the obtained reduced basis is (v(j) )dj=1 . We specify a constant c > 0 by s LLL K ( ` )js` max cd = N (I ) jNj (K )j : ` Let j = jj ( `c)js`=2 for 1 j d. We de ne a linear transformation that maps P any v = di=1 vi !i 2 I` to v = [v1 ; : : : ; vd ; 1 1 (v); : : : ; d d (v)]t : This is when K is totally real. If f has complex roots: for any complex conjugate pairs i and p i , we replacepi (v) and i (v) in the de nition of by respectively, is huge, we use a hash-table and represent any normal prime ideal by its corresponding (p; r) pair. Exceptional prime ideals require more place, but there are very few exceptional primes. 2. It is only during the approximation process (namely, to obtain the Hermite normal form of I` ) that one needs to compute a system of O-generators for normal prime ideals. Such a computation is however very fast.
14 3. To avoid over ows, we do not compute jj ( ` )j, c and j but their logarithms. P One checks that dj=1 ln jj ( ` )j = ln jNK ( ` )j if one is in doubt about the precision. 4. To choose the constant LLLmax, one can compute the C constant from the formulas given in the proof of Theorem 3, but one can also perform some LLL reductions to obtain the practical value of C . Notice that when one knows C and LLLmax, one can estimate the number of iterations. 5. To know how many good primes are sucient to compute the last algebraic integer, one can compute the C 0 constant as shown in the proof of Proposition 4, which gives a bound for the coecients of the integral representation. 6. The last algebraic integer is often a small root of unity. This is because the last ideal I` is principal, and we know an approximation to the embeddings of one of its generators. This generator has unusual short norm in the corresponding lattice, therefore it is no surprise that the LLL algorithm nds this generator, making H`+1 equal to < 1 > . In the latter case, the last algebraic integer is often equal to 1: one should try to bypass the computation of the error and apply directly to nd some factors of n. The algorithm has been implemented using version 1.39 of the PARI library [1] developed by Henri Cohen et al. In December, 1996, it completed the factorization of the 100-digit cofactor of 17186 + 1, using the quadratic polynomials 5633687910X 2 ? 4024812630168572920172347X +482977515620225815833203056197828591062 and ?77869128383X 2 ? 2888634446047190834964717X + 346636133525639208946167278118238554489. Each dependency had about 1.5 million relations. It took the square root code about 10 hours to do both square roots on a 75Mhz Sparc 20.
7 Conclusion We presented an algorithm suitable for implementation to solve the square root problem of the number eld sieve. This algorithm is a variant of Montgomery's square root. We modi ed the square root approximation process by using an integral basis instead of the power basis: this allows to work with integers instead of rationals, and to search the algebraic integer ` in the whole ideal I` , not in some of its submodules. We introduced the simulated annealing method in the ideal simpli cation process. From results of [3], we proposed an ecient ideal square root process and proved its validity. We postponed the computation of the error to avoid useless computations. The present running time of the algorithm is negligible compared to other stages of the number eld sieve. In practice, the algorithm behaves as if it had linear complexity, but one should note that this is only heuristic as few things are proved about the complexity. It is an open problem to determine precisely the complexity of the algorithm. Acknowledgements. I am particularly grateful to both Arjen and Hendrik Lenstra for many explanations about the number eld sieve. I wish to thank Jean-Marc Couveignes and Peter Montgomery for enlightening discussions. I also thank Philippe Hoogvorst for his helpful comments, and for carrying out experiments.
15
A Proof of Theorem 3 This theorem is related to the classical result of the geometry of numbers which states that for any integral ideal I , there exists an algebraic integer 2 I such that jNK ()j M(K )N (I ) where M(K ) denotes the Minkowski constant of K . It relies on Minkowski's convex body theorem which can be viewed as a generalization of the pigeon-hole principle. Following an idea of Montgomery [14], we use the pigeon-hole principle to estimate precisely each component of ` . The only thing we need to know about LLL-reduced bases is that if (b1 ; : : : ; bd ) is an LLL-reduced basis of a lattice , then det()
Yd
kbi k 2d d? (
= det()
(1)
1) 4
i=1 kb1 k 2(d?1)=2kxk if x 2 ; x 6= 0
(2)
where det denotes the lattice determinant and k:k denotes the Euclidean norm. In the following, we will use the notation k:k even for vectors with di erent P d numberqof coordinates. Here, if x = i=1 xi !i is an algebraic number of K , then kxk = Pdi=1 x2i . We will use the notation (x)i to denote the i-th coordinate of x. From now on (all along the proof), we assume that K is totally real to simplify the de nition of , but a similar reasoning applies to other cases with a di erent choice of constants.
Lemma 5. There exists a computable constant C depending only on K such that for every x 2 K , and for any integer j = 1; : : : ; d: 1
jj (x)j C kxk j( x)d j j j C kxk
(3) (4)
1
+
1
P P Proof. We have x = di=1 xi !i where xi 2 Q . Therefore j (x) = di=1 xi j (!i ). Using triangle inequality and Cauchy-Schwarz, we obtain: jj (x)j
d X i=1
v u d uX jxi jjj (!i )j t jxi j
where C1 = max1jd de nition of .
i=1
2
v u d uX t jj (!i )j kxkC ; 2
1
i=1
qPd i jj (!i )j . This proves (3), which implies (4) by =1
2
ut
Lemma 6. There exists two computable constants C and C depending only 2
3
on K such that for any integral ideal I` , there exists a real M and an algebraic
16 integer z 2 I` ; z 6= 0 satisfying:
M d C2
Y j 2J
j
(5)
kzk M N (I` ) =d 8j 2 J j kzk M N (I` ) =d k z k C M N (I` ) =d
(6) (7) (8)
1
1
1
3
where J = fj = 1; : : : ; d = j > 1g. Proof. Let C2 = 2d(d?1)=4dd 2d+1. Since 2d(d?1)=4dd
Y Y dj e < C j by de j 2J j 2J Y dY d 2
nition of J , there exists M > 0 such that 2d(d?1)=4d
dj e < M C2 j : j 2J j 2J = (n1 ; : : : ; nd ) 2 N d such that each ni
This M satis es (5). The number of n satis es ni kv(i) k Md N (I` )1=d is at least
Yd M N (I`) =d Yd M N (I` ) =d d e 1
1
i=1
dkv i
( )
k
dkv i
( )
i=1
k
M d by (1) > Y d e: j d d d 2 (d?1)=4 j 2J
(i)
k c is a positive integer less than j . By the pigeonFor such an n, bj MnNi dk(Iv` )1=d hole principle, there therefore exists two distinct n = (n1 ; : : : ; nd) and n0 = (n01 ; : : : ; n0d) both in N d such that for all i = 1; : : : ; d:
P
ni kv(i) k Md N (I` )1=d n0i kv(i) k Md N (I` )1=d (i) 0 (i) 8j 2 J bj MniNdk(vI )1k=d c = bj MniNdk(vI )1k=d c `
`
(9) (10) (11)
De ne z = di=1 (ni ? n0i )v(i) . Then z 2 I` ; z 6= 0 and by (9) and (10), we have for all i = 1; : : : ; d:
jni ? n0i j:kv i k Md N (I` ) =d : ( )
1
This proves (6) by triangle inequality . Furthermore, for all j 2 J and for all i = 1; : : : ; d, the quantity j jni ? n0i j:kv(i) k is equal to
M N (I )1=d ni dkv(i) k ? n0i dkv(i) k ; d ` j M N (I` )1=d j M N (I` )1=d
17 which is, by (11), less than Md N (I` )1=d . This proves (7) by triangle inequality. Finally:
k z k = 2
d X j =1
j( z )j j + 2
kzk + 2
X
d X j =1
j( z )d j j +
j C1 kzk2 +
X
2
j C1 kzk2 by (4)
j 2J 0 j62J 1 h i X X @1 + C 1 + C 1A M N (I` ) =d 1
1
j 62J
1
j 2J
2
p
by (6), (7) and the de nition of J . This proves (8) with C3 = 1 + dC1 .
ut
Now, if is the algebraic integer output by the second LLL reduction, (2) implies that k k2 2d?1k z k2. Since kk k k, (8) implies that
kk 2 d? = C M N (I` ) =d : (
1) 2
1
3
Q Q Q Moreover, jNK ()j = dj=1 jj ()j = j2J jj ()j) j62J jj ()j : On the one hand, by (3): Y
j 62J
h
jj ()j (C kk)d?jJ j 2 d? = C C M N (I` ) =d 1
(
1) 2
1
1
3
id?jJ j
:
Qj2J j d+j j Qj2J j , where by the arithmetic-
Q
On the other hand, j2J jj ()j = geometric mean inequality:
(
)
0 1jJ j X j( )d j j @ j( )d j j A (k k )jJ j (2d? k z k )jJ j j 2J j 2J i h d? = =d jJ j Y
+
2
+
2
(
1) 2
2
2
C3 M N (I` )1
We collect these two inequalities: d?jJ j
h
1
by (8).
id?jJ j
jNK ()j QC 2 d? = C M N (I` ) =d j 2J j d Q ; C ) 2d d? = C dM dN (I` ) max(1 1
1
(
1) 2
(
3
1
jJ j
+
1) 2
3 j 2J j max(1; C1d )2d(d?1)=2C3d C2 N (I` ) by (5).
This completes the proof with C = 2d(d?1)=2 max(1; C1d )C2 C3d :
2
18
References 1. Batut, C., Bernardi, D., Cohen, H., and Olivier, M. Pari-gp computer package. Can be obtained by ftp at megrez.math.u-bordeaux.fr. 2. Buchmann, J. A., and Lenstra, Jr., H. W. Approximating rings of integers in number elds. J. Theor. Nombres Bordeaux 6, 2 (1994), 221{260. 3. Buhler, J. P., Lenstra, H. W., and Pomerance, C. Factoring integers with the number eld sieve. pages 50-94 in [8]. 4. Cohen, H. A course in computational algebraic number theory. Springer, 1993. 5. Couveignes, J.-M. Computing a square root for the number eld sieve. pages 95-102 in [8]. 6. Cowie, J., Dodson, B., Elkenbracht-Huizing, R. M., Lenstra, A. K., Montgomery, P. L., and Zayer, J. A world wide number eld sieve factoring record: On to 512 bits. In Proceedings of ASIACRYPT'96 (1996), vol. 1163 of Lecture Notes in Computer Science, Springer-Verlag, pp. 382{394. 7. Elkenbracht-Huizing, M. An implementation of the number eld sieve. Experimental Mathematics 5, 3 (1996), 231{253. 8. Lenstra, A. K., and Lenstra, Jr., H. W. The development of the Number Field Sieve, vol. 1554 of Lecture Notes in Mathematics. Springer-Verlag, 1993. 9. Lenstra, A. K., Lenstra, Jr., H. W., and Lovasz, L. Factoring polynomials with rational coecients. Math. Ann. 261 (1982), 515{534. 10. Lenstra, A. K., Lenstra, Jr., H. W., Manasse, M. S., and Pollard, J. M. The number eld sieve. pages 11-42 in [8]. 11. Lenstra, A. K., Lenstra, Jr., H. W., Manasse, M. S., and Pollard, J. M. The factorization of the ninth fermat number. Math. Comp. 61 (1993), 319{349. 12. Lenstra, Jr., H. W. Factoring integers with elliptic curves. Ann. of Math. 126 (1987), 649{673. 13. Lenstra, Jr., H. W. Algorithms in algebraic number theory. Bull. Amer. Math. Soc. 26 (1992), 211{244. 14. Montgomery, P. L. Square roots of products of algebraic numbers. Draft of June, 1995. Available at ftp://ftp.cwi.nl/pub/pmontgom/sqrt.ps.gz. 15. Montgomery, P. L. Square roots of products of algebraic numbers. In Mathematics of Computation 1943-1993: a Half-Century of Computational Mathematics (1994), W. Gautschi, Ed., Proceedings of Symposia in Applied Mathematics, American Mathematical Society, pp. 567{571. 16. Pohst, M., and Zassenhaus, H. Algorithmic algebraic number theory. Cambridge University Press, 1989. 17. Pollard, J. M. Factoring with cubic integers. pages 4-11 in [8]. 18. Reeves, C. R. Modern Heuristic Techniques for Combinatorial Problems. Blackwell Scienti c Publications, 1993.
