I'm resurrecting this thread because moving forward with degree 6 polynomial selection in msieve requires understanding the multidimensional generalization of this problem.

Suppose you want a degree 6 GNFS algebraic polynomial for a large input number, for example RSA200. Suppose further that you manage to find a polynomial with decent size properties, for example

Code:

R(x) =
-733353745969457329647324759819358
+10146121641768070157 x
A(x) =
+196658068845691396854154694784799560081209151290611987559
-593991137294951464528877714492146951852356505880 x
-1848998357142768993289988278919373127654 x^2
+1599893759674383491380658833465 x^3
+2423179599042378391811 x^4
-873569620551 x^5
+180 x^6

(Nobody get excited, this was just 45 minutes work with a GPU)

The traditional way of optimizing the root properties of A(x) is to choose (exquisitely carefully) a rotation polynomial Rot(x) such that A(x) + Rot(x) * R(x) has an unusually large number of roots modulo powers of many different small primes p^k.

You do that using a sieve; for a single p^k, look at each of the p^k residue classes r_i. If Rot(x) is a constant C, then we want C such that

A(r_i) + C * R(r_i) = 0 mod p^k

If a solution exists, the set of such values lies on an arithmetic progression

R0 + j * p^k

where R0 = -A(r_i) * modinv(R(r_i), p^k) mod p^k.

So collect all of the arithmetic progressions for all the residue classes modulo the first few p^k, then sieve to find an C that maximizes a weighted sum of all the arithmetic progressions that it falls on, and yet is small enough so that the size of the new A(x) is not adversely affected.

If Rot(x) is non-constant, things are only a little harder. If Rot(x) has a linear term Y*x, then just do a separate sieve for each Y, and for each Y modify each arithmetic progression to use a start value of (R0 - Y * r_i) mod p^k for residue class r_i. With quadratic Rot(x), use (R0 - Y * r_i - Z * r_i^2) mod p^k.

The amount of work to do in the root sieve is a power of the skew S to use for the polynomial. For degree 6 polynomials the range of admissible C is a bit more than the ratio of the constant coefficients in the algebraic polynomials, so that a single line of admissible C is of size about 10^26 (or ~S^2.5). That's a gigantic sieve! Even worse, for quadratic rotations the admissible range of Y is ~S^1.5 and the range for Z is about S^.5. So to fully explore the available search space of rotation polynomials one has to examine S^(2.5+1.5+.5) = S^4.5 rotations. For the polynomial above, S is about 770e6 so this is more sieve lines than we can ever hope to exhaustively search, each of which is larger than we could hope to exhaustively search.

There are two ways to deal with this; we could restrict the range that is actually searched, but this should be a method of last resort. The other possibility is what I'm working on: moving away from lines and sieving over lattices. For degree 5 and linear rotation polynomials this works pretty well; the range of C is ~S^1.5 and the range of Y is ~S^.5, so even for very large inputs it's feasible to search all the lines but over a lattice for each line, and the latter just involves modifying the initialization for the sieve. That approach won't work for degree 6; we need to sieve over a 3-D lattice that is massively skewed.

Which brings me to my problem: for residue class r_i mod p^k, pick an integer triplet (i,j,k); then the corresponding rotation is

Code:

| C | | p^k -r_i -r_i^2 | | i | | R0 |
| Y | = | 0 1 0 | | j | + | 0 |
| Z | | 0 0 1 | | k | | 0 |

I need to pick a lattice to sieve over that is the composition of several of these r_i for several different p^k. So how does one compute a basis for the intersection of two shifted 3x3 lattices? This is analogous to finding a reduced basis for NFS lattice sieving, but in the 3x3 case rather than the 2x2 case.

If the shift was not an issue, one could use an orthonormal transformation to produce basis vectors for the intersection, but to make the shifts not matter I need to find a common (C,Y,Z) triplet that is on both lattices. I only know how to do that by brute force, which works fine for degree 5 but is intractable when the product of the p^k is going to be 10^20 or so.

Bonus question: The feasible region for rotations looks like an ellipsoid that is probably not aligned with the unit vectors in th (C,Y,Z) plane. If the above can be done, and I then want to use lattice reduction to find a better basis for sieving, how can I twist the basis around so that small values of (i,j,k) are aligned to the axes of rotation of the ellipsoid? Doing this would allow the maximum number of points in the feasible region to sieve over.