    # Bad Primes in Computational Algebraic Geometry

Computations over the rational numbers often suffer from intermediate coefficient swell. One solution to this problem is to apply the given algorithm modulo a number of primes and then lift the modular results to the rationals. This method is guaranteed to work if we use a sufficiently large set of good primes. In many applications, however, there is no efficient way of excluding bad primes. In this note, we describe a technique for rational reconstruction which will nevertheless return the correct result, provided the number of good primes in the selected set of primes is large enough. We give a number of illustrating examples which are implemented using the computer algebra system Singular and the programming language Julia. We discuss applications of our technique in computational algebraic geometry.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Many exact computations in computer algebra are carried out over the rationals and extensions thereof. Modular techniques are an important tool to improve the performance of such algorithms since intermediate coefficient growth is avoided and the resulting modular computations can be done in parallel. For this, we require that the algorithm under consideration is also applicable over finite fields and returns a deterministic result. The fundamental approach is then as follows: Compute the result modulo a number of primes. Then reconstruct the result over from the modular results.

###### Example 1.

To compute

 12+13=56

using modular methods, the first step is to apply the Chinese remainder isomorphism:

 Z/5 × Z/7 × Z/101 ≅ Z/3535\vskip6.0ptplus2.0ptminus2.0pt 12 ⟼ ( ¯¯¯3 , ¯¯¯4 , ¯¯¯¯¯¯51 ) + 13 ⟼ ( ¯¯¯2 , ¯¯¯5 , ¯¯¯¯¯¯34 ) ∥ ( ¯¯¯0 , ¯¯¯2 , ¯¯¯¯¯¯85 ) ⟼ ¯¯¯¯¯¯¯¯590

The second step is to reconstruct a rational number from .

## 2. Rational Reconstruction

###### Theorem 2.

 For every integer , the -Farey map

is injective.

There are efficient algorithms for computing preimages of the Farey map, see, for example, [8, Sec. 5].

###### Example 3.

We use the computer algebra system Singular  to compute the preimage of the Farey map the setting of Example 1:

ring r = 0, x, dp;

farey(590,3535);

5/6

The basic concept for modular computations is then as follows:

1. Compute the result over for distinct primes .

2. Use the Chinese remainder isomorphism

 Z/N≅Z/p1×…×Z/pr

to lift the modular results to where .

3. Compute the preimage of the lift with respect to the -Farey map.

4. Verify the correctness of the lift.

This will yield the correct result, provided is large enough (that is, the -result is contained in the domain of the -Farey map), and provided none of the is bad.

###### Definition 4.

A prime is called bad (with respect to a fixed algorithm and input) if the result over does not reduce modulo to the result over .

By convention, this includes the case where, modulo , the input is not defined or the algorithm in consideration is not applicable.

### 3.1. Bad primes in Gröbner basis computations

Consider a set of variables and a monomial ordering on the monomials in . For a set of polynomials , write for its set of lead monomials. For and prime, write for the image of in .

###### Theorem 5.

 Suppose with all primitive and homogeneous. Let be the reduced Gröbner basis of , the reduced Gröbner basis of , and a minimal strong Gröbner basis of . Then

does not divide any lead coefficient in .

###### Example 6.

Using Singular, we determine the bad primes for a Gröbner basis computation of the Jacobian ideal of a projective plane curve. We compute a minimial strong Gröbner basis over :

option("redSB");
ring R = integer,(x, y, z),lp;
poly f = x7y5 + x2yz9 + xz11 + y3z9;
ideal I = groebner(ideal(diff(f, x), diff(f, y), diff(f,z)));
13781115527868730344777310464613260 83521912290113517241074608876444 60 12 4 12 12 45349632 12 1473863040 12 22674816 12 3888 12 12 12 13608 12 108 54 6 2 27 3 1 4 2 2 1 216 1 2 3 1 540 12 108 27 3 1 9 3 1 1 1 1 1 7 1 5 1 1

The bad primes, that is, the primes with , are then the prime factors

 p=2,3,5,7,11,13,257,247072949,328838088993550682027

of the lead coefficients. In contrast, the lead coefficients of the Gröbner basis over involve only the prime factors , and hence not all bad primes. As shown by the following computation, is indeed a bad prime:

ring R0 = 0,(x, y, z),lp;

ring R1 = 257,(x, y, z),lp;

### 3.2. Classification of Bad Primes

Bad primes can be classified as follows, see

[3, Sec. 3] for details:

• [leftmargin=*]

• Type 1: The input modulo is not valid (this poses no problem).

• Type 2: There is a failure in the course of the algorithm (for example, a matrix may not be invertible modulo ; this wastes computation time if it happens).

• Type 3: A computable invariant with known expected value (for example, a Hilbert polynomial) has a wrong value in a modular computation (to detect this we have to do expensive tests for each prime, although the set of bad primes usually is finite, and hence bad primes rarely occur).

