Sieving

Once our implementation of the Gaussian integer method was working, we began searching for (c₁,c₂) pairs of integers whose residues were smooth. All values of c₁ in the interval [1,95000] were considered. We chose to look at only positive values of c₁ for simplicity; negative values of c₁ are acceptable as far as the algorithm is concerned. Our sieve size (the set of considered c₂ values for each value of c₁) was $3 \times 10^5$. Again, only positive c₂ values were considered. These choices were definitely not optimal.

Since -1 is a quadratic nonresidue modulo p, while -2 is a residue, we chose r = -2. (This simplified the algorithm, since $Q(\sqrt{-2})$has class number 1.) By factoring x²+2 modulo p, we obtained a value of S such that $S^2 \equiv -2 \bmod p$. We computed the convergents to the continued fraction of S/p, and selected Vand T as follows:

$$\begin{eqnarray*}V & = & 48454067936694480117959482723, \\ T & = & 22760185083691921160273336139. \end{eqnarray*}$$

These specific values of T and V were chosen from several candidate pairs (T_i,V_i), $T_i,V_i \lesssim \sqrt{p}$ because they satisfy the desirable property T² + 2V² = p.

For each fixed value of c₁, the sieve initialized an integer array C of size $3 \times 10^5$, representing the possible values of c₂. For each prime power q^h, value $d \equiv c_1VT^{-1} \bmod q^h$ was computed. The contents of the array were incremented by $\lfloor 1000 \cdot \ln q \rfloor$ at every location $c_2 \equiv d \bmod q^h$. (The factor of 1000 and the floor function $\lfloor \cdot \rfloor$ are used so that we may maintain three decimal digits of accuracy yet perform integer, instead of floating point, calculations.) When the loop over prime powers q^h had finished, C[c₂] contained the (approximate) real logarithm of the smooth part of the residue c₁V-c₂T.

After all of the small primes had been sieved for a fixed value of c₁, any c₂ for which C[c₂] (i.e., the sum of the logarithms of its small prime factors) was close to the actual (real) logarithm of the residue c₁V - c₂T was considered ``interesting.'' If

$$\left\lfloor 1000 \cdot \ln (c_1V-c_2T) \right\rfloor - C[c_2] < 18420,$$ (10)

then the (c₁,c₂) pair was considered partially-smooth and assumed to have at most one ``moderate-size'' prime factor. (The magic number 18420 is $1000 \cdot \ln 10^8$. We considered any prime between 814,279 (the largest ``small prime'' in Q) and 10⁸ to be of moderate size.) If

$$\left\lfloor 1000 \cdot \ln (c_1V-c_2T) \right\rfloor - C[c_2] < 10000,$$ (11)

the (c₁,c₂) pair was considered fully-smooth and assumed to have no prime factors outside of the factor base. (We need to use a large bound here as a result of an implementation decision concerning how real logarithms of c₁V-c₂T residues are calculated.) Although these assumptions were heuristic, they narrowed down the list of candidate (c₁,c₂) pairs so that time was not wasted on factoring c₁V - c₂T residues with multiple large prime factors.

Even though we had collected both partially- and fully-smooth pairs of (c₁,c₂) values, we ended up using only the fully-smooth pairs to generate equations. We found enough fully-smooth pairs in the search space to generate the needed equations; it was not necessary to resort to a ``large prime variation'' and allow equations with one or two large prime factors [27]. (The large prime variation would have increased the density of the set of equations which have to be solved, which was undesirable.) On average, for each value of c₁ we found 17.27 values of c₂ in the interval [1,300000] such that (c₁,c₂)was fully-smooth. Overall, we found about $1.6 \times 10^6$fully-smooth (c₁,c₂) pairs (out of $2.7 \times 10^{10}$ checked values), or about one fully-smooth pair per $1.65 \times 10^4$ possible pairs.

For each fully-smooth pair (c₁,c₂), the corresponding equation:

$$c_1V - c_2T \equiv V(c_1 + c_2s) \bmod p'$$

was factored (if possible) over the extended (complex) factor base. It was likely that c₁V-c₂Twas smooth with respect to the small real primes; it was not always the case, however, that the right-hand side of the equation could be factored over the small complex primes. Only one-third of the fully-smooth (c₁,c₂) pairs of integers actually yielded equations which could be completely factored over the complex factor base.

