1. primary decomposition

Any ideal $I \subset R$ in a Noetherian ring can be written as $I = \overset{r}{\underset{i=1}{\cap}} q_i$ with $q_i$ primary ideals (that is, $q_i \not= R$ and $gf \in q_i$ implies $g \in q_i$ or $f^p \in q_i$ for some $p > 0$).

This generalises the unique factorisation (valid in factorial rings) $f = f_1^{p_1} \cdot \ldots \cdot f_r^{p_r}$ with $f_i$ irreducible, from elements to ideals. In $K[x]$ we have both, unique factorisation and primary decomposition and any algorithm for primary decomposition needs factorisation (because a primary decomposition of a principal ideal $I = \langle f \rangle$ is equivalent to a factorisation of $f$).

In contrast to factorisation, primary decomposition is, in general, not unique, even if we consider minimal decompositions, that is, the associated primes $p_i = \sqrt{q_i}$ are all distinct and none of the $q_i$ can be omitted in the intersection. However, the minimal (or isolated) primes, that is, the minimal elements of Ass$\,(I) = \{p_1, \dots, p_r\}$ with regard to inclusion, are uniquely determined. The minimal primes are the only ``geometrically visible'' primes in the sense that

$$V(I) = \bigcup_{p_j \in \text{minAss}\,(I)} V(p_j)$$

is the decomposition of $V(I)$ into irreducible components. A non-minimal associated prime $p_i \not\in$   minAss$\,(I)$ is called embedded, because there exists a $p_j \in$   minAss$\,(I),\; p_j \subset p_i$. This means geometrically $V(p_i) \subset V(p_j)$, that is, the irreducible component of $V(I)$ corresponding to $p_i$ is embedded in some bigger irreducible component.

As an example we compute the primary decomposition of the ideal $I = \langle x^2y^3-x^3yz,\; y^2z - xz^2\rangle$ in SINGULAR, the output being slightly changed in order to save space.

LIB "primdec.lib";              //calling library for primary decomposition
ring R  = 0,(x,y,z),dp;
ideal I = x2y3-x3yz,y2z-xz2;
primdecGTZ(I);
==> [1]: [1]:               [2]: [1]:               [3]: [1]:
           _[1]=-y2+xz              _[1]=z2                 _[1]=z
         [2]:                       _[2]=y                  _[2]=x2
           _[1]=-y2+xz           [2]:                    [2]:
                                    _[1]=z                  _[1]=z
                                    _[2]=y                  _[2]=x

The result is a list of three pairs of ideals (for each pair, the first ideal is the primary component, the second ideal the corresponding prime component). The second prime component [2] : [2] is embedded in the first [1] : [2]. The first primary component [1] : [1] is already prime, the other two are not.

Hence, $I = (y^2 - xz) \cap (y,z^2) \cap (x^2,z)$ and we obtain:

$$V(I) = \{y^2 - xz = 0\} \cup \underset{\text{\small (embedded component)}}{\{y=z^2=0\}} \cup \{x^2 = z = 0\}$$

All known algorithms for primary decompositions in $K[x]$ are quite involved and use many different sub-algorithms from various parts of computer algebra, in particular Gröbner bases, resp. characteristic sets, and multivariate polynomial factorisation over some (algebraic or transcendental) extension of the field $K$. For an efficient implementation which can treat examples of interest in algebraic geometry, a lot of extra small additional algorithms have to be used. In particular one should use ``easy'' splitting as soon and as often as possible, see DGP.

In SINGULAR the algorithms of GTZ (which was the first practical and general primary decomposition algorithm), the recent algorithm of SY and some of the homological algebra algorithms for primary decomposition of EHV have been implemented. For detailed and improved versions of these algorithms, together with extensive comparisons, see DGP.

Here are some major ingredients for primary decomposition:

1. Reduction to zero-dimensional primary decomposition (GTZ);
   
   maximal independent sets;
   ideal quotient, saturation, intersection.

2. Zero-dimensional primary decomposition (GTZ);
   
   lexicographical Gröbner basis;
   factorisation of multivariate polynomials;
   generic change of variables;
   primitive element computation.

Related algorithms:

1. Computation of the radical;
   
   square-free part of univariate polynomials;
   find (random) regular sequences (EHV).

2. Computation of the equidimensional part (EHV);
   
   Ext-annihilators;
   ideal quotients, saturation and intersection.

To see how homological algebra comes into play, let us compute the equidimensional part of $V(I)$, that is, the union of all maximal dimensional components of $V(I)$, or, algebraically, the intersection of all minimal primes. Following EHV, we can calculate the equidimensional part of a variety via Ext-groups:

If $c =$   codim$\,_{K[x]}(I)$, then the equidimensional part of $I$ is the annihilator ideal of the module Ext $^c_{K[x]} (K[x]/I, K[x])$ by EHV.

For example, the equidimensional part of $V = \{xz = yz = 0\}$ is given by the ideal $\langle z \rangle =$   ann$\,\bigl($Ext$^1 (K[x,y,z]/\langle xz,yz\rangle, K[x,y,z])\bigr)$.

Using SINGULAR, we obtain this via:

LIB "homolog.lib";  
ring  r  = 0,(x,y,z),dp;
ideal I  = xz, yz;
module M = Ext_R(1,I);
quotient(M,freemodule(nrows(M)));
==> _[1] =z

Note that module M = Ext_R(i,I) computes a presentation matrix of Ext$^i(R/I,R)$. Hence, identifying a matrix with its column space in the free module of rank equal to the number of rows, Ext$^1(R/I,R) = R^n/M$ with $R^n =$ freemodule(nrows(M)) and, therefore, Ann$\,\bigl($Ext$^1(R/I,R)\bigr) = M : R^n =$ quotient(M,freemodule(nrows(M))).

Above, we used the procedure Ext_R(-,-) from homolog.lib. Below we show that the Ext groups can easily be computed directly in a system which offers free resolutions, respectively syzygies, transposition of matrices and presentations of sub-quotients of a free module (modulo in SINGULAR). Indeed, the Ext-annihilator can be computed more directly (and faster) without computing the Ext group itself:

Take a free resolution of $R/I$:

$$0 \longleftarrow R/I \longleftarrow R \longleftarrow R^{n_1} \longleftarrow \cdots.$$

Then consider the dual sequence:

$$0 \longrightarrow$$   Hom$$(R,R) \overset{d^0}{\longrightarrow }$$   Hom$$(R^{n_1}, R) \overset{d^1}{\longrightarrow } \cdots .$$

This leads to:

Ext$^i(R/I,R) =$   Ker$\,(d^i)/$Im$\,(d^{i-1})$    and    Ann$\,\bigl($Ext$^i(R/I,R)\bigr) =$   Im$\,(d^{i-1}) :$   Ker$\,(d^i)$.

The corresponding SINGULAR commands are:

int i = 1;
resolution L = res(I,i+1);
module Im    = transpose(L[i]);
module Ker   = syz(transpose(L[i+1]));
module ext   = modulo(Ker,Im);           //the Ext group
ideal ann    = quotient(Im,Ker);         //the Ext-annihilator

Since the resolution can be computed by iterated syzygy computation, this is a beautiful example of geometric use of syzygies. However, the algorithm is not at all obvious, but based on the non-trivial theorem of Eisenbud, Huneke and Vasconcelos.