n = p * q 的因數分解之研究

全文

(1)國立交通大學資訊科學與工程研究所碩士論文. n = p. ×. q. 的因數分解之研究. Study on Factorization of n = p × q. 研究生：張煜晧指導教授：葉義雄. 教授. 中華民國九十五年九月.

(2) n = p. ×. q 的因數分解之研究. Study on Factorization of n = p × q. 研究生：張煜晧. Student：Yu-Hao Chang. 指導教授：葉義雄. Advisor：Yi-Shiung Yeh. 國立交通大學資訊科學與工程研究所碩士論文. A Thesis Submitted to Institute of Computer Science and Engineering College of Computer Science National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Master in. Computer Science September 2006 Hsinchu, Taiwan, Republic of China. 中華民國九十五年九月.

(3) n = p. ×. q 的因數分解之研究. 學生：張煜晧. 指導教授：葉義雄博士. 國立交通大學資訊工程學系碩士班. 摘. 要. RSA 密碼系統(RSA Cryptosystem)是使用最為廣泛的公鑰密碼系統之一，其安全性乃建立在大整數難以分解為其質因數乘積的事實之上，此一事實並被稱為 RSA 假定(RSA assumption)。一般相信沒有確定型的圖靈機(deterministic Turing machine，簡稱 DTM)可在多項式時間內破解 RSA 假定，多項式時間的演算法若被發現，RSA 密碼系統將變得不再安全。因為如此，許多科學家致力於研究有效率的分解演算法。目前所知，分解小於 110 位數的大數時，「二次篩選法(quadratic sieve factoring algorithm，簡稱 QS)」是最快的通用演算法。受限於時間與硬體資源，我們主要著眼於 QS 的一種變型，稱之為「複數多項式二次篩選法(multiple polynomial quadratic sieve，簡稱 MPQS)」。為了確認 RSA 假定的強度，我們提出一個方法來加速 MPQS 的篩選程序，其實驗結果將有助於分析 RSA 抵抗現行分解技術的強度，同時可被納入實作 RSA 密碼系統時的考量。. 關鍵字：RSA 密碼系統、因數分解、二次篩選、複數多項式二次篩選法. i.

(4) Study on Factorization of n = p × q. student：Yu-Hao Chang. Advisor：Dr. Yi-Shiung Yeh. Institute of Computer Science and Information Engineering National Chiao Tung University. ABSTRACT. The RSA Cryptosystem is one of the most used public-key cryptosystems. The security it rests on the fact that it is computationally infeasible to factor a large integer into its component primes. This fact is referred to as the RSA assumption. It is believed that there is no deterministic Turing machines (DTM) that can break the RSA assumption in polynomial time. If a polynomial-time algorithm is found, the RSA Cryptosystem would be insecure. Owing to this, many scientists have devoted themselves to researching efficient factoring algorithms. So far, the quadratic sieve factoring algorithm (abbreviated to QS) is the fastest known general-purpose method for factoring numbers having less than about 110 digits. Restricted by time and computer hardware, we focus on one of the variants of the QS, called the multiple polynomial quadratic sieve (MPQS). To ensure the strength of the RSA assumption, we propose a scheme to enhance the sieving procedure of the MPQS. The experimental results are contributive to the analyses of the strength of the RSA assumption against the modern factoring technology and should be taken into consideration on future cryptographic implementations based on the RSA cryptosystem.. Keywords：RSA Cryptosystem, factoring integers, quadratic sieve, multiple polynomial quadratic sieve. ii.

(5) 誌. 謝. 能夠完成這篇論文，首先必須感謝我的指導教授葉義雄老師。這兩年來，承蒙老師悉心指導，使我在研究和待人處世方面獲益匪淺；此外，老師一直容忍我各種任性的要求，如今學業有成，我衷心感謝他。再者，要感謝定宇學長，在我論文寫作的期間，不時給我鼓勵，並與我討論研究方向；這篇論文得以順利付梓，定宇學長居功厥偉。還有特別感謝中山大學的官大智教授，願意提供他珍貴的程式，供我們研究。也感謝韓禹學長，幫忙整理我的論文，不吝指導我寫作的技巧。感謝實驗室的同窗好友，宗哥、鴻祥、阿甘、昇哥，研究所的兩年有幸與你們同行，讓我習得了不少特殊技能，這將是我一生寶貴的資產。碩二的學弟妹，Qting、 Gobby、伯昕，雖然我大你們一屆，卻總是受你們照顧，在此謝謝你們。還有新進的學弟妹，佳君、小強、家明，與你們相處的時日雖短，但非常開心，祝你們接下來的兩年一帆風順。也謝謝以德學長、鎮宇學長、靜紋學姊，一直以來對我照顧有加。最後，要謝謝我的爺爺、奶奶、爸爸、媽媽，多虧你們苦心栽培，今日我才能拿到這個學位，謝謝你們。僅以此文獻給所有支持和關心我的人，謝謝!. 張煜晧中華民國九十五年九月 iii.

(6) Contents 中文摘要........................................................................................................................1 Abstract.........................................................................................................................2 誌謝.............................................................................................................................. iii Contents .......................................................................................................................iv List of Tables.................................................................................................................v Chapter 1 Introduction..........................................................................................1 1.1 Elementary Number Theory.......................................................................1 1.2 The RSA Cryptosystem..............................................................................2 1.3 RSA and Factoring Integers .......................................................................4 Chapter 2 Factoring Algorithms...........................................................................5 2.1 The Dixon’s Random Squares Algorithm ..................................................5 2.2 The Quadratic Sieve Factoring Algorithm.................................................9 2.2.1 Setting Up the Factor Base ..........................................................10 2.2.2 The Sieving Procedure.................................................................12 2.2.3 Improvements on the QS .............................................................14 2.3 The Multiple Polynomial Quadratic Sieve ..............................................15 2.3.1 Polynomials Selection..................................................................16 2.3.2 The Details of Choosing the Coefficients ....................................18 2.3.3 Sieving .........................................................................................21 Chapter 3 3.1 3.2 3.3 3.4. The Modified Multiple Polynomial Quadratic Sieve......................23 Motivation for the Modified Multiple Polynomial Quadratic Sieve .......24 Square Roots of n Modulo pok..................................................................25 Modified Sieving Procedure ....................................................................28 Parallel Sieving ........................................................................................35. Chapter 4 Experimental Results.........................................................................37 4.1 Environment.............................................................................................37 4.2 Results......................................................................................................37 Chapter 5 Conclusion ..........................................................................................40 References ...................................................................................................................41. iv.

(7) List of Tables Table 1..........................................................................................................38 Table 2..........................................................................................................39. v.

(8) Chapter 1 Introduction The RSA Cryptosystem [ 1 ] is one of the most important public-key cryptosystems, and the security of it rests on the fact that it is computationally infeasible to factor a large integer into its component primes. If an efficient algorithm is found that can factor any large integer in polynomial time, the RSA Cryptosystem would be insecure. In this chapter, we will describe some important number-theoretic results, the RSA Cryptosystem, the details of setting up it, etc.. 1.1 Elementary Number Theory In the beginning of this section, we first introduce some basic definitions from elementary group theory.. Definition 1: [ 1 ] For a finite multiplicative group G, define the order of an element g ∈ G to be the smallest positive integer m such that gm = 1. If there are n elements of G, then we say that G is a multiplicative group of order n.. We then proceeds to mention a very important theorem, called the Lagrange’s theorem [ 1 ].. Theorem 1 Suppose G is a multiplicative group of order n, and g ∈ G. Then the order of g divides n.. 1.

(9) From Theorem 1, it is clear that gn = (gm)n/m = 1 for any element g ∈ G. For any positive integer n, let Z n* denote the set of residues modulo n that are relatively prime to n. It can be easily verified that Z n* is a (finite) multiplicative group. The Euler phi-function φ ( n) [ 1 ] is defined to be the number of positive integers not exceeding n and relatively prime to n. That is, Z n* = φ ( n). Given the prime-power factorization of n, a well-known theorem provides a formula to evaluate the value of φ ( n) [ 2 ]:. Theorem 2 Let n = p1a1 p2 a2 ... pk ak be the prime-power factorization of the positive integer n. Then. ⎛. φ ( n) = n ⎜ 1 − ⎝. ⎛ 1 ⎞⎛ 1 ⎞ 1 ⎞ ⎟⎜1 − ⎟ ... ⎜1 − ⎟ . p1 ⎠⎝ p2 ⎠ pk ⎠ ⎝. (1). By using the results above, it is easy to see that g φ ( n ) ≡ 1 (mod n). (2). for any element g ∈ Z n* . This fact is fairly important and essentially relevant to the RSA Cryptosystem.. 1.2 The RSA Cryptosystem The RSA Cryptosystem is one of the most important public-key cryptosystems, which is invented by Ronald Rivest, Adi Shamir, and Leonard Adleman in 1977. In this section, we will describe how it works. Let n = p × q, where p and q are two large 2.