Once factored, the equations were sorted, and duplicate appearances of equations were removed. We found two causes of duplicate equations. The first cause was two distinct (c₁,c₂) pairs of integers which happened to have the same residue $c_1V-c_2T \bmod p$. The second (and much more frequent) cause was two pairs (c₁,c₂) and (c₁',c₂')with $c_1' = k \cdot c_1$, $c_2' = k \cdot c_2$, where k is a small integer. In this case,

$$\begin{eqnarray*}c_1'V - c_2'T = V(c_1' + c_2's) & \Longrightarrow & k(c_1V - c_2T) = kV(c_1 + c_2s) \end{eqnarray*}$$

and the factors of k cancel, leaving the same equation generated directly from the (c₁,c₂) pair. (Future implementations should probably only consider (c₁,c₂) pairs where c₁ and c₂ are relatively prime.) When all the duplicate equations were removed, we were left with a system of 288,017 distinct equations in 119,775 unknowns.

Before this system of equations was reduced using structured Gaussian elimination (see Section 6), it was first ``squeezed'' to eliminate dependent variables. Consider a real prime $q \in Q$ which may be factored into complex primes (a+bs) and (a-bs). All three of these variables could appear as unknowns in any of the 288,017 equations. Since it is sufficient to compute the discrete logarithm for any two of the three, we may reduce the density of the system by rewriting all appearances of one of the three variables in terms of the other two. For each triple of related primes q_i,(a_i+b_is),(a_i-b_is), the prime which appeared the fewest number of times in the system of equations was rewritten in terms of the other two. In most cases we replaced (a_i-b_is) with q_i(a_i+b_is)^-1; in the remaining instances (a_i+b_is) was replaced with q_i(a_i-b_is)^-1. This process reduced the number of unknown variables appearing in the system to 96,321.

The sieving program was written in C and utilized the portable multiprecision arithmetic package distributed by A. Lenstra and M. Manasse with their multiple polynomial quadratic sieve program. The sieve was run on a Silicon Graphics 4D-220 computer. It has 4 R3000 MIPS Computers, Inc., 25 MHz processors, each rated at about 18 mips, or 18 times the speed of a DEC VAX 11/780 computer (and about 1.3 times the speed of a DECstation 3100). This machine has 128 Mbytes of main memory. Since our program occupied about 20 Mbytes of memory, all four processors could sieve simultaneously over different c₁ intervals.

It took our implementation approximately one hour of CPU time (running on a single processor) to successfully sieve 1000 possible c₁ values, for a total running time of about 100 hours. Factoring the equations derived from fully-smooth (c₁,c₂) pairs took another 104 minutes, on average, per 1000 possible c₁ values. Neither of these running times is particularly impressive; performance could probably be improved significantly without too much effort. However, it should be noted that even our implementation was able to collect the necessary equations using only a few days worth of CPU time.

Aside from improving the efficiency of the algorithm, one could achieve about an order of magnitude improvement by working with a smaller factor base and using the single and double large prime variations. The reason we did not do this is that we wished to extrapolate our result to estimate the security of much larger systems. The efficiency of sieving is very easy to estimate, as one can run experiments on small ranges and also use the data from the large factoring experiments. On the other hand, this kind of extrapolation is harder to do in the case of solving linear equations, and so we decided to generate a very large system that we could experiment with.

In order to gain more information concerning the performance of the Gaussian integer scheme, experiments were also performed on a 224-bit prime q using essentially the same algorithm. This prime was chosen to be similar to the one used in the Sun system in that $\frac{q-1}{2}$is a prime and q = T² + 2V² for particular integers $T,V \lesssim \sqrt{q}$. The factor base was the same as that used in the case of the 192-bit prime. Sieving occurred over a substantially larger range; all values of c₁ in the interval $[-3 \times 10^5,9 \times 10^5]$ were considered. Our sieve size (the range of possible c₂ values) was kept at $3 \times 10^5$. Sieving took under 1200 hours of CPU time, and produced about $1.1 \times 10^6$ candidate fully-smooth (c₁,c₂)pairs. After removing those pairs for which c₁ and c₂ were not relatively prime, factoring both sides of the corresponding $c_1V - c_2T \equiv V(c_1+c_2s)$ equations, and ``squeezing'' the resulting system, we were left with a system of 164,841 equations in 94,398 unknowns.