• Type 4: A computable invariant with unknown expected value (for example, the lead ideal in a Gröbner basis computation) is wrong (this can be handled by a majority vote, however we have to compute the invariant for each modular result and store the modular results).

• Type 5: otherwise.

The Type 5 case in fact occurs, as is shown by the following example. For an ideal and a prime define .

###### Example 7.

Consider the algorithm computing the radical of the Jacobian ideal for the curve

 I=\scalebox1$⟨x6+y6+7x5z+x3y2z−31x4z2−224x3z3+244x2z4+1632xz5+576z6⟩$.

Note that, with respect to the degree reverse lexicographic order, , that is, is not bad with respect to the input. The following computation in Singular first determines the minimal associated primes of and .

LIB "primdec.lib";
ring R0 = 0, (x, y, z), dp;
poly f = x6+y6+7x5z+x3y2z-31x4z2-224x3z3+244x2z4+1632xz5+576z6;
ideal U0 = radical(ideal(f, diff(f, x), diff(f, y), diff(f, z)));
minAssGTZ(U0);
: w=y         : w=y
:  w=x+6z           w=x-4z
ring R5 = 5, (x, y, z), dp;
poly f =imap(R0,f);
ideal U5 = radical(ideal(f, diff(f, x), diff(f, y), diff(f, z)));
minAssGTZ(U5);
: w=y         : w=y
:  w=x-z             w=x+z
minassGTZ(imap(R0,U0));
: w=y
: w=x+z

This shows that , but .

## 4. Error-Tolerant Reconstruction

Our goal is to reconstruct the -result from the modular result in the presence of bad primes. Our basic strategy will be to find an element with in the lattice

 Λ=⟨(N,0),(r,1)⟩⊂Z2.
###### Lemma 8.

[3, Lem. 4.2] All with are collinear.

Now suppose with . We assume that is the product of the good primes with correct result , and is the product of the bad primes with wrong result .

###### Theorem 9.

[3, Lem. 4.3] If

and

 abmodN′=s

then . So if , then (by Lemma 8)

 xy=ab for all (x,y)∈Λ % with (x2+y2)

and such vectors exist. Moreover, if

and is a shortest vector in , we also have .

Hence, if , the Gauss-Lagrange-Algorithm for finding a shortest vector gives independently of , provided . We use the programming language Julia111See http://julialang.org/., to illustrate the resulting algorithm.

function ErrorTolerantReconstruction(r::Integer, N::Integer)
a1 = [N, 0]
a2 = [r, 1]
while dot(a1, a1) > dot(a2, a2)
q = dot(a1, a2)//dot(a2, a2)
a1, a2 = a2, a1 - Integer(round(q))*a2
end
if dot(a1, a1) < N
return a1//a1
else
return false
end
end

The following table shows timings (in seconds), for and of bit-length , comparing the Julia-function with implementations in the Singular-kernel (optimized C/C++ code) and the current Singular-interpreter:

Singular-kernel Julia Singular-interpreter
0.001 0.005 0.055

Building on Julia as a fast mid-level language, a backwards-compatible just-in-time compiled Singular-interpreter is under development.

###### Example 10.

In the setting of Example 1, we obtain from by

julia> ErrorTolerantReconstruction(590, 3535)
5//6

which computes the sequence

 (3535,0) =6⋅(590,1)+(−5,−6), (590,1) =−48⋅(−5,−6)+(350,−287).
###### Example 11.

Now we introduce an error in the modular results:

Error tolerant reconstruction computes

 (3535,0) =1⋅(2711,1)+(824,−1), (2711,1) =3⋅(824,−1)+(239,4) (824,−1) =3⋅(239,4)+(107,−13) (239,4) =2⋅(107,−13)+(25,30) (107,−13) =1⋅(25,30)+(82,−43)

hence yields

 2530=5⋅55⋅6=56.

Note that

 (52+62)⋅5=305<707=7⋅101.

## 5. General Reconstruction Scheme for Commutative Algebra

For a given ideal , we want to compute some ideal (or module) associated to by a deterministic algorithm. We proceed along the following lines:

1. Over compute from for in a suitable finite set of primes.

2. Replace by a subset according to a majority vote on (see also [3, Rmk. 5.7]).

3. For compute the coefficient-wise CRT–lift to , identifying generators by their lead monomials.

4. Lift by error tolerant rational reconstruction to .

5. Test for some random extra prime .

6. Verify .

7. If the lift, test or verification fails, then enlarge and repeat.

###### Theorem 12.

[3, Lem. 5.6] If the bad primes form a Zariski closed proper subset of , then this strategy terminates with the correct result.

We discuss an application from algebraic geometry. The goal is to compute adjoint curves, that is, curves which pass with sufficiently high multiplicity through the singularities of a given curve, see Figure 1. We consider an integral, non-degenerate projective curve with normalization map , and a saturated homogeneous ideal with . We write for the singular locus of . Let