(10) primes. By Theorem 2, it is clear that φ ( n) = (p – 1)(q – 1). An integer d is chosen such that gcd(d, φ ( n) ) = 1. We next compute e = d–1 mod φ ( n).. (3). (Since gcd(d, φ ( n) ) = 1, the inverse of d modulo φ ( n) must exist.) Then, the private key is pair (d, n), and the public key is pair (e, n). To encrypt a message M (where M is a nonnegative integer less than n), the cipher C is computed as C = Me mod n.. (4). To decrypt the cipher C, we compute M' = Cd mod n.. (5). We now verify that M' = M. Since e = d–1 mod φ ( n), we have that ed ≡ 1 (mod φ ( n) ). ⇒ ed = k φ ( n) + 1, for some k ∈ N.. (6). We first consider the case that M ∈ Z n* . Using the result from Section 1.1, it follows that M' = Cd mod n = (Me)d mod n = Med mod n = M kφ ( n ) +1 mod n = ( M φ ( n ) ) M mod n k. = (1)k M mod n = M.. (7). For M ∉ Z n* , if M = 0, it is clear that M' = M. If M ≠ 0, without loss of generality, suppose that M = kp for some k ∈ N. Since M < n, it must be the case that gcd(k, q) = 1, namely gcd(M, q) = 1. Then it follows from the Fermat’s Little Theorem [ 2 ] that 3.

(11) M(q – 1) ≡ 1 (mod q). ⇒ M(q – 1) = k'q + 1, for some k' ∈ N.. (8). Thus we have M' = Cd mod n = M kφ ( n ) +1 mod n = M kφ ( n ) M mod n = (M(q – 1))k(p – 1) M mod n = (k'q + 1)k(p – 1) M mod n k ( p −1). =M. ∑ i =0. Cik ( p −1) ( k ' q ) mod n i. = M + (kp ) q. k ( p −1). ∑ i =1. = M + kn. k ( p −1). ∑ i =1. Cik ( p −1) ( k ') q i −1 mod n i. Cik ( p −1) ( k ') q i −1 mod n i. = M,. (9). as desired.. 1.3 RSA and Factoring Integers The security of the RSA Cryptosystem rests on the fact that it is computationally infeasible to factor a large integer into its component primes. Obviously, if n = p × q can be factored, it is easy to compute φ ( n) = (p – 1)(q – 1) and then compute d = e–1 mod φ ( n) exactly. Therefore, to ensure the security of the RSA Cryptosystem, it is necessary to set n large enough. Nowadays, it is believed that there is no efficient algorithm that can factor any large integer in polynomial time. If a polynomial-time algorithm is found, the RSA Cryptosystem would be insecure.. 4.

(12) Chapter 2 Factoring Algorithms Throughout this chapter, we suppose that n = p × q is the composite integer that we want to factor, where p, q are two large primes, and p and q are roughly the same size. To attempt to factor n, the straightforward method is trial division, which divides n by each prime less than or equal to. n until p or q is found. This method is guaranteed to. find p, q. However, it is computationally infeasible to factor large enough n by using this method. For very large n, we need to use more effective algorithms. Mathematicians have been attempting to find more efficient factoring algorithms for a long time, and a lot of powerful algorithms have been proposed, such as the well-known Pollard’s rho-algorithm and p – 1 algorithm, the continued fraction algorithm, the elliptic curve factoring algorithm, the quadratic sieve factoring algorithm (abbreviated to QS) [ 3 ] and the number field sieve (abbreviated to NFS) [ 4 ]. Because of the restriction of time and computer hardware, we will focus on the quadratic sieve algorithm. The rest of this chapter is organized as follows. Section 2.1 introduces the Dixon’s random squares algorithm, which consists of several essential concepts still used in the QS and NFS (specifically, the concepts of a factor base, being smooth over a factor base, and finding dependencies among vectors over Z 2 ). In Section 2.2, we will give a brief overview of the QS. Finally, Section 2.3 presents the multiple polynomial quadratic sieve (abbreviated to MPQS) [ 3 ], one of the most useful variants of the QS, which is widely employed in practice.. 2.1 The Dixon’s Random Squares Algorithm The basic idea many factoring algorithms use is pretty simple and is described as 5.

(13) follows. Suppose we can find two integers x and y such that x ! ≡ ± y (mod n) and x 2 ≡ y 2 (mod n).. ( 10 ). (x + y) (x – y) = x2 – y2 ≡ 0 (mod n),. ( 11 ). Then. but neither (x + y) nor (x – y) is divisible by n. Therefore gcd(x + y, n) and gcd(x – y, n) must be non-trivial factors of n. This means that n is successfully factored. If integers x and y satisfying ( 10 ) are produced randomly, then there is no guarantee that x ! ≡ ± y (mod n), and the factorization of n may not be yielded. However, what is the probability that x ≡ ± y (mod n)? It can be proved that x ! ≡ ± y (mod n) with probability ≤ 1/2 . In other words, there is at least 1/2 chance that gcd(x + y, n) and gcd(x – y, n) will be nontrivial. By producing enough x and y satisfying ( 10 ), the probability of success can be increased above any desired threshold. The Dixon’s random squares algorithm is a method used to find two integers x and y satisfying ( 10 ). It begins by choosing several random integers ri such that ri2 > n, and then proceeds to compute the values f(ri) = ri2 mod n.. ( 12 ). f(ri) ≡ ri2 (mod n),. ( 13 ). It is clear that for all ri,. and f(ri) ≠ ri2. Therefore the right side of the congruence ( 13 ) is already a perfect square for any ri, and of course multiplying arbitrary ones of the ri2’s will yield a perfect square. The idea is to then find a subset S of these ri’s such that. ∏ f (r ) = y , 2. ri ∈S. i. for some y.. ( 14 ). If this can be done, then by letting. x = ∏ ri , ri ∈S. 6. ( 15 ).

(14) a congruence of the desired type follows 2. ⎛ ⎞ x ≡ ⎜⎜ ∏ ri ⎟⎟ (mod n) ⎝ ri ∈S ⎠ 2. ≡ ∏ ri 2 (mod n) ri ∈S. ≡ ∏ f (ri ) (mod n) ri ∈S. ≡ y2. (mod n).. ( 16 ). Notice that the equation ( 14 ) holds if and only if every prime factor of. ∏ f (r ) is used an even number of times. This then gives us an idea to find S: if ri ∈S. i. we have known the complete factorization of each of the f(ri)’s, it is easy to check to see if the product of some specific f(ri)’s is a square. However, it is clearly difficult to factor each of the f(ri)’s. Therefore, instead of factoring each of the f(ri)’s, we just retain those f(ri)’s, which can be “easily” factored, and use them. The details of doing this will be explained below. For simplicity, we first give the definitions of a factor base and being smooth over a factor base as follows:. Definition 2: A factor base β is a nonempty set of prime integers. An integer α is said to be smooth over the factor base β if all the prime factors of α occur in β (in other words, α factors completely over β ).. Here is an example to illustrate.. Example 1: Suppose that β = {2, 3, 7, 13} is the factor base and α = 504 = 23 × 32 × 7. 7.

(15) Then α is smooth over β because all the prime factors of α (namely, 2, 3, 7) occur in β .. The method of Dixon uses a factor base β = {p1, p2, …, pb}, which is a set of the b smallest primes, for an appropriate value b (it is generally recommended that r log r. 2 2 ). For all ri, we then check to see if f(ri) is smooth over β . If it is, b≈ ln 2 r log 2 r. this ri is said to be “useful”, and is reserved; otherwise we throw this ri out, and try the next one. Suppose W = { rα1 , rα 2 , …, rα m } is a set of ri’s with the property that f( rα j ) is smooth over β for 1 ≤ j ≤ m, and b. f (rα j ) = ∏ pk. ek , j. ( 17 ). k =1. with ek , j ≥ 0, 1 ≤ j ≤ m, 1 ≤ k ≤ b. We then attempt to find a set S satisfying ( 14 ) from the subsets of W. Observe that every subset U of W can be mapped to a r vector z = ( z1 , z2 , ..., zm ) ∈ (Z 2 ) m as follows (where (Z 2 ) m denotes the m-dimensional. vector space over the finite field Z 2 of 2 elements):. ⎧⎪1 if rα j ∈ U zj = ⎨ ⎪⎩0 if rα j ∉U. ( 18 ). for 1 ≤ j ≤ m. It is clear that this mapping is one-to-one and onto, and. ∏ ri ∈U. m. (. f (ri ) = ∏ f (rα j ) j =1. ). zj. ⎛ b e ⎞ = ∏ ⎜ ∏ pk k , j ⎟ j =1 ⎝ k =1 ⎠ m. zj. b ⎛ m e z ⎞ = ∏ ⎜ ∏ pk k , j j ⎟ k =1 ⎝ j =1 ⎠. 8.

(16) m. b. ∑ ek , j z j. = ∏ pk j=1. .. ( 19 ). k =1. As described previously,. ∏ f (r ) is a perfect square if and only if ri ∈U. i. m. ∑e j =1. k, j. z j ≡ 0 (mod 2). ( 20 ). for 1 ≤ k ≤ b. This homogeneous linear system can be written in matrix form as. ⎡ e1,1 e1,2 ⎢e ⎢ 2,1 e2,2 ⎢ M M ⎢ ⎣ eb ,1 eb ,2. L e1,m ⎤ ⎡ z1 ⎤ ⎡0 ⎤ L e2,m ⎥⎥ ⎢ z2 ⎥ ⎢0 ⎥ ⎢ ⎥ ≡ ⎢ ⎥ (mod 2). O M ⎥ ⎢ M ⎥ ⎢M ⎥ ⎥⎢ ⎥ ⎢ ⎥ L eb ,m ⎦ ⎣ zm ⎦ ⎣0 ⎦. ( 21 ). r The question then becomes one of solving the equation ( 21 ). If a solution s of r the equation ( 21 ) is found, the set S can then be constructed according to s . It is a standard result from linear algebra [ 4 ] that if m > b then the equation ( 21 ) has at least Z 2. m −b. > 2 solutions. This means that there must be at least one non-trivial. solution of the equation ( 21 ),which can be used to construct a nonempty set S satisfying ( 14 ). Since the equation ( 21 ) is solved only modulo 2, it can be simplified by replacing the ek , j with (ek , j mod 2) for 1 ≤ j ≤ m, 1 ≤ k ≤ b. There are many efficient algorithms for solving a homogeneous linear system over a finite field, such as Gauss-Jordan elimination [ 5 ], block Lanczos algorithm [ 6 ], and Wiedemann algorithm [ 7 ]. In fact, it spends most of time determining whether f(ri) is smooth over β for all ri, instead of solving the linear system. Therefore, the real question is how to find enough ri with f(ri) smooth over β in an efficient way.. 2.2 The Quadratic Sieve Factoring Algorithm The quadratic sieve factoring algorithm is a well-known algorithm invented by Carl 9.

(17) Pomerance in 1981. It was the fastest known general-purpose factoring algorithm until the number field sieve was proposed, and has been widely used in practice for a long time. Generally speaking, the QS is faster than the number field sieve for numbers having less than about 110 digits. Up to now, the QS is still the algorithm of choice for factoring large integers between 50 and 110 digits. In reality the QS extends the ideas of the Dixon’s random squares algorithm. At its kernel, the QS is essentially the same as the Dixon’s method. There are two major differences between them. The first one is that instead of using the function f(ri) = ri2 mod n, the function f(ri) = ri2 – n. ( 22 ). is used. It is easy to see that for all ri the congruence f(ri) ≡ ri2 (mod n). ( 13 ). still holds even though the function f(ri) has been replaced. Hence the new f(ri) can play the same role the old f(ri) plays. The second difference is in how to obtain integers ri. In the Dixon’s method, we simply choose the ri’s at random. In contrast, the QS uses successive integers as ri’s, such as ri = ⎢⎣ n ⎥⎦ + i, i = 1, 2, …. It looks like that the QS is not much different from the Dixon’s method. But through these slight modifications, some special tricks can be used and the running time becomes dramatically faster. In this section, we describe the details of doing this.. 2.2.1 Setting Up the Factor Base As with the Dixon’s method, the QS also begins by fixing a factor base β = {p1, p2, …, pb}. Then we search for integers ri with f(ri) is smooth over β . However, notice that not any prime can be put into β . For any pk ∈ β , it must be satisfied that there exists at least one ri such that f(ri) is divisible by pk; otherwise there is no f(ri) 10.

(18) divisible by this pk, and putting it into β doesn’t make sense at all. Therefore, for any pk ∈ β , pk | f(ri), for some ri. ⇔ pk | (ri2 – n), for some ri. ⇔ ri2 ≡ n (mod pk), for some ri. ⇔ n is a quadratic residue modulo pk.. ⎛ n ⎞ ⇔ ⎜ ⎟ = 1. ⎝ pk ⎠. ( 23 ). ⎛ n ⎞ Where ⎜ ⎟ denotes the Legendre symbol, which can be evaluated by using the ⎝ pk ⎠ following theorem [ 8 ].. Theorem 3. Suppose po is an odd prime. Then. ⎛ n ⎞ ( po −1) / 2 mod po ⎜ ⎟=n ⎝ po ⎠. ( 24 ). The modular exponentiation of ( 24 ) can be computed efficiently by using the well-known Square-and-Multiply algorithm [ 8 ]. Thus we can decide which odd prime po should be put into β by easily determining whether ( n ( po −1) / 2 mod po) = 1. On the other hand, we should choose the primes of β as small as possible, because the f(ri)’s are intuitively thought more likely smooth over β when the primes of β are smaller. At this point, we can set up our factor base as follows. First, we set β to be an empty set. Then we should put the prime 2 into β since f(ri) = ri2 – n is even as ri is 11.

(19) odd. We then proceed to start at po = 3 and check to see if ( n ( po −1) / 2 mod po) = 1. If it does, then po is added to β , otherwise it is discarded. In either case, the next prime is assigned to po, and the process continues until | β |= b, for an appropriate value b.. 2.2.2 The Sieving Procedure Once β has been set up completely, we begin to determine whether f(ri) is smooth over β for all ri. As described previously, this procedure is the most time-consuming part of this kind of algorithms. Let’s consider how to determine which f(ri) is smooth over β . Obviously the straightforward method is trial division, which divides f(ri) by every prime of β . However, this method is incredibly inefficient. In general, a specific f(ri) is not divisible by most primes of β . Therefore, a lot of time is wasted attempting to divide a specific f(ri) by those primes which don’t actually divide it. In Dixon’s method, it seems that we have no alternative but to do trial division. In fact, the key breakthroughs occur when we change the viewpoint of the operations. Instead of focusing on one fixed f(ri) at a time and trying to divide it by all the primes of β , we fix a prime of β and determine which f(ri) are divisible by it. It is easy to see which f(ri) = ri2 – n is divisible by 2 by determining if ri is odd (because ri2 – n is divisible by 2 if and only if ri is odd). On the other hand, for a fixed odd prime po ∈ β , we need to find all the ri’s with po | (ri2 – n). ⇔ ri2 ≡ n (mod po). ⇔ ri is a solution to the congruence r2 ≡ n (mod po).. ( 25 ). We already know that n is a quadratic residue modulo po and po is an odd prime, so the congruence r2 ≡ n (mod po) has exactly two solutions in Z po , say so,1 and so,2. (Moreover, these two solutions are negatives of each other modulo po, namely so,2 = 12.

(20) po – so,1.) Let so ∈ {so,1, so,2}. Then it is clear that ri is a solution to the congruence r2 ≡ n (mod po). ⇔ ri = so + tpo, t ∈ Z.. ( 26 ). Hence it remains to consider how to compute so,1 and so,2 in a reasonable manner. Fortunately, there is an efficient method called the Shanks-Tonelli algorithm [ 1 ], which can be used to compute these modular square roots efficiently. Since so,1 and so,2 only depend on n and po, when we set up β , we also compute (and store) them for each po in β . Although all the ri’s satisfying ( 26 ) can be found, it is obviously impossible to use all of them. In practice, we pick an interval and just consider the ri’s in this interval. Such an interval is called the sieving interval. To simplify matters, suppose the sieving interval is ⎡ ⎢⎣ n ⎥⎦ + 1, ⎢⎣ n ⎥⎦ + δ ⎤ , and ri = ⎢⎣ n ⎥⎦ + i, for 1 ≤ i ≤ δ . The ⎣ ⎦ bound δ is selected such that it is expected more than b f(ri)’s which correspond to the ri’s within this range will be smooth over β . Then an array of computer memory is allocated, and for i = 1, 2, …, δ , f(ri) = ri2 – n is calculated and stored in the array. Since the ri’s are successive instead of being random, every ri can be mapped to the index of the array element which saves the corresponding f(ri). Suppose the array elements are M[1], M[2], …, M [δ ]. We can store the f(ri)’s in such a way: for each ri, M[i] is assigned to f(ri), namely M[ri – ⎢⎣ n ⎥⎦ ] = f(ri). Therefore, given an ri, we can easily determine which M[l] = f(ri), 1 ≤ l ≤ δ . In the next step of the algorithm, the congruence r2 ≡ n (mod po) is solved for each odd prime po ∈ β . All the ri’s satisfying ⎢ n ⎥ + 1 ≤ ri = so + tpo ≤ ⎢ n ⎥ + δ , t ∈ Z ⎣ ⎦ ⎣ ⎦ are then picked out, and the corresponding M[ri – ⎢⎣ n ⎥⎦ ]’s are divided by po 13. ( 27 ).

(21) repeatedly until their quotients are not divisible by po any more. This procedure is performed for every odd prime po ∈ β . Similarly, for every odd ri, f(ri) is divided by 2 repeatedly until it is not divisible by 2 any more. (Even we can easily divide f(ri) by 2c by doing bitwise right shifts if f(ri) is divisible by 2c.) In the end all the M[l]’s are scanned for which M[l] = 1, 1 ≤ l ≤ δ . M[l] = 1 if and only if f( ⎢⎣ n ⎥⎦ + l) is smooth over β . Consequently, we can find out all the ri’s within the sieving interval with f(ri)’s smooth over β . By using this technique, every division executed is “meaningful”. That is to say, f(ri) is divided by pk if and only if f(ri) is divisible by pk for every prime pk ∈ β . Any blind division trying to divide an f(ri) by the pk which doesn’t evenly divide it. Moreover, the divisions that divide an integer by its prime factor are much faster than the other divisions. Therefore, through omitting the useless divisions, the running time is dramatically speeded up. The approach described in this subsection is called the sieving procedure, which yields the so-called quadratic sieve algorithm.. 2.2.3 Improvements on the QS Although the algorithm has been dramatically improved, the sieving procedure is still the most time-consuming part of the algorithm. There are several methods of accelerating the speed of sieving. One way is simply to set the size of each f(ri) as small as possible. In order to do this, observe that replacing the sieving. δ δ⎤ ⎡ interval ⎡ ⎢⎣ n ⎥⎦ + 1, ⎢⎣ n ⎥⎦ + δ ⎤ by ⎢ ⎢⎣ n ⎥⎦ − , ⎢⎣ n ⎥⎦ + ⎥ can effectively decrease ⎣ ⎦ 2 2⎦ ⎣ the sizes of half the f(ri)’s. Although the f(ri)’s corresponding to the ri’s. δ ⎡ ⎤ within ⎢ ⎣⎢ n ⎦⎥ − , ⎣⎢ n ⎦⎥ ⎥ are negative, we can still factor them (by especially 2 ⎣ ⎦ regarding (–1) as a factor). However, condition ( 14 ) must be still satisfied for some S. 14.

(22) In other words, except that every prime factor of. ∏ f (r ) is used an even number of ri ∈S. times,. i. ∏ f (r ) is necessarily positive, i.e., (–1) of ∏ f (r ) is also used an even ri ∈S. i. ri ∈S. i. number of times. Therefore, the question can be easily solved by adding (–1) to our factor base, and the approach of finding S just works like the Dixon’s method. Besides the method described above, another technique usually used is to predict which f(ri) is smooth over β by using logarithmic operations. Observe that b. f (ri ) = ∏ pk. ek ,i. k =1. b. ⇒ log( f (ri )) = ∑ ek ,i log( pk ) k =1 b. ⇒ log( f ( ri )) − ∑ ek ,i log( pk ) = 0. ( 28 ). k =1. with ek ,i ≥ 0. Thus we can probably predict whether f(ri) is smooth over β as follows. First, we compute log(f(ri)) for each ri in the sieving interval. For every pk ∈ β , we then proceed to subtract log(pk) from log(f(ri)) for those f(ri)’s are divisible by pk. This can be done efficiently because all the ri’s satisfying ( 25 ) can be easily found. If the log(f(ri)) is reduced to 0 by this procedure, the corresponding f(ri) is necessarily smooth over β . However, this event only happens when ek,i = 0, 1 for 1 ≤ k ≤ b. If ek,i > 1, this procedure can not yield the accurate predictions. But if we specify a reasonable threshold and only preserve the f(ri)’s whose log(f(ri))’s are reduced below this threshold, we can eliminate a lot of f(ri)’s which are not smooth over β . We only try to factor the remained f(ri)’s. On the other hand, some f(ri)’s smooth over β may also be eliminated. Therefore, the size of the threshold is a trade-off between eliminating too many “useful” f(ri)’s and reserving too many “useless” f(ri)’s.. 2.3 The Multiple Polynomial Quadratic Sieve 15.

(23) The multiple polynomial quadratic sieve was suggested by Peter Montgomery and is one of the variants of the QS. As the name implies, it uses several polynomial functions instead of just one f(ri) = ri2 – n in the QS. A big problem in the QS is that as ri gets large, f(ri) = ri2 – n also becomes large. Of course, the larger f(ri) is, the less likely it is that f(ri) is smooth over β . For fighting the drift to infinity of f(ri), the MPQS uses several polynomial functions g1(ri), g2(ri), …. Once the values of one polynomial get “too” large, we discard it and use a new one. This procedure not only makes the values of gh(ri) smaller, but also makes the sieving interval and the factor base much smaller. Of course, all this is done to increase the speed of finding the gh(ri)’s smooth over β . In the MPQS, the polynomials must be chosen according to certain conditions. In the subsection below, we then proceed to describe the details of doing this.. 2.3.1 Polynomials Selection Observe that if we use polynomial functions of the form gh(ri) = (ri + bh)2 – n, the values of different gh(ri)’s actually overlap. Hence selecting polynomials in such way doesn’t make sense. The MPQS uses the polynomial functions of the form gh(ri) = ah ri2 + 2bh ri + ch,. ( 29 ). where the coefficients ah, bh, ch are chosen according to the guidelines below. 1. ah is a perfect square, say ah = dh2. 2. Choose 0 ≤ bh < ah such that bh2 ≡ n (mod ah). 3. Choose ch such that bh2 – ah ch = n. (Such a ch must exist because of our choice of bh.) If these can be done, then ah × gh(ri) 16.

(24) = (ah ri)2 + 2(ah ri) bh + ah ch = (ah ri)2 + 2(ah ri) bh + (bh2 – n) = (ah ri + bh)2 – n.. ( 30 ). ah × gh(ri) ≡ (ah ri + bh)2 (mod n).. ( 31 ). gh(ri) ≡ [ d%h (ah ri + bh)]2 (mod n),. ( 32 ). Thus. Moreover. where d%h = dh–1 mod n (assume dh and n are relatively prime). As with the QS, gh(ri) is congruent to a perfect square modulo n, and this is what we want. On the other hand, what about the factor base? Suppose the factor base β = {p1, p2, …, pb}. For any pk ∈ β , the condition must be still satisfied that there exists at least one ri such that gh(ri) is divisible by pk. That is, for any prime pk ∈ β , pk | gh(ri), for some gh and ri.. ( 33 ). For pk = 2, the condition ( 33 ) can always hold by restricting the values of ah and ch. Consider that for any odd prime po ∈ β , if gcd(ah, po) = 1, then po | gh(ri), for some ri. ⇔ po | ah × gh(ri), for some ri. ⇔ po | [(ah ri + bh)2 – n], for some ri. ⇔ (ah ri + bh)2 ≡ n (mod po), for some ri. ⇔ n is a quadratic residue modulo po.. ⎛ n ⎞ ⇔ ⎜ ⎟ = 1. ⎝ po ⎠ ⇔ n ( po −1) / 2 mod po = 1.. ( 34 ). If gcd(ah, po) ≠ 1 (namely gcd(ah, po) = po), there may not exist gh and ri such that gh(ri) is divisible by po. However, this can be avoided by choosing ah such that for every 17.

(25) odd prime po ∈ β , ah is not divisible by po. Besides this method, if we choose ah to be a power of a prime, there is at most one odd prime in β such that gh(ri) is never divisible by it. Therefore, the procedure used to set up the factor base in the QS can be also used in the MPQS.. 2.3.2 The Details of Choosing the Coefficients The MPQS chooses gh(ri)’s to custom fit not only the number n, but also the length of the sieving interval. Suppose we use the sieving interval [ −δ , δ ] of length 2δ before we change gh(ri). Consider gh(ri) = ah ri2 + 2bh ri + ch 2 ⎛ ⎛ bh ⎞ ⎛ bh ⎞ ⎞ bh 2 − ah ch 2 = ah ⎜ ri + 2ri ⎜ ⎟ + ⎜ ⎟ ⎟ − ⎜ ah ⎝ ah ⎠ ⎝ ah ⎠ ⎟⎠ ⎝ 2. ⎛ b ⎞ n = ah ⎜ ri + h ⎟ − . ah ⎠ ah ⎝. ( 35 ). We would like to make the values of g h (ri ) to be as small as possible on the sieving interval. One way to do this is to have the minimum and maximum values of gh(ri) over [ −δ , δ ] be roughly the same in absolute values, but be opposite in sign. It is clear that the minimum value of gh(ri) is g h (− ah, i.e., –1 < −. bh n ) = − . Since we choose 0 ≤ bh < ah ah. bh ≤ 0, the minimum value of gh(ri) ah. over [ −δ , δ ] is g h (−. bh n ) = − . Moreover, the maximum value of gh(ri) ah ah. over [ −δ , δ ] appears at ri = δ , and it is 2. ⎛ b ⎞ n g h (δ ) = ah ⎜ δ + h ⎟ − ah ⎠ ah ⎝ 18.

(26) ≈ ah δ 2 −. n ah. (ahδ ) 2 − n = . ah. ( 36 ). As described above, we expect. bh (ahδ ) 2 − n n = g h (− ) ≈ g h (δ ) = . ah ah ah ⇒ n ≈ ( ah δ ) 2 − n .. ⇒ ahδ ≈ 2n . ⇒ ah ≈ ⇒ dh ≈. 2n. δ. .. 2n. δ. .. ( 37 ). This then helps us to select a suitable dh. Recall that in subsection 2.3.1, the coefficients ah, bh, ch must be chosen according to three guidelines. The condition 1 can be easily satisfied. If the condition 2 has been satisfied, the condition 3 can be also satisfied by choosing ch bh 2 − n . Therefore, the real question is how to choose bh according to the condition = ah 2. To do this, n must be a quadratic residue modulo ah. This is true if and only if n is a quadratic residue modulo d for every prime factor d of ah [ 8 ], i.e., for every prime d with d | ah, ⎛n⎞ ⎜ ⎟ = 1. ⎝d ⎠. ( 38 ). Hence, we would like to choose ah with its factorization known (namely, choose dh with its factorization known, because ah = dh2). For convenience, we choose dh as a prime close to. 2n. δ. ⎛ n ⎞ such that ⎜ ⎟ = 1. ⎝ dh ⎠ 19.

(27) Once dh has been chosen, we then proceed to solve the congruence r2 ≡ n (mod dh2),. ( 39 ). and set bh to be one of the modular square roots. If the congruence r2 ≡ n (mod dh). ( 40 ). can be solved, we can also compute the solutions of the congruence ( 39 ) by the following theorem [ 9 ].. Theorem 4. (Hensel’s Lemma). Suppose that f(x) is a polynomial with integer coefficients and that k is an integer with. k ≥ 2. Suppose further that r is a solution of the congruence f(x) ≡ 0 (mod pk – 1). Then, (i). if f ′( r )! ≡ 0 (mod p), then there is a unique integer t, 0 ≤ t < p, such that f(r + tpk – 1) ≡ 0 (mod pk), given by ⎛ f (r ) ⎞ t ≡ – f% ′(r ) ⎜ k −1 ⎟ (mod p), ⎝ p ⎠ where f% ′(r ) is an inverse of f ′(r ) modulo p;. (ii) if f ′( r ) ≡ 0 (mod p) and f(r) ≡ 0 (mod pk), then f(r + tpk – 1) ≡ 0 (mod pk) for all integers t; (iii) if f ′( r ) ≡ 0 (mod p) and f(r)! ≡ 0 (mod pk), then f(x) ≡ 0 (mod pk) has no solutions with x ≡ r (mod pk – 1).. Suppose f(r) = r2 – n, sh is a solution of the congruence ( 40 ) (namely, the congruence f(x) ≡ 0 (mod dh)). By the Theorem 4, we can easily calculate one solution of congruence ( 39 ) as follows. First, compute. ⎛ f ( sh ) ⎞ th = – f% ′( sh ) ⎜ ⎟ mod dh, ⎝ dh ⎠ 20. ( 41 ).

(28) where f% ′( sh ) is an inverse of f ′( sh ) modulo dh, i.e.,. f% ′( sh ) = (2sh)– 1 mod dh.. ( 42 ). sh' = sh + th dh mod dh2. ( 43 ). Then. is a solution of the congruence ( 39 ) (namely, the congruence f(x) ≡ 0 (mod dh2)). As described previously, the Shanks-Tonelli algorithm can be used to compute the modular square roots of the congruence ( 40 ). However, if we choose dh by using the tricks below, this work can be done more efficiently. Suppose we choose dh as a. ⎛ n ⎞ prime with ⎜ ⎟ = 1 and ⎝ dh ⎠ dh ≡ 3 (mod 4). If this can be done, then n ( d h −1) / 2 ≡ 1 (mod dh) and. (n. ). ( d h +1) / 4 2. ( 44 ). dh + 1 is an integer. Thus, 4. ≡ n ( dh +1) / 2 (mod dh) ≡ n n ( d h −1) / 2 (mod dh) ≡ n (mod dh).. ( 45 ). That is, ( n ( d h +1) / 4 mod dh) is a modular square root of the congruence ( 40 ). Therefore, we can set sh to be ( n ( d h +1) / 4 mod dh) and use it to compute sh'.. 2.3.3 Sieving Just as the QS, we need to solve the congruence gh(r) ≡ 0 (mod po) for each odd prime po in the factor base β . Nothing but whenever we use a new polynomial as gh(ri), we need to do this work again for the new polynomial. Fortunately, the congruence ah r2 + 2bh r + ch ≡ 0 (mod po) 21. ( 46 ).

(29) can be easily solved by using the standard formula for solving a quadratic polynomial. (Recall that there is at most one po with gcd(ah, po) ≠ 1, and we would not solve gh(r) ≡ 0 (mod po) for this po.) r = (2ah)–1[–2bh ± ((2bh)2 – 4ah ch)1/2] mod po = 2–1ah–1[–2bh ± 2(bh2 – ah ch)1/2] mod po = ah–1[–bh ± n1/2] mod po.. ( 47 ). Since gcd(ah, po) = 1, (ah–1 mod po) exists and we can always find it. Moreover, the square roots of n modulo po (n1/2 mod po) can be computed by using the Shanks-Tonelli algorithm. Therefore, the sieving procedure of the MPQS just works the same way as the QS, besides using multiple polynomials instead of a single one.. 22.

(30) Chapter 3 The Modified Multiple Polynomial Quadratic Sieve As described above, sieving procedure is the most time-consuming part of the MPQS. Specifically, it spends most of time doing trial division in order to determine which g(ri) is smooth over β . Trial division must be applied because we don’t know how many times pj divides a given g(ri) (if g(ri) is divisible by pj), for each pj ∈ β . However, if we can explicitly compute the number of times (pj divides a given g(ri)) without doing any trial division, is it possible to improve the MPQS? Notice that if this can be done, we can determine whether a given g(ri) is smooth over β by doing logarithmic operations. We illustrate this technique with a small example.. Example 2:. Suppose that g(ri) = 504 = 23 × 32 × 7 and β = {2, 3, 7, 13}. Then g(ri) is smooth over β because. 504 = 1. On the other hand, we can conclude the same 2 × 32 × 7 3. result according to the reason that 504 is divisible by 23, 32, 7 and log(504) – [3 × log(2) + 2 × log(3) + log(7)] = 0.. ( 48 ). This idea is fairly simple, and we will particularly mention it latter in this chapter. Generally speaking, trial division (of large numbers) spends more time than logarithmic operations. Therefore, it remains to consider how to compute the number of times pj divides a given g(ri) without doing any trial division. In this Chapter, we 23.

(31) will describe our methods of doing this. The remaining sections of this chapter are organized as follows. Section 3.1 introduces the basic ideas of our methods. In Section 3.2, we discuss how to compute the square roots of n modulo pok for each odd prime po ∈ β . The results of doing this are very important and will be used in the following steps. In Section 3.3, we describe how to solve the congruence g(r) ≡ 0 (mod pok) for the g(r) in the MPQS; and how to apply these results to the sieving procedure. Finally, in order to make the MPQS more practical, Section 3.4 provides a scheme to parallelize the sieving procedure.. 3.1 Motivation for the Modified Multiple Polynomial Quadratic Sieve Recall from Subsection 2.3.3 that in the sieving procedure we first solve the congruence g(r) = a r2 + 2b r + c ≡ 0 (mod po). ( 49 ). for each odd prime po in the factor base β . By doing this, we can find all the ri with g(ri) divisible by po. At this point, we already know which g(ri) is divisible by po, but how do we know the exponent of po in the prime power factorization of g(ri)? It might appear to be necessary to divide g(ri) by po repeatedly until its quotient is not divisible by po (i.e. do trial division). Recall that the maximum values of g (ri ) are about δ as. n . Thus, it is intuitively reasonable that many of g (ri ) are almost as large 2. n , and it would spend a lot of time to divide each g(ri) by its prime factors.. However, if the exponent of po (in the factorization of g(ri)) can be derived without doing any trial division, this shift may lead to a speed-up. 24.

(32) Now suppose that we can find the solutions of the congruence g(r) = a r2 + 2b r + c ≡ 0 (mod pok). ( 50 ). for any positive integer k, and So,k = {ri | g(ri) ≡ 0 (mod pok)} = {ri | g(ri) is divisible by pok}.. ( 51 ). Notice that if g(ri) is divisible by pok + 1, it is also divisible by pok. Thus it is clear that So,k + 1 ⊆ So,k, k = 1, 2, …. Consider Do,k = So,k – So,k + 1 = {ri | g(ri) = t pok, gcd(t, po) = 1}.. ( 52 ). For a particular k, if Do,k can be found (i.e. So,k, So,k + 1 can be found), we can find all the g(ri) divisible exactly by pok but not divisible by pok + 1. In other words, we can find all the g(ri) in whose prime power factorization pok appears. Of course, the prerequisite is that the congruence ( 50 ) can be solved for any positive integer k. In the next section, we briefly discuss how to solve the simplest case of the congruence ( 50 ).. 3.2 Square Roots of n Modulo pok Suppose g(r) = r2 – n,. ( 53 ). a r2 + 2b r + c.. ( 54 ). the simplest polynomial of the form. (Of course, a, b, c must be chosen according to the guidelines in the MPQS.) In fact g(r) is a special form of these polynomials, and it plays an important role in our method. We now consider how to solve the congruence g(r) = r2 – n ≡ 0 (mod pok) 25. ( 55 ).

(33) for any positive integer k.. ⎛ n ⎞ Recall that for every odd prime po ∈ β , ⎜ ⎟ = 1. Therefore, for any positive ⎝ po ⎠ integer k there are two square roots of n modulo pok according to the theorem below [ 8 ].. Theorem 5. Suppose that p is an odd prime, e is a positive integer, and gcd(a, p) = 1. Then the ⎛a⎞ congruence y2 ≡ a (mod pe) has no solutions if ⎜ ⎟ = −1, and two solutions (modulo ⎝ p⎠ ⎛a⎞ pe) if ⎜ ⎟ = 1. ⎝ p⎠. Since there are exactly two modular square roots, it is clear that they are negatives of each other modulo pok, and we need to compute just one of them. When k = 1, as described previously the square roots of n modulo po can be computed efficiently by using the Shanks-Tonelli algorithm. When k ≥ 2, the Hensel’s Lemma is applied. Suppose uk – 1 is a solution of the congruence g(r) ≡ 0 (mod pok – 1). Then uk – 1! ≡ 0 (mod po).. ( 56 ). To see this, consider that (uk – 1)2 – n ≡ 0 (mod pok – 1) ⇒. (uk – 1)2 – n ≡ 0 (mod po).. ( 57 ). If uk – 1 ≡ 0 (mod po), it implies that n ≡ 0 (mod po), which is a contradiction since gcd(n, po) = 1. Hence uk – 1! ≡ 0 (mod po), and g ′(uk −1 ) = 2 uk −1 ! ≡ 0 (mod po). (Notice that po is an odd prime.) Therefore, case (i) of Hensel’s Lemma always 26. ( 58 ).

(34) applies. That is, uk = (uk – 1 + tk – 1 pok – 1) is a solution of the congruence g(r) ≡ 0 (mod pok), given by. ⎛ g (uk −1 ) ⎞ tk – 1 ≡ – g% ′(uk −1 ) ⎜ (mod po), k −1 ⎟ ⎝ po ⎠. ( 59 ). where g% ′(uk −1 ) is an inverse of g ′(uk −1 ) modulo po. We now consider the solution uk + 1 = (uk + tk pok) of the congruence g(r) ≡ 0 (mod pok + 1). First, g ′(uk ) = 2 uk = 2 (uk – 1 + tk – 1 pok – 1) ≡ 2 uk – 1 (mod po). ≡ g ′(uk −1 ) (mod po).. ( 60 ). Thus g% ′(uk ) = g% ′(uk −1 ), and we don’t need to compute g% ′(uk ) repeatedly once g% ′(uk −1 ) is computed. By extending this result, it is clear that g% ′(uk ) = g% ′(u1 ) for any k ≥ 1 (where u1 is a solution of the congruence g(r) ≡ 0 (mod po)). Suppose qk −1 =. g (uk −1 ) . po k −1. ( 61 ). Then we can compute qk as follows: qk =. g (u k ) po k. =. uk 2 − n po k. =. (uk −1 + tk −1 po k −1 ) 2 − n po k. =. 2 uk −1 tk −1 po k −1 + (tk −1 po k −1 ) 2 + ((uk −1 ) 2 − n) po k. 2 uk −1 tk −1 po k −1 + g (uk −1 ) po k 2 uk −1 tk −1 + qk −1 = (tk −1 ) 2 po k − 2 + . po. = (tk −1 ) 2 po k − 2 +. When k ≥ 3, 27. ( 62 ).

(35) qk ≡. 2 uk −1 tk −1 + qk −1 (mod po). po. ( 63 ). These results then provide an efficient method to evaluate qk and (qk mod po) through qk – 1 when k ≥ 3. In the rest of this section, we discuss the size of uk we computed. Of course, we wish to make each uk as small as possible, i.e. 1 ≤ uk ≤ po k − 1,. ( 64 ). 1 ≤ uk −1 ≤ po k −1 − 1.. ( 65 ). 0 ≤ tk −1 ≤ po − 1,. ( 66 ). for k ≥ 1. Assume. If we choose tk – 1 with. then 0 ≤ tk −1 po k −1 ≤ ( po − 1) po k −1.. ⇒ 1 + 0 ≤ uk −1 + tk −1 po k −1 ≤ ( po k −1 − 1) + ( po − 1) po k −1. ⇒ 1 ≤ uk ≤ po k − 1.. ( 67 ). Therefore we take u1 ∈ Z po and. where qk =. tk = – g% ′(u1 ) qk mod po,. ( 68 ). 1 ≤ u k ≤ po k − 1. ( 69 ). g (u k ) . Then po k. follows for any k ≥ 1.. 3.3 Modified Sieving Procedure Suppose 28.

(36) g(r) = a r2 + 2b r + c,. ( 70 ). where coefficients a, b, c satisfy the guidelines in the MPQS. In the beginning of this section, we first consider the solutions to the congruence g(r) = a r2 + 2b r + c ≡ 0 (mod pok),. ( 71 ). for any positive integer k. Throughout this section, we will suppose that gcd(a, po) = 1. (Recall that there is at most one po with gcd(a, po) ≠ 1.) Since gcd(a, po) = 1, the congruence ( 71 ) can be solved by using the standard formula for solving a quadratic polynomial. That is, r = (2a)–1[–2b ± ((2b)2 – 4a c)1/2] mod pok = 2–1a–1[–2b ± 2(b2 – a c)1/2] mod pok = a–1[–b ± n1/2] mod pok = a–1[–b ± no( k ) ] mod pok,. ( 72 ). where no( k ) denotes the square root of n modulo pok. Recall that n is a quadratic residue. ⎛ n ⎞ modulo po and the Legendre symbol ⎜ ⎟ = 1. According to Theorem 5, there are ⎝ po ⎠ exactly two square roots of n modulo pok, namely no( k ) and pok – no( k ) . Therefore, the congruence ( 71 ) has exactly two solutions modulo pok, say. so( k,1) = a–1[–b + no( k ) ] mod pok. ( 73 ). so( k,2) = a–1[–b – no( k ) ] mod pok.. ( 74 ). and. We already know from Section 3.2 that no( k ) can be computed efficiently by using the Hensel’s Lemma. We now consider 29.

(37) So,k = {r | g(r) is divisible by pok}.. ( 75 ). It is clear that So,k = {r | g(r) ≡ 0 (mod pok)} = {so( k ) + t pok | so( k ) ∈{so( k,1) , so( k,2) }, t ∈ Z}.. ( 76 ). Once So,k, So,k + 1 are found, Do,k = So,k – So,k + 1. ( 77 ). is found. As described in Section 3.1, we can find all the g(r) divisible exactly by pok but not divisible by pok + 1. If we wish to find Do,1, Do,2, …, Do , ko −1 , we need to find So,1, So,2, …, So ,ko in advance. It might appear to be necessary to first compute so( k,1) , so( k,2) , for k = 1, 2, …, ko. Fortunately, we just need to compute so( k,1o ) , so( k,2o ) . To see this, notice that. so( k,1o ) ≡ so( k,1) (mod pok). ( 78 ). so( k,2o ) ≡ so( k,2) (mod pok),. ( 79 ). and. for k ≤ ko . We will give a brief proof below. For clarity, suppose that ao( k ) is an inverse of a modulo pok. Then for k ≤ ko , a ao( ko ) ≡ 1 (mod po ko ). ⇒ a ao( ko ) ≡ 1 (mod pok). ⇒ ao( ko ) ≡ ao( k ) (mod pok).. Moreover,. (n ). ( ko ) 2 o. ≡ n (mod po ko ).. 30. ( 80 ).

(38) ⇒. (n ). ( ko ) 2 o. ≡ n (mod pok).. ( 81 ). In other words, no( ko ) is a square root of n modulo pok. Without loss of generality, suppose no( ko ) ≡ no( k ) (mod pok).. ( 82 ). From the discussion above, it follows that. so( k,1o ) = ao( ko ) [–b + no( ko ) ] mod po ko ≡ ao( ko ) [–b + no( ko ) ] (mod pok) ≡ ao( k ) [–b + no( k ) ] (mod pok). ≡ so( k,1) (mod pok).. ( 83 ). so( k,2o ) ≡ so( k,2) (mod pok).. ( 84 ). Similarly, it can be proved that. Therefore, we obtain the following result: for k ≤ ko , So,k = {so( k ) + t pok | so( k ) ∈{so( k,1) , so( k,2) }, t ∈ Z} = {so( ko ) + t pok | so( ko ) ∈{so( k,1o ) , so( k,2o ) }, t ∈ Z}.. ( 85 ). That is to say, if so( k,1o ) , so( k,2o ) have been computed, So,1, So,2, …, So ,ko are found (so are Do,1, Do,2, …, Do ,ko −1 ). so( k,1o ) , so( k,2o ) can be computed by using the formula ( 72 ). Since no( ko ) only depends on n, po and ko, when we set up β , we also compute (and store) it for each po in β . Whenever we compute so( k,1o ) , so( k,2o ) for new polynomials, the Hensel’s Lemma would not be used. We now describe the proposed scheme. To simplify matters, suppose that the sieving interval is [1 − δ , δ ] , ri = i – δ , and all the g(ri)’s which correspond to the ri’s 31.

(39) within this range are divisible by po at most ko times. We first evaluate so( k,1o ) , so( k,2o ) for each odd prime po ∈ β . By doing this, So,1, So,2, …, So ,ko can be found. Therefore, we already know which g(ri) is divisible by po exactly k times, for any k ≤ ko . Then we can determine which g(ri) is smooth over β without doing any trial division. To see this, suppose that g(ri) is divisible by pj exactly ei,j times with ei , j ≥ 0, for each pj ∈ β . Then it is clear that g(ri) is smooth over β . b. ⇔ g (ri ) = ∏ p j i , j . e. j =1. b. ⇔ log( g (ri )) = ∑ ei , j log( p j ). j =1 b. ⇔ log( g (ri )) − ∑ ei , j log( p j ) = 0.. ( 86 ). j =1. Now we can apply this result to the sieving procedure as follows. First, an array of size 2δ is allocated, and suppose the array elements are M[1], M[2], …, M [2δ ]. We then evaluate the logarithm of g(ri) for each ri, and assign it to M [i ], namely M [ri + δ ] = log(g(ri)). For every ri ∈ So,k and 1 – δ ≤ ri ≤ δ , log(po) is subtracted from M [ri + δ ], where k ≤ ko . (That is, for every ri ∈ Do,k and 1 – δ ≤ ri ≤ δ , k × log(po) is subtracted from M [ri + δ ], where k < ko. Moreover, since all the g(ri)’s we consider are divisible by po at most ko times, for every ri ∈ S o ,ko and 1 – δ ≤ ri ≤ δ , ko × log(po) is subtracted from M [ri + δ ].) This procedure is performed for every odd prime po ∈ β . Similarly, if g(ri) is divisible by 2 at most k' times, k' × log(2) must be subtracted from M [ri + δ ]. In fact, k' can be easily determined by scanning g(ri) from the least significant bit towards more significant bits, until the first 1 bit is found. Finally, all the M[i]’s are scanned for which M[i] = 0, 1 ≤ i ≤ 2δ . Clearly, M[i] = 0 if and only if g(ri) is 32.

(40) smooth over β . By using this technique, we can determine which g(ri) is smooth over β without doing any trial division (of large numbers). At this point, we know that after sieving M[i] = 0 if and only if g(ri) is smooth over β . However, if g(ri) is not smooth over β , how large M[i] is? The answer of this question is very useful. Generally, a logarithm of a positive integer is not a finite decimal. Therefore, when we do logarithmic operations, inaccuracy may appear. Hence, we need a reasonable bound to estimate which M[i] actually equals 0. The following result is fairly simple, and we will give a brief proof of it.. Theorem 6. Suppose that the factor base β = {p1, p2, …, pb}, pb is the maximum prime in β , pb < d=. a (in general, this condition holds), and g(ri) is divisible by pj at most ei,j times,. where ei , j ≥ 0, j = 1, 2, …, b. Then, if g(ri) is not smooth over β ,. ⎛ b e ⎞ qi = g(ri) / ⎜ ∏ p j i , j ⎟ > pb. ⎝ j =1 ⎠. Proof. ( 87 ). We will prove this by assuming that qi ≤ pb,. ( 88 ). and obtain a contradiction. Since g(ri) is not smooth over β , it must be the case that qi ≠ 1. Suppose q is a prime factor of qi. Since g(ri) is divisible by qi, it is also divisible by q. As described in Section 2.3 (because q ≤ pb < d and d is prime, gcd(a, q) = ⎛n⎞ gcd(d2, q) = 1), n is a quadratic residue modulo q, i.e. ⎜ ⎟ = 1. Recall from Section ⎝q⎠ 2.2 that ⎛ n ⎞ β = {pj | pj is prime, pj ≤ pb and ⎜ ⎟ = 1 }. ⎜p ⎟ ⎝ j⎠ 33. ( 89 ).

(41) Hence, it must be the case that q ∈ β . Without loss of generality, suppose q = p1. Since. ⎛ b e ⎞ g(ri) = qi × ⎜ ∏ p j i , j ⎟ ⎝ j =1 ⎠. ( 90 ) e +1. and q = p1 is a prime factor of qi, it follows that g(ri) is divisible by p1 i ,1 . Thus, we obtain a contradiction, because g(ri) is divisible by p1 at most ei,1 times by our assumption.. From Theorem 6, it is not difficult to see that after sieving M[i] > log(pb) if and only if g(ri) is not smooth over β . In the rest of this section, we briefly discuss how large ko should be for each odd prime po ∈ β . Recall that all the g(ri)’s which correspond to the ri’s within the sieving interval are divisible by po at most ko times, and the maximum values of g (ri ) are about δ. n . Hence, it seems to be reasonable to set ko to be 2. ⎢ ⎛ n ⎞⎥ ⎢1 ⎢log po ⎜⎜ δ ⎟⎟ ⎥ ≅ ⎢⎣log po (δ ) ⎥⎦ + ⎢ log po ( n ) − log po ( 2 ) 2 ⎠ ⎥⎦ ⎣2 ⎢⎣ ⎝. (. )⎥⎥⎦ .. ( 91 ). However, this value of ko is not appropriate. It goes without saying that the smaller ko is, the less time it spends to compute any quantity related to ko. (e.g. po ko and a–1 mod po ko ). But if ko is too small, a lot of g(ri)’s smooth over β would be considered not smooth over β . Thus there is a trade-off. Our idea is to take ko = ⎢⎣ log po ( 2δ ) ⎥⎦. ( 92 ). po ko ≤ 2δ < po ko +1.. ( 93 ). such that. Then there is at most one ri with 1 – δ ≤ ri ≤ δ satisfies 34.

(42) ri ≡ so( k,1o +1) (mod po ko +1 ).. ( 94 ). Similarly, there is at most one ri with 1 – δ ≤ ri ≤ δ satisfies ri ≡ so( k,2o +1) (mod po ko +1 ).. ( 95 ). Therefore, there are at most two ri’s with g(ri)’s divisible by po ko +1. Even though there are two g(ri)’s divisible by po ko +1 , they may not be smooth over β . Hence, most of the g(ri)’s smooth over β would not be eliminated. In reality, ko can be set to be smaller according to the properties of g(ri)’s smooth over β . For example, if the exponents of most g(ri)’s (smooth over β ) are always small, we can set ko to be much smaller.. 3.4 Parallel Sieving In order to make the MPQS more practical, we parallelize the sieving procedure. In general, the way this is done is to partition the sieving interval into several subintervals, and then each processor sieves over a different subinterval. To make the implementation simple, we propose another scheme in this section. Our idea is to make each processor use different quadratic polynomial functions. In fact, this can be easily done if each processor uses different coefficients dh. Here we will use the schemes described in the last of Subsection 2.3.2. Recall that dh ≡ 3 (mod 4). Suppose the MPQS sieves in parallel by using t computers. Then, the jth computer uses the dh satisfying that dh is prime and dh = 4(k t + j) + 3,. ( 96 ). where j = 1, 2, …, t, k ∈ Z. It is clear that dh ≡ 3 (mod 4). Moreover, it is very easy to prove that different computers never use the same dh’s, and this is what we want. The drawback of this method is that t can not be a multiple of 3. Since the tth computer uses the dh satisfying 35.

(43) dh = 4(k t + t) + 3 = 4t(k + 1) + 3,. ( 97 ). if t is a multiple of 3, dh is necessarily a multiple of 3 for any k ∈ Z. However, dh needs to be prime. Therefore, t can not be a multiple of 3 in this method. Fortunately, it is not difficult to prevent t from being a multiple of 3. Hence this method is indeed practical.. 36.

(44) Chapter 4 Experimental Results In our research, we use the program called the GQS [ 10 ], which is developed by Professor D. J. Guan in order to implement the MPQS. The GQS is written with the C language and based on the GMP [ 11 ]. (GMP is a library for arbitrary precision arithmetic, and performs very well on most computers.) We modified the GQS with parallel sieving, and successfully factored a 100-digit number n (where n = p × q, p, q are prime, and p and q are roughly the same size), by distributing the computations to 32 workstations in the department of CS in the NCTU (National Chiao Tung University). However, we do not implement the main idea described in Chapter 3. In the rest of this chapter, we will present the experimental results we obtained.. 4.1 Environment Throughout our experiments, the GQS was performed on the workstations in the department of CS in the NCTU. Each workstation is equipped with AMD Athlon XP 2700+ CPU (running at 2.2 GHz on average), 2 GB memory and 512 MB disk space total. Moreover, they are running operating system RedHat Linux 9.0. The sieving procedure needs to allocate many memories to sieve, and take a lot of disk space to save the sieving results. Therefore, it is important to reserve large enough memory and disk space.. 4.2 Results The asymptotic running time of the MPQS is [ 1 ]:. (. O e(1+ o (1)). ln( n )ln(ln( n )). 37. ).. ( 98 ).

(45) The notation o(1) denotes a function of n that approaches 0 as n → ∞. Formula ( 98 ) can be used to estimate the time required for factoring n using one personal computer. By using one workstation, we have successfully factored the numbers having less than about 50-70 digits. We tabulate the estimated times and the execution times for some values of n in Table 1. Table 1. Estimated running times of. Execution times involving. the MPQS in one PC. one workstation. log10(n) / log2(n). 50 / 166. 0.01 (hours). 0.16 (hours). 60 / 199. 0.16 (hours). 0.11 (hours). 70 / 233. 2.00 (hours). 6.00 (hours). Execution times involving one workstation. The results show that the execution time is not close to the estimation when n has more than about 70 digits. We have also factored larger numbers in the range 80-100 digits. These were done by using parallel sieve on 32 workstations. The execution times are given as follows:. 38.

(46) Table 2. Estimated running times of. Execution times involving. the MPQS in one PC. 32 workstations. log10(n) / log2(n). 80 / 266. 22 (hours). 1.4 (hours). 82 / 272. 36 (hours). 3.4 (hours). 85 / 282. 72 (hours). 6.7 (hours). 90 / 299. 210 (hours). 11.3 (hours). 100 / 332. 80 (days). 6.6 (days). Execution times involving 32 workstations. The results show that sieving by using 32 workstations is not 32 times faster than using one PC. Moreover, the larger n is, the less the speedup is.. 39.

(47) Chapter 5 Conclusion In this paper, we present the methods to enhance the sieving procedure of the MPQS. The advantage of our methods is that it doesn’t need to do a lot of trial division for large numbers. Conversely, our methods need to do a lot of addition and multiplication for smaller numbers. Therefore, this shift may improve the MPQS. For factoring RSA moduli, the NFS is recently the most-used algorithm. A lot of RSA Challenge Numbers were successfully factored by using the NFS. The asymptotic running time of the NFS is [ 1 ]:. (. O e(1.92+ o (1))(ln( n )). 1/ 3. (ln(ln( n )))2 / 3. ),. ( 99 ). which is faster than the MPQS for numbers having more than about 125-130 digits. In fact, the NFS is improved upon by the MPQS, and it still uses the essential concepts of the MPQS. Hence, our ideas can be also applied to the NFS. If the MPQS can be improved by using our methods, the NFS can be also improved. On the other hand, the parameters of our methods (such as the size of the factor base and the length of the sieving interval) can be optimized to reduce the running time. Thus the complexity of the algorithm may be actually smaller.. 40.

(48) References [1] [2] [3] [4]. [5] [6] [7]. [8] [9]. [ 10 ] [ 11 ]. Douglas R. Stinson, Cryptography: Theory and Practice, 2nd Edition, Chapman & Hall/CRC, 2002. Kenneth H. Rosen, Elementary Number Theory and Its Applications, 5th Edition, Pearson Addison Wesley, 2004. Carl Pomerance, “The Quadratic Sieve Factoring Algorithm”, University of Georgia, 1998. A. K. Lenstra, H. W. Lenstra, Jr., et al, “The Development of the Number Field Sieve”, Lecture Notes in Mathematics, vol. 1554, Springer-Verlag, 1993. Stephen H. FriedBerg, et al, Linear Algebra, 4th Edition, Pearson Prentice Hall, 2002. Otto Bretscher, Linear Algebra with Applications, 3rd Edition, Pearson Prentice Hall, United States of America, 2004. P. L. Montgomery, “A block Lanczos Algorithm for finding dependencies over GF(2)”, Advances in Cryptology – EUROCRYPT ’95, vol. 921, pp. 106-120, 1995. D. H. Wiedemann, “Solving Sparse Linear Equations over Finite Fields”, IEEE Trans. Information Theory, IT-32, pp. 54-62, 1986. Ramanujachary Kumanduri, Christina Romero, Number Theory with Computer Applications, 1st Edition, Prentice Hall, Upper Saddle River, New Jersey, 1997. D. J. Guan, “Experience in Factoring Large Integers Using Quadratic Sieve”, 2003. The GNU MP Bignum Library, http://www.swox.com/gmp/.. 41.

(49)