be the pullback of a hyperplane, and

the pullback of . Then the exact sequence

 0→˜IOΓ→π∗(˜IO¯¯¯Γ)→F→0

induces, for , an exact sequence

 0→Im/I(Γ)m¯¯¯¯¯¯ϱm→H0(¯¯¯¯Γ,O¯¯¯Γ(mH−Δ(I)))→H0(Γ,F)→0.
###### Definition 13.

The ideal is an adjoint ideal of if is surjective for .

Since , we obtain:

###### Theorem 14.

 With notation as above:

 I is an adjoint ideal of Γ ⟺ IP¯¯¯¯¯¯¯¯¯¯¯OΓ,P=IP for all P∈Sing(Γ).

The conductor of is the largest ideal of which is also an ideal in .

###### Definition 15.

The Gorenstein adjoint ideal of is the largest homogeneous ideal with

 GP=COΓ,P for all P∈Sing(Γ).

The Gorenstein adjoint ideal has many applications in the geometry of curves.

###### Example 16.

If be an irreducible plane algebraic curve of degree , then cuts out the canonical linear series.

###### Example 17.

If is a rational plane curve of degree , then maps to a rational normal curve of degree in .

###### Example 18.

The Gorenstein adjoint ideal can be used in the Brill-Noether-Algorithm to compute Riemann-Roch spaces for singular curves.

The Gorenstein adjoint ideal can be computed via a local-to-global strategy.

###### Definition 19.

The local adjoint ideal of at is the largest homogeneous ideal with

###### Lemma 20.

[5, Prop. 5.4] With notation as above,

 G=⋂P∈SingΓG(P).
###### Definition 21.

Let be the coordinate ring of an affine model of and let . A ring is called a minimal local contribution to at if and for all .

The minimal local contribution to at is unique and can be computed using Grauert-Remmert-type normalization algorithms, see . It can be written as with an ideal and a common denominator .

###### Algorithm 22.

[5, Alg. 4] With notation as above, is the homogenization of the preimage of under .

## 7. Modular version of the algorithm

Applying the general modular strategy gives an algorithm which is two-fold parallel (taking Lemma 20 into account). We use primes such that the algorithm is applicable to the variety defined by . Efficient verification can be realized through a semi-continuity argument, see [5, Theorem 8.14]. Table 1 gives timings (in seconds on a 2.2 GHz processor) for plane curves of degree with singularities of type .

Rows LA and IQ refer to global computations of the Gorenstein adjoint ideal via linear algebra  and ideal quotients, respectively. The row Maple-IB shows timings for the normalization of the curve via a computation of an integral basis in Maple . The row locIQ gives timings for the local-to-global (Lemma 20), and modLocIQ for the modular local-to-global strategy. We also give timings for parallel computations and for the modular probabilistic algorithm obtained by omitting the verification. In square brackets, the number of primes in the modular strategy is shown, in round brackets the number of cores used simultaneously in a parallel computation.

Observe that, in the example, a local-to-global strategy does not give any benefit when computing over the rationals, since the singular locus does not decompose. However, by Chebotarev’s density theorem, the singular locus is likely to decompose when passing to a finite field, as illustrated by the last two rows of the table.

## References

•  Arbarello, E.; Ciliberto, C.: Adjoint hypersurfaces to curves in following Petri, in Commutative Algebra, Lecture Notes in Pure and Applied Mathematics, vol. 84, Dekker, New York, 1-21 (1983).
•  Arnold, E. A.: Modular algorithms for computing Gröbner bases, J. Symb. Comput. 35, 403–419 (2003).
•  J. Böhm, W. Decker, C. Fieker, G. Pfister. The use of bad primes in rational reconstruction, Math. Comp. 84, 3013–3027 (2015).
•  J. Böhm, W. Decker, S. Laplagne, G. Pfister, A. Steenpaß, S. Steidel. Parallel algorithms for normalization, J. Symb. Comp. 51, 99–114 (2013).
•  J. Böhm, W. Decker, G. Pfister, S. Laplagne. Local to global algorithms for the Gorenstein adjoint ideal of a curve, Preprint (2015), arXiv:1505.05040.
•  Decker, W., Greuel, G.-M., Pfister, G., Schönemann, H., 2015. Singular 4-0-2 – A computer algebra system for polynomial computations. http://www.singular.uni-kl.de
•  G.-M. Greuel, S. Laplagne, S. Seelisch, Normalization of rings, J. Symb. Comp. 45(9), 887–901 (2010).
•  P. Kornerup, R. T. Gregory, Mapping integers and Hensel codes onto Farey fractions, BIT 23, 9–20 (1983).
•  Mnuk, M.: An algebraic approach to computing adjoint curves. J. Symbolic Comput., 23(2-3), 229-240 (1997).
•  van Hoeij, M.: An algorithm for computing an integral basis in an algebraic function field. J. Symbolic Comput. 18, no. 4, 353-363 (1994).