This is NoetherianOperators.info, produced by Macaulay2, version 1.24.11

NoetherianOperators -- algorithms for computing local dual spaces and sets of Noetherian operators
**************************************************************************************************

Description
===========

The NoetherianOperators package includes algorithms for computing Noetherian
operators and local dual spaces of polynomial ideals, and related local
combinatorial data about its scheme structure. In addition, the package
provides symbolic methods for computing Noetherian operators and Noetherian
multipliers of polynomial modules.

The problem of characterizing ideal membership with differential conditions was
first addressed by Gröbner ("Uber eine neue idealtheoretische Grundlegung der
algebraischen Geometrie", Math. Ann. 115 (1938), no. 1, 333–358). Despite this
early algebraic interest by Gröbner, a complete description of primary ideals
and modules in terms of differential operators was first obtained by analysts
in the Fundamental Principle of Ehrenpreis and Palamodov. At the core of the
Fundamental Principle, one has the notion of Noetherian operators to describe a
primary module.

In case of an ideal supported at one point a set of Noetherian operators forms
a Macaulay inverse system that spans the dual space of the ideal. These notions
relate to the work of Macaulay ("The algebraic theory of modular systems",
Cambridge Press, (1916)).

In this package, we implement several (exact symbolic and approximate
numerical) algorithms for the computation of sets of Noetherian operators.

Methods and types for computing and manipulating Noetherian operators:

  * DiffOp
  * noetherianOperators
      * Strategy => "PunctualQuot"
      * Strategy => "MacaulayMatrix"
      * Strategy => "Hybrid"
  * specializedNoetherianOperators
  * numericalNoetherianOperators
  * getIdealFromNoetherianOperators
  * getModuleFromNoetherianOperators
  * coordinateChangeOps
  * noethOpsFromComponents
  * solvePDE
  * differentialPrimaryDecomposition

Methods for numerically computing and manipulating local dual spaces:

  * truncatedDual
  * zeroDimensionalDual
  * eliminatingDual
  * localHilbertRegularity
  * gCorners
  * pairingMatrix
  * isPointEmbedded
  * isPointEmbeddedInCurve
  * colon

Auxiliary numerical linear algebra methods:

  * rationalInterpolation

For the task of computing Noetherian operators, here we implement the
algorithms developed in the papers Noetherian Operators and Primary
Decomposition (see https://arxiv.org/abs/2006.13881 ), Primary ideals and their
differential equations (see https://arxiv.org/abs/2001.04700 ), and Primary
decomposition of modules: a computational differential approach (see
https://arxiv.org/abs/2104.03385 ). These include both symbolic and numerical
algorithms, and a hybrid algorithm, where numerical data is used to speed up
the symbolic algorithm.

To compute the initial ideal and Hilbert regularity of positive dimensional
ideals we use the algorithm of R. Krone ("Numerical algorithms for dual bases
of positive-dimensional ideals." Journal of Algebra and Its Applications,
12(06):1350018, 2013.). These techniques are numerically stable, and can be
used with floating point arithmetic over the complex numbers. They provide a
viable alternative in this setting to purely symbolic methods such as standard
bases.

Authors
=======

  * Robert Krone <krone@math.gatech.edu>
  * Justin Chen <justin.chen@math.gatech.edu>
  * Marc Harkonen <harkonen@gatech.edu>
  * Yairon Cid-Ruiz <Yairon.CidRuiz@UGent.be>
  * Anton Leykin <anton.leykin@gmail.com>

Certification a gold star
=========================

Version 2.2.1 of this package was accepted for publication in volume 12 (see
https://msp.org/jsag/2022/12-1/ ) of The Journal of Software for Algebra and
Geometry (see https://msp.org/jsag/ ) on 26 September 2022, in the article
Noetherian operators in Macaulay2 (see https://msp.org/jsag/2022/12-1/p05.xhtml
) (DOI: 10.2140/jsag.2022.12.33 (see https://doi.org/10.2140/jsag.2022.12.33
)). That version can be obtained from the journal (see
https://msp.org/jsag/2022/12-1/jsag-v12-n1-x05-NoetherianOperators.m2 ) or from
the Macaulay2 source code repository (see
https://github.com/Macaulay2/M2/blob/28faaabb38111c186c23ef6e6a3d487b0823390e/
M2/Macaulay2/packages/NoetherianOperators.m2 ).

Version
=======

This documentation describes version 2.2.1 of NoetherianOperators.

Source code
===========

The source code from which this documentation is derived is in the file
NoetherianOperators.m2 (see common/share/Macaulay2/NoetherianOperators.m2 ).

Exports
=======

  * Types
      * DiffOp -- differential operator
  * Functions and commands
      * amult -- Computes the arithmetic multiplicity of a module
      * coordinateChangeOps -- induced Noetherian
        operators under coordinate change
      * differentialPrimaryDecomposition
        -- compute a differential primary
        decomposition
      * diffOp -- create a differential operator
      * diffOpRing -- create and cache the ring of
        differential operators
      * eliminatingDual -- eliminating dual space of a
        polynomial ideal
      * gCorners -- generators of the initial ideal of a
        polynomial ideal
      * getIdealFromNoetherianOperators,
        -- Computes a primary ideal corresponding to a list of Noetherian
        operators and a prime ideal
      * getModuleFromNoetherianOperators
        -- Computes a primary submodule
        corresponding to a list of Noetherian operators and a prime ideal
      * "isPointEmbedded" -- see
        isPointEmbedded(AbstractPoint,Ideal,List),
        -- numerically
        determine if the point is an embedded component of the scheme
      * "isPointEmbeddedInCurve" -- see
        isPointEmbeddedInCurve(AbstractPoint,Ideal),
        -- numerically
        determine if the point is an embedded component of a 1-dimensional
        scheme
      * joinIdeals -- Computes the join of two ideals
      * localHilbertRegularity -- regularity of
        the local Hilbert function of a polynomial ideal
      * mapToPunctualHilbertScheme -- maps
        an ideal into a point in a certain punctual Hilbert scheme
      * noetherianOperators -- Noetherian operators
      * noethOpsFromComponents -- merge
        Noetherian operators for non-primary ideals
      * normalize -- rescale a differential operator
      * numericalNoetherianOperators, --
        Noetherian operators via numerical interpolation
      * orthogonalInSubspace -- Orthogonal of a
        space
      * pairingMatrix -- Applies dual space functionals
        to polynomials
      * rationalInterpolation -- numerically
        interpolate rational functions
      * "solvePDE" -- see solvePDE(Module), --
        solve linear systems of PDE with constant coefficients
      * specializedNoetherianOperators,
        -- Noetherian operators evaluated at a point
      * truncatedDual -- truncated dual space of a
        polynomial ideal
      * zeroDimensionalDual -- dual space of a
        zero-dimensional polynomial ideal
  * Methods
      * "amult(Ideal)" -- see amult, -- Computes the arithmetic
        multiplicity of a module
      * "amult(Module)" -- see amult, -- Computes the arithmetic
        multiplicity of a module
      * "coordinateChangeOps(Matrix,DiffOp)" -- see coordinateChangeOps,
        -- induced Noetherian operators under coordinate
        change
      * "coordinateChangeOps(Matrix,List)" -- see coordinateChangeOps,
        -- induced Noetherian operators under coordinate
        change
      * "coordinateChangeOps(RingMap,DiffOp)" -- see coordinateChangeOps,
        -- induced Noetherian operators under coordinate
        change
      * "coordinateChangeOps(RingMap,List)" -- see coordinateChangeOps,
        -- induced Noetherian operators under coordinate
        change
      * "differentialPrimaryDecomposition(Ideal)" -- see
        differentialPrimaryDecomposition, --
        compute a differential primary decomposition
      * "differentialPrimaryDecomposition(Module)" -- see
        differentialPrimaryDecomposition, --
        compute a differential primary decomposition
      * DiffOp Matrix -- apply a differential operator
      * "DiffOp RingElement" -- see DiffOp Matrix, --
        apply a differential operator
      * diffOp(Matrix) -- create a differential
        operator
      * "diffOp(RingElement)" -- see diffOp(Matrix),
        -- create a differential operator
      * "eliminatingDual(AbstractPoint,Ideal,List,ZZ)" -- see
        eliminatingDual, -- eliminating dual space of a
        polynomial ideal
      * "eliminatingDual(AbstractPoint,Matrix,List,ZZ)" -- see
        eliminatingDual, -- eliminating dual space of a
        polynomial ideal
      * evaluate(DiffOp,AbstractPoint)
        -- evaluate coefficients of a
        differential operator
      * "evaluate(DiffOp,Matrix)" -- see evaluate(DiffOp,AbstractPoint),
        -- evaluate coefficients of a
        differential operator
      * "gCorners(AbstractPoint,Ideal)" -- see gCorners, --
        generators of the initial ideal of a polynomial ideal
      * "gCorners(AbstractPoint,Matrix)" -- see gCorners, --
        generators of the initial ideal of a polynomial ideal
      * "getIdealFromNoetherianOperators(List,Ideal)" -- see
        getIdealFromNoetherianOperators, --
        Computes a primary ideal corresponding to a list of Noetherian
        operators and a prime ideal
      * "getModuleFromNoetherianOperators(Ideal,List)" -- see
        getModuleFromNoetherianOperators, --
        Computes a primary submodule corresponding to a list of Noetherian
        operators and a prime ideal
      * "hilbertFunction(DualSpace)" -- see
        hilbertFunction(ZZ,DualSpace),
      * "hilbertFunction(List,DualSpace)" -- see
        hilbertFunction(ZZ,DualSpace),
      * hilbertFunction(ZZ,DualSpace)
      * isPointEmbedded(AbstractPoint,Ideal,List)
        -- numerically
        determine if the point is an embedded component of the scheme
      * isPointEmbeddedInCurve(AbstractPoint,Ideal)
        -- numerically
        determine if the point is an embedded component of a 1-dimensional
        scheme
      * "joinIdeals(Ideal,Ideal)" -- see joinIdeals, --
        Computes the join of two ideals
      * "localHilbertRegularity(AbstractPoint,Ideal)" -- see
        localHilbertRegularity, -- regularity of the
        local Hilbert function of a polynomial ideal
      * "localHilbertRegularity(AbstractPoint,Matrix)" -- see
        localHilbertRegularity, -- regularity of the
        local Hilbert function of a polynomial ideal
      * "mapToPunctualHilbertScheme(Ideal)" -- see
        mapToPunctualHilbertScheme, -- maps an
        ideal into a point in a certain punctual Hilbert scheme
      * noetherianOperators(Ideal), --
        Noetherian operators of a primary ideal
      * noetherianOperators(Ideal,Ideal)
        -- Noetherian operators of a
        primary component
      * noetherianOperators(Module), --
        Noetherian operators of a primary submodule
      * noetherianOperators(Module,Ideal)
        -- Noetherian operators of a
        primary component
      * "noethOpsFromComponents(List)" -- see noethOpsFromComponents,
        -- merge Noetherian operators for non-primary
        ideals
      * "normalize(DiffOp)" -- see normalize, -- rescale a
        differential operator
      * "numericalNoetherianOperators(Ideal)" -- see
        numericalNoetherianOperators, --
        Noetherian operators via numerical interpolation
      * "orthogonalInSubspace(DualSpace,PolySpace,Number)" -- see
        orthogonalInSubspace, -- Orthogonal of a space
      * "orthogonalInSubspace(PolySpace,PolySpace,Number)" -- see
        orthogonalInSubspace, -- Orthogonal of a space
      * "pairingMatrix(PolySpace,DualSpace)" -- see pairingMatrix,
        -- Applies dual space functionals to polynomials
      * "pairingMatrix(PolySpace,PolySpace)" -- see pairingMatrix,
        -- Applies dual space functionals to polynomials
      * "pairingMatrix(PolySpace,RingElement)" -- see pairingMatrix,
        -- Applies dual space functionals to polynomials
      * "pairingMatrix(RingElement,DualSpace)" -- see pairingMatrix,
        -- Applies dual space functionals to polynomials
      * "pairingMatrix(RingElement,RingElement)" -- see pairingMatrix,
        -- Applies dual space functionals to polynomials
      * "rationalInterpolation(List,List,Matrix)" -- see
        rationalInterpolation, -- numerically
        interpolate rational functions
      * "rationalInterpolation(List,List,Matrix,Matrix)" -- see
        rationalInterpolation, -- numerically
        interpolate rational functions
      * rationalInterpolation(List,List,Ring)
        -- numerically
        interpolate rational functions
      * "solvePDE(Ideal)" -- see solvePDE(Module),
        -- solve linear systems of PDE with constant coefficients
      * "solvePDE(Matrix)" -- see solvePDE(Module),
        -- solve linear systems of PDE with constant coefficients
      * solvePDE(Module) -- solve linear systems
        of PDE with constant coefficients
      * "specializedNoetherianOperators(Ideal,AbstractPoint)" -- see
        specializedNoetherianOperators, --
        Noetherian operators evaluated at a point
      * "specializedNoetherianOperators(Ideal,Matrix)" -- see
        specializedNoetherianOperators, --
        Noetherian operators evaluated at a point
      * truncate(DualSpace,List,ZZ),
        -- truncate a polynomial space or dual space
      * "truncate(PolySpace,List,ZZ)" -- see truncate(DualSpace,List,ZZ),
        -- truncate a polynomial space or
        dual space
      * truncate(DualSpace,ZZ) -- truncate
        a polynomial space or dual space
      * "truncate(PolySpace,ZZ)" -- see truncate(DualSpace,ZZ),
        -- truncate a polynomial space or dual
        space
      * "truncatedDual(AbstractPoint,Ideal,ZZ)" -- see truncatedDual,
        -- truncated dual space of a polynomial ideal
      * "truncatedDual(AbstractPoint,Matrix,ZZ)" -- see truncatedDual,
        -- truncated dual space of a polynomial ideal
      * "truncatedDual(Matrix,Ideal,ZZ)" -- see truncatedDual,
        -- truncated dual space of a polynomial ideal
      * "truncatedDual(Matrix,Matrix,ZZ)" -- see truncatedDual,
        -- truncated dual space of a polynomial ideal
      * "zeroDimensionalDual(AbstractPoint,Ideal)" -- see
        zeroDimensionalDual, -- dual space of a
        zero-dimensional polynomial ideal
      * "zeroDimensionalDual(AbstractPoint,Matrix)" -- see
        zeroDimensionalDual, -- dual space of a
        zero-dimensional polynomial ideal
      * "zeroDimensionalDual(Matrix,Ideal)" -- see zeroDimensionalDual,
        -- dual space of a zero-dimensional polynomial
        ideal
      * "zeroDimensionalDual(Matrix,Matrix)" -- see zeroDimensionalDual,
        -- dual space of a zero-dimensional polynomial
        ideal
  * Symbols
      * DependentSet -- option for computing Noetherian
        operators
      * "StandardBasis" -- see gCorners, -- generators of the
        initial ideal of a polynomial ideal
      * "AllVisible" -- see isPointEmbedded(AbstractPoint,Ideal,List),
        -- numerically
        determine if the point is an embedded component of the scheme
      * "Regularity" -- see isPointEmbeddedInCurve(AbstractPoint,Ideal),
        -- numerically
        determine if the point is an embedded component of a 1-dimensional
        scheme
      * "Rational" -- see noetherianOperators(Ideal,Ideal),
        -- Noetherian operators of a
        primary component
      * "InterpolationDegreeLimit" -- see numericalNoetherianOperators,
        -- Noetherian operators via numerical
        interpolation
      * "InterpolationTolerance" -- see numericalNoetherianOperators,
        -- Noetherian operators via numerical
        interpolation
      * "NoetherianDegreeLimit" -- see numericalNoetherianOperators,
        -- Noetherian operators via numerical
        interpolation
      * "TrustedPoint" -- see numericalNoetherianOperators,
        -- Noetherian operators via numerical
        interpolation
      * Sampler -- optional sampler function
      * "IntegralStrategy" -- see Strategy => "MacaulayMatrix",
        -- strategy for computing Noetherian operators
      * "KernelStrategy" -- see Strategy => "MacaulayMatrix",
        -- strategy for computing Noetherian operators

For the programmer
==================

The object NoetherianOperators, is a package (Macaulay2Doc).

amult -- Computes the arithmetic multiplicity of a module
*********************************************************

Synopsis
========

  * Usage:
        m = amult(M)
        m = amult(I)
  * Inputs:
      * M, a module (Macaulay2Doc), a submodule of a free module
      * I, an ideal (Macaulay2Doc), an ideal
  * Outputs:
      * m, an integer (Macaulay2Doc), the arithmetic multiplicity of
        the module or the ideal

Description
===========

Given a submodule $M$ of a free module $F$, one computes the arithmetic
multiplicity of $M$ as the sum, along the associated primes of $F/M$, of the
length of the largest submodule of finite length of the quotient $M/F$
localized at the associated prime. The arithmetic multiplicity is a fundamental
invariant from a differential point of view as it yields the minimal size of a
differential primary decomposition. For more details the reader is referred to
the paper Primary Decomposition with Differential Operators (see
https://arxiv.org/abs/2101.03643 ).

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x1,x2,x3,x4]                                                     |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : U = image matrix{{x1*x2,x2*x3,x3*x4,x4*x1}, {x1^2,x2^2,x3^2,x4^2}}      |
|                                                                             |
|o2 = image | x1x2 x2x3 x3x4 x1x4 |                                           |
|           | x1^2 x2^2 x3^2 x4^2 |                                           |
|                                                                             |
|                             2                                               |
|o2 : R-module, submodule of R                                                |
+-----------------------------------------------------------------------------+
|i3 : amult U                                                                 |
|                                                                             |
|o3 = 22                                                                      |
+-----------------------------------------------------------------------------+
|i4 : I = ideal( x1^3*x3^2-x2^5, x2^2*x4^3-x3^5, x1^5*x4^2-x2^7, x1^2*x4^5-x3^|
|                                                                             |
|                5     3  2      5     2  3      7     5  2      7     2  5   |
|o4 = ideal (- x2  + x1 x3 , - x3  + x2 x4 , - x2  + x1 x4 , - x3  + x1 x4 )  |
|                                                                             |
|o4 : Ideal of R                                                              |
|-----------------------------------------------------------------------------|
|7 )                                                                          |
+-----------------------------------------------------------------------------+
|i5 : amult I                                                                 |
|                                                                             |
|o5 = 207                                                                     |
+-----------------------------------------------------------------------------+

Ways to use amult:
==================

  * "amult(Ideal)"
  * "amult(Module)"

For the programmer
==================

The object amult, is a method function (Macaulay2Doc).

colon -- colon of a (truncated) dual space
******************************************

Synopsis
========

  * Usage:
        Dg = colon(D,g)
        DJ = colon(D,J)
  * Inputs:
      * D, an instance of the type DualSpace (NAGtypes),
      * g, a ring element (Macaulay2Doc),
      * J, an ideal (Macaulay2Doc),
  * Optional inputs (Macaulay2Doc),:
      * Tolerance => ..., default value .000001
  * Outputs:
      * an instance of the type DualSpace (NAGtypes), Dg, DJ

Description
===========

This computes (a part of) the dual space of the dual. See

R. Krone and A. Leykin, "Numerical algorithms for detecting embedded
components.", arXiv:1405.7871

for a description.

For the programmer
==================

The object colon, is a method function with options (Macaulay2Doc).

coordinateChangeOps -- induced Noetherian operators under coordinate change
***************************************************************************

Synopsis
========

  * Usage:
        coordinateChangeOps(phi, D)
        coordinateChangeOps(phi, L)
  * Inputs:
      * phi, a matrix (Macaulay2Doc), or a ring map (Macaulay2Doc),
      * D, a differential operator,
      * L, a list (Macaulay2Doc), of differential operators,
  * Outputs:
      * a differential operator, (resp. a list of differential
        operators)

Description
===========

Let $I$ be an ideal in a polynomial ring $K[x_1, ..., x_n]$, and $\phi \in
GL_n(K)$ a matrix representing a $K$-linear automorphism of $R$. Then there is
an automorphism $\psi$ of the Weyl algebra $K[x_i, dx_i]$ such that if $D_1,
..., D_r$ is a set of Noetherian operators for $I$ then $\psi(D_1), ...,
\psi(D_r)$ is a set of Noetherian operators for $\phi(I)$. This function
computes the induced operators for a given $\phi$. The action of $\psi$ on
polynomial variables $x_i$ is given by $\phi$, while the action of $\psi$ on
differential variables $dx_i$ is given by the inverse transpose of $\phi$.

+---------------------------------------------------------------------+
|i1 : R = QQ[x,y,t]                                                   |
|                                                                     |
|o1 = R                                                               |
|                                                                     |
|o1 : PolynomialRing                                                  |
+---------------------------------------------------------------------+
|i2 : I = ideal(x^2, y^2 - x*t)                                       |
|                                                                     |
|             2   2                                                   |
|o2 = ideal (x , y  - x*t)                                            |
|                                                                     |
|o2 : Ideal of R                                                      |
+---------------------------------------------------------------------+
|i3 : P = radical I                                                   |
|                                                                     |
|o3 = ideal (y, x)                                                    |
|                                                                     |
|o3 : Ideal of R                                                      |
+---------------------------------------------------------------------+
|i4 : N = noetherianOperators I                                       |
|                                                                     |
|o4 = {| 1 |, | dy |, | tdy^2+2dx |, | tdy^3+6dxdy |}                 |
|                                                                     |
|o4 : List                                                            |
+---------------------------------------------------------------------+
|i5 : phi = map(R, R, diagonalMatrix apply(numgens R, i -> random QQ))|
|                                                                     |
|                 9   1   9                                           |
|o5 = map (R, R, {-x, -y, -t})                                        |
|                 2   2   4                                           |
|                                                                     |
|o5 : RingMap R <-- R                                                 |
+---------------------------------------------------------------------+
|i6 : N' = coordinateChangeOps_phi N                                  |
|                                                                     |
|o6 = {| 1 |, | 2dy |, | 4tdy^2+4/9dx |, | 8tdy^3+8/3dxdy |}          |
|                                                                     |
|o6 : List                                                            |
+---------------------------------------------------------------------+
|i7 : I' = phi I                                                      |
|                                                                     |
|            81 2  1 2   81                                           |
|o7 = ideal (--x , -y  - --x*t)                                       |
|             4    4      8                                           |
|                                                                     |
|o7 : Ideal of R                                                      |
+---------------------------------------------------------------------+
|i8 : P' = phi P                                                      |
|                                                                     |
|            1   9                                                    |
|o8 = ideal (-y, -x)                                                  |
|            2   2                                                    |
|                                                                     |
|o8 : Ideal of R                                                      |
+---------------------------------------------------------------------+
|i9 : I' == getIdealFromNoetherianOperators(N', P')                   |
|                                                                     |
|o9 = false                                                           |
+---------------------------------------------------------------------+

See also
========

  * noetherianOperators -- Noetherian operators

Ways to use coordinateChangeOps:
================================

  * "coordinateChangeOps(Matrix,DiffOp)"
  * "coordinateChangeOps(Matrix,List)"
  * "coordinateChangeOps(RingMap,DiffOp)"
  * "coordinateChangeOps(RingMap,List)"

For the programmer
==================

The object coordinateChangeOps, is a method function (Macaulay2Doc).

DependentSet -- option for computing Noetherian operators
*********************************************************

Description
===========

Let $Q \subseteq R := \mathbb{K}[x_1,\dots,x_n]$ be an $d$-dimensional primary
ideal. Then there exists a set of $d$ variables in $R$ which is algebraically
independent in $R/I$. We refer to these as the independent variables, and the
remaining variables are the dependent variables. The function
independentSets (Macaulay2Doc), can compute sets of independent
variables for symbolic ideals.

The functions computing Noetherian operators, namely

  * noetherianOperators -- Noetherian operators
  * specializedNoetherianOperators, --
    Noetherian operators evaluated at a point
  * numericalNoetherianOperators, --
    Noetherian operators via numerical interpolation

pass to a polynomial ring in the dependent variables, with the coefficient
field being the fraction field of a polynomial ring in the independent
variables. Because of this, computing Noetherian operators requires a knowledge
of a dependent set of variables, which can be set using the option
DependentSet. Note that the $dx$-monomials will only involve dependent
variables.

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x,y];                                                            |
+-----------------------------------------------------------------------------+
|i2 : I = ideal((x+y)^2);                                                     |
|                                                                             |
|o2 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i3 : P = radical I;                                                          |
|                                                                             |
|o3 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i4 : A = noetherianOperators(I, P, DependentSet => {x})                      |
|                                                                             |
|o4 = {| 1 |, | dx |}                                                         |
|                                                                             |
|o4 : List                                                                    |
+-----------------------------------------------------------------------------+
|i5 : B = noetherianOperators(I, P, DependentSet => {y})                      |
|                                                                             |
|o5 = {| 1 |, | dy |}                                                         |
|                                                                             |
|o5 : List                                                                    |
+-----------------------------------------------------------------------------+
|i6 : getIdealFromNoetherianOperators(A, P) ==                                |
|                                                                             |
|o6 = true                                                                    |
|-----------------------------------------------------------------------------|
|getIdealFromNoetherianOperators(B, P)                                        |
+-----------------------------------------------------------------------------+

The symbolic method noetherianOperators, will
usually be able to figure out a dependent set of variables automatically. On
the other hand, numerical computations using
specializedNoetherianOperators, and
numericalNoetherianOperators, will usually
require the user to set the option DependentSet.

Caveat
======

The option DependentSet is ignored when calling noetherianOperators,
with Strategy => "PunctualQuot".
Note that this is the default strategy for
noetherianOperators(Ideal).

See also
========

  * noetherianOperators -- Noetherian operators
  * specializedNoetherianOperators, --
    Noetherian operators evaluated at a point
  * numericalNoetherianOperators, --
    Noetherian operators via numerical interpolation

For the programmer
==================

The object DependentSet, is a symbol (Macaulay2Doc).

differentialPrimaryDecomposition -- compute a differential primary decomposition
********************************************************************************

Synopsis
========

  * Usage:
        differentialPrimaryDecomposition U
  * Inputs:
      * U, a module (Macaulay2Doc), a submodule of a free module,
        or an ideal (Macaulay2Doc),
  * Outputs:
      * a list (Macaulay2Doc), a minimal differential primary
        decomposition of U

Description
===========

Let $R$ be a polynomial ring over a field $K$. Given a submodule $U$ of an
$R$-module $M$, a differential primary decomposition of $U$ in $M$ is a list of
pairs $(p_1, A_1), ..., (p_k, A_k)$ where $p_1, ..., p_k$ are the associated
primes of $M/U$ and $A_i \subseteq \operatorname{Diff}_{R/K}(M, R/p_i)$ are
differential operators satisfying $$U_p = \bigcap_{p_i \subseteq p} \{ w \in
M_p : \delta(w) = 0 , \ \forall \delta \in A_i \}.$$ This notion was introduced
in [2] (cf. Definition 4.1), in which it was shown that the size of a
differential primary decomposition (which is defined to be $\sum_{i=1}^k
|A_i|$) is at least amult(U), and moreover differential primary
decompositions of size equal to amult(U) exist (and are called minimal).

This method contains an implementation of Algorithm 4.6 in [2].

The following example appears as Example 6.2 in [1]:

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x_1,x_2,x_3]                                                     |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : U = image matrix {{x_1^2,x_1*x_2,x_1*x_3}, {x_2^2,x_2*x_3,x_3^2}}       |
|                                                                             |
|o2 = image | x_1^2 x_1x_2 x_1x_3 |                                           |
|           | x_2^2 x_2x_3 x_3^2  |                                           |
|                                                                             |
|                             2                                               |
|o2 : R-module, submodule of R                                                |
+-----------------------------------------------------------------------------+
|i3 : differentialPrimaryDecomposition U                                      |
|                                                                             |
|                                   2                                         |
|o3 = {{ideal x , {| 1 |}}, {ideal(x  - x x ), {| -x_3 |}}, {ideal (x , x ),  |
|              1   | 0 |            2    1 3    |  x_1 |             3   2    |
|     ------------------------------------------------------------------------|
|     {|   0  |}}}                                                            |
|      | dx_3 |                                                               |
|                                                                             |
|o3 : List                                                                    |
+-----------------------------------------------------------------------------+

References
==========

[1] Chen, J., Cid-Ruiz, Y. (2021). Primary decomposition of modules: a
computational differential approach (see https://arxiv.org/abs/2104.03385 )

[2] Cid-Ruiz, Y., Sturmfels, B. (2021). Primary Decomposition with Differential
Operators (see https://arxiv.org/abs/2101.03643 )

See also
========

  * solvePDE -- solve linear systems of PDE with
    constant coefficients
  * noetherianOperators -- Noetherian operators

Ways to use differentialPrimaryDecomposition:
=============================================

  * "differentialPrimaryDecomposition(Ideal)"
  * "differentialPrimaryDecomposition(Module)"

For the programmer
==================

The object differentialPrimaryDecomposition,
is a method function with options (Macaulay2Doc).

DiffOp -- differential operator
*******************************

Description
===========

A differential operator on the ring $R = \mathbb{K}[x_1,\dots,x_n]$ can be
thought of as $k$-vectors of polynomials in $S = R[dx_1, \dotsc, dx_n]$, with
coefficients in $R$, and monomials in variables $dx_1, \dots, dx_n$, where
$dx_i$ corresponds to the partial derivative with respect to $x_i$. Hence a
differential operator is an element of the free module $S^k$. These operators
form an $R$-vector space, and operate on elements of $R^k$. The result of the
operation lies in $R$, and is equal to the sum of the entrywise operations.

The ring $S$ can be obtained from $R$ using diffOpRing.

+---------------------------------------------+
|i1 : R = QQ[x,y]                             |
|                                             |
|o1 = R                                       |
|                                             |
|o1 : PolynomialRing                          |
+---------------------------------------------+
|i2 : S = diffOpRing R                        |
|                                             |
|o2 = S                                       |
|                                             |
|o2 : PolynomialRing                          |
+---------------------------------------------+
|i3 : D = diffOp((x+y)*dx + (3+x) * dx*dy^2)  |
|                                             |
|o3 = | (x+3)dxdy^2+(x+y)dx |                 |
|                                             |
|                1                            |
|o3 : DiffOp in S                             |
+---------------------------------------------+
|i4 : (x^2+3) * D                             |
|                                             |
|o4 = | (x3+3x2+3x+9)dxdy^2+(x3+x2y+3x+3y)dx ||
|                                             |
|                1                            |
|o4 : DiffOp in S                             |
+---------------------------------------------+
|i5 : D + D                                   |
|                                             |
|o5 = | (2x+6)dxdy^2+(2x+2y)dx |              |
|                                             |
|                1                            |
|o5 : DiffOp in S                             |
+---------------------------------------------+
|i6 : D(x^5*y^2)                              |
|                                             |
|       5 2     4 3      5      4             |
|o6 = 5x y  + 5x y  + 10x  + 30x              |
|                                             |
|o6 : R                                       |
+---------------------------------------------+
+---------------------------------------+
|i7 : D = diffOp(matrix{{x*dx}, {y*dy}})|
|                                       |
|o7 = | xdx |                           |
|     | ydy |                           |
|                                       |
|                2                      |
|o7 : DiffOp in S                       |
+---------------------------------------+
|i8 : f = matrix{{x^2}, {y^2}}          |
|                                       |
|o8 = | x2 |                            |
|     | y2 |                            |
|                                       |
|             2      1                  |
|o8 : Matrix R  <-- R                   |
+---------------------------------------+
|i9 : D f                               |
|                                       |
|       2     2                         |
|o9 = 2x  + 2y                          |
|                                       |
|o9 : R                                 |
+---------------------------------------+

See also
========

  * diffOp(Matrix) -- create a differential operator
  * diffOp(RingElement) -- create a differential
    operator

Functions and methods returning a differential operator:
========================================================

  * "coordinateChangeOps(Matrix,DiffOp)" -- see coordinateChangeOps,
    -- induced Noetherian operators under coordinate
    change
  * "coordinateChangeOps(RingMap,DiffOp)" -- see coordinateChangeOps,
    -- induced Noetherian operators under coordinate
    change
  * "normalize(DiffOp)" -- see normalize, -- rescale a
    differential operator

Methods that use a differential operator:
=========================================

  * DiffOp Matrix -- apply a differential operator
  * "DiffOp RingElement" -- see DiffOp Matrix, -- apply a
    differential operator
  * evaluate(DiffOp,AbstractPoint),
    -- evaluate coefficients of a differential operator
  * "evaluate(DiffOp,Matrix)" -- see evaluate(DiffOp,AbstractPoint),
    -- evaluate coefficients of a
    differential operator

For the programmer
==================

The object DiffOp, is a type (Macaulay2Doc), with
ancestor classes Vector (Macaulay2Doc), < BasicList (Macaulay2Doc)
< Thing (Macaulay2Doc).

diffOp -- create a differential operator
****************************************

Ways to use diffOp:
===================

  * diffOp(Matrix) -- create a differential operator
  * "diffOp(RingElement)" -- see diffOp(Matrix), --
    create a differential operator

For the programmer
==================

The object diffOp, is a method function (Macaulay2Doc).

DiffOp Matrix -- apply a differential operator
**********************************************

Synopsis
========

  * Operator: SPACE (Macaulay2Doc),
  * Usage:
        D f
  * Inputs:
      * D, a differential operator,
      * f, a matrix (Macaulay2Doc), consisting of one column,
        with entries in the coefficient ring of ring D, or a ring element (Macaulay2Doc).
  * Outputs:
      * a ring element (Macaulay2Doc), in the same ring as f

Description
===========

Let $R = \mathbb{F}[x_1,\dots,x_n]$ and $S = R[dx_1,\dotsc,dx_n]$. The elements
of $S$ operate naturally on elements of $R$. The operator $dx_i$ acts as a
partial derivative with respect to $x_i$, i.e., $dx_i \bullet f =
\frac{\partial f}{\partial x_i}$, and a polynomial acts by multiplication,
i.e., $x_i \bullet f = x_i f$.

Suppose $D \in S^k$ and $f \in R^k$. Then the operation of $D$ on $f$ is
defined as $D\bullet f := \sum_{i=1}^k D_i \bullet f_i \in R$.

+---------------------------------------------+
|i1 : R = QQ[x,y]                             |
|                                             |
|o1 = R                                       |
|                                             |
|o1 : PolynomialRing                          |
+---------------------------------------------+
|i2 : S = diffOpRing R                        |
|                                             |
|o2 = S                                       |
|                                             |
|o2 : PolynomialRing                          |
+---------------------------------------------+
|i3 : D = diffOp matrix{{x*dx}, {(y+1)*dx*dy}}|
|                                             |
|o3 = |    xdx    |                           |
|     | (y+1)dxdy |                           |
|                                             |
|                2                            |
|o3 : DiffOp in S                             |
+---------------------------------------------+
|i4 : f = matrix{{x+y}, {x*y*(y+1)}}          |
|                                             |
|o4 = | x+y    |                              |
|     | xy2+xy |                              |
|                                             |
|             2      1                        |
|o4 : Matrix R  <-- R                         |
+---------------------------------------------+
|i5 : D f                                     |
|                                             |
|       2                                     |
|o5 = 2y  + x + 3y + 1                        |
|                                             |
|o5 : R                                       |
+---------------------------------------------+

As with diffOp(Matrix), a $1\times 1$ matrix may be
replaced by a ring element (Macaulay2Doc).

+--------------------+
|i6 : D = diffOp dx^2|
|                    |
|o6 = | dx^2 |       |
|                    |
|                1   |
|o6 : DiffOp in S    |
+--------------------+
|i7 : D(x^3+y*x^2)   |
|                    |
|o7 = 6x + 2y        |
|                    |
|o7 : R              |
+--------------------+

Ways to use this method:
========================

  * DiffOp Matrix -- apply a differential operator
  * "DiffOp RingElement"

diffOp(Matrix) -- create a differential operator
************************************************

Synopsis
========

  * Function: diffOp,
  * Usage:
        diffOp M
  * Inputs:
      * M, a matrix (Macaulay2Doc), consisting of a single
        column, or a ring element (Macaulay2Doc).
  * Outputs:
      * a differential operator,

Description
===========

Creates a differential operator from a vector of polynomials in $S =
\mathbb{K}[x_1,\dotsc,x_n][dx_1,\dotsc,dx_n]$. The ring $S$ is obtained from
the ring $R = \mathbb{K}[x_1,\dotsc,x_n]$ by using diffOpRing.

+------------------------------------------------------------+
|i1 : R = QQ[x_1,x_2]                                        |
|                                                            |
|o1 = R                                                      |
|                                                            |
|o1 : PolynomialRing                                         |
+------------------------------------------------------------+
|i2 : S = diffOpRing R                                       |
|                                                            |
|o2 = S                                                      |
|                                                            |
|o2 : PolynomialRing                                         |
+------------------------------------------------------------+
|i3 : diffOp matrix {{(x_1 * x_2 + 3)*dx_1*dx_2^2}, {dx_2^2}}|
|                                                            |
|o3 = | (x_1x_2+3)dx_1dx_2^2 |                               |
|     |        dx_2^2        |                               |
|                                                            |
|                2                                           |
|o3 : DiffOp in S                                            |
+------------------------------------------------------------+

A ring element (Macaulay2Doc), can be used instead of a $1
\times 1$ matrix.

+--------------------------+
|i4 : diffOp (x_1^2*dx_1^2)|
|                          |
|o4 = | x_1^2dx_1^2 |      |
|                          |
|                1         |
|o4 : DiffOp in S          |
+--------------------------+

See also
========

  * DiffOp -- differential operator
  * diffOpRing -- create and cache the ring of differential
    operators

Ways to use this method:
========================

  * diffOp(Matrix) -- create a differential operator
  * "diffOp(RingElement)"

diffOpRing -- create and cache the ring of differential operators
*****************************************************************

Synopsis
========

  * Usage:
        diffOpRing R
  * Inputs:
      * R, a polynomial ring (Macaulay2Doc), in $n$
        variables
  * Outputs:
      * a polynomial ring (Macaulay2Doc), in $2n$
        variables
  * Consequences:
      * the ring is cached in {tt R} under the key "DiffOpRing".

Description
===========

Takes a polynomial ring $R = \mathbb{K}[x_1,\dotsc,x_n]$ and creates the ring
$S = R[dx_1,\dotsc,dx_n]$.

+----------------------------+
|i1 : R = QQ[x_1..x_3, a,b]; |
+----------------------------+
|i2 : S = diffOpRing R;      |
+----------------------------+
|i3 : gens S                 |
|                            |
|o3 = {dx , dx , dx , da, db}|
|        1    2    3         |
|                            |
|o3 : List                   |
+----------------------------+
|i4 : coefficientRing S      |
|                            |
|o4 = R                      |
|                            |
|o4 : PolynomialRing         |
+----------------------------+

Differential operators on $R$ have entries in $S$.

+------------------------------+
|i5 : ring diffOp(dx_3^2) === S|
|                              |
|o5 = true                     |
+------------------------------+
|i6 : ring diffOp(a_R) === S   |
|                              |
|o6 = true                     |
+------------------------------+

Subsequent calls to diffOpRing will not create new rings

+-----------------------+
|i7 : diffOpRing R === S|
|                       |
|o7 = true              |
+-----------------------+

Caveat
======

the created ring is not a Weyl algebra, it is a commutative ring

See also
========

  * diffOp(Matrix) -- create a differential operator
  * diffOp(RingElement) -- create a differential
    operator

For the programmer
==================

The object diffOpRing, is an a cache function (Macaulay2Doc).

eliminatingDual -- eliminating dual space of a polynomial ideal
***************************************************************

Synopsis
========

  * Usage:
        S = eliminatingDual(p, I, v, d)
  * Inputs:
      * p, a point (NAGtypes),
      * I, an ideal (Macaulay2Doc), or a one-row matrix (Macaulay2Doc),
        of generators
      * v, a list (Macaulay2Doc), a list of the integers
        designating which variables to bound
      * d, an integer (Macaulay2Doc), the degree bound for the
        designated variables
  * Optional inputs (Macaulay2Doc),:
      * Tolerance, => ..., default
        value null, optional argument for numerical tolerance
  * Outputs:
      * S, an instance of the type DualSpace (NAGtypes),

Description
===========

Given a list of variable indices, this method computes a basis for all dual
elements orthogonal to I which have total degree in the variables on the list
bounded by d.

+---------------------------------------------+
|i1 : R = CC[x,y];                            |
+---------------------------------------------+
|i2 : I = ideal{x^2-y^3}                      |
|                                             |
|              3    2                         |
|o2 = ideal(- y  + x )                        |
|                                             |
|o2 : Ideal of R                              |
+---------------------------------------------+
|i3 : eliminatingDual(origin R, I, {0}, 2)    |
|                                             |
|o3 = | y5+x2y2 y4+x2y y3+x2 xy2 xy x y2 y 1 ||
|                                             |
|o3 : DualSpace                               |
+---------------------------------------------+

This function generalizes truncatedDual, in that if v
includes all the variables in the ring, then its behavior is the same.

+-------------------------------------------+
|i4 : eliminatingDual(origin R, I, {0,1}, 2)|
|                                           |
|o4 = | xy y2 x y 1 |                       |
|                                           |
|o4 : DualSpace                             |
+-------------------------------------------+

See also truncatedDual.

Caveat
======

The space of dual elements satisfying the conditions is not in general of
finite dimension. If the dimension is infinite, this function will not
terminate.  This is not checked.  To ensure termination, the local dimension of
I at p should not exceed the length of v, and certain genericity constraints on
the coordinates must be met.

Ways to use eliminatingDual:
============================

  * "eliminatingDual(AbstractPoint,Ideal,List,ZZ)"
  * "eliminatingDual(AbstractPoint,Matrix,List,ZZ)"

For the programmer
==================

The object eliminatingDual, is a method function with options (Macaulay2Doc).

evaluate(DiffOp,AbstractPoint) -- evaluate coefficients of a differential operator
**********************************************************************************

Synopsis
========

  * Function: evaluate (NAGtypes),
  * Usage:
        evaluate(D, p)
  * Inputs:
      * D, a differential operator,
      * p, a point (NAGtypes), or a one-row matrix (Macaulay2Doc),
  * Outputs:
      * a differential operator,

Description
===========

Evaluates the polynomial coefficients of a differential operator,
at a point. This can be used to obtain a set of specialized Noetherian
operators.

+------------------------------+
|i1 : R = QQ[x,y];             |
+------------------------------+
|i2 : S = diffOpRing R;        |
+------------------------------+
|i3 : D = diffOp(x*dx+y*dy^3)  |
|                              |
|o3 = | ydy^3+xdx |            |
|                              |
|                1             |
|o3 : DiffOp in S              |
+------------------------------+
|i4 : evaluate(D, point{{1,2}})|
|                              |
|o4 = | 2dy^3+dx |             |
|                              |
|                1             |
|o4 : DiffOp in S              |
+------------------------------+

Ways to use this method:
========================

  * evaluate(DiffOp,AbstractPoint),
    -- evaluate coefficients of a differential operator
  * "evaluate(DiffOp,Matrix)"

gCorners -- generators of the initial ideal of a polynomial ideal
*****************************************************************

Synopsis
========

  * Usage:
        G = gCorners(p, I)
  * Inputs:
      * p, a point (NAGtypes),
      * I, an ideal (Macaulay2Doc), or a one-row matrix (Macaulay2Doc),
        of generators
  * Optional inputs (Macaulay2Doc),:
      * StandardBasis => ..., default value false
      * Tolerance, => ..., default
        value null, optional argument for numerical tolerance
  * Outputs:
      * G, a matrix (Macaulay2Doc), generators of the initial
        ideal in a one-row matrix

Description
===========

This method computes the generators of the initial ideal of an ideal, with
respect to a local order. These are precisely the monomials in the corners of
the staircase diagram of the initial ideal. The ring of the ideal should be
given a (global) monomial order and the local order will be taken to be the
reverse order. The point p is moved to the origin, so the monomial generators
represent terms of the Taylor expansion at p.

+----------------------------------------------------------------+
|i1 : R = CC[x,y];                                               |
+----------------------------------------------------------------+
|i2 : I = ideal{x^2-y^2}                                         |
|                                                                |
|            2    2                                              |
|o2 = ideal(x  - y )                                             |
|                                                                |
|o2 : Ideal of R                                                 |
+----------------------------------------------------------------+
|i3 : p = point matrix{{1,1}};                                   |
+----------------------------------------------------------------+
|i4 : gCorners(p, I)                                             |
|-- warning: experimental computation over inexact field begun   |
|--          results not reliable (one warning given per session)|
|                                                                |
|o4 = | y |                                                      |
|                                                                |
|             1      1                                           |
|o4 : Matrix R  <-- R                                            |
+----------------------------------------------------------------+

If the optional argument StandardBasis is set to true, the output is instead a
matrix of elements of the ideal with the point p translated to the origin such
that the lead terms generate the initial ideal, i.e., a standard basis. Note
that the coordinates of the standard basis elements are translated to be
centered at the point p.

+-------------------------------------------------+
|i5 : S = gCorners(p, I, StandardBasis=>true)     |
|                                                 |
|o5 = | -.5x2+.5y2-1x+y |                         |
|                                                 |
|             1      1                            |
|o5 : Matrix R  <-- R                             |
+-------------------------------------------------+
|i6 : R = CC[x,y,z];                              |
+-------------------------------------------------+
|i7 : J = ideal{z*(x*y-4), x-y}                   |
|                                                 |
|o7 = ideal (x*y*z - 4z, x - y)                   |
|                                                 |
|o7 : Ideal of R                                  |
+-------------------------------------------------+
|i8 : q = point matrix{{1.4142136, 1.4142136, 0}};|
+-------------------------------------------------+
|i9 : gCorners(q, J, Tolerance=>1e-5)             |
|                                                 |
|o9 = | y z |                                     |
|                                                 |
|             1      2                            |
|o9 : Matrix R  <-- R                             |
+-------------------------------------------------+
|i10 : gCorners(q, J, StandardBasis=>true)        |
|                                                 |
|o10 = | -1x+y -.5x2z-1.41421xz+z |               |
|                                                 |
|              1      2                           |
|o10 : Matrix R  <-- R                            |
+-------------------------------------------------+

Ways to use gCorners:
=====================

  * "gCorners(AbstractPoint,Ideal)"
  * "gCorners(AbstractPoint,Matrix)"

For the programmer
==================

The object gCorners, is a method function with options (Macaulay2Doc).

getIdealFromNoetherianOperators -- Computes a primary ideal corresponding to a list of Noetherian operators and a prime ideal
*****************************************************************************************************************************

Synopsis
========

  * Usage:
        Q = getIdealFromNoetherianOperators(L, P)
  * Inputs:
      * L, a list (Macaulay2Doc), of differential operators,
      * P, an ideal (Macaulay2Doc), a prime ideal
  * Outputs:
      * Q, an ideal (Macaulay2Doc), the primary ideal
        corresponding to L and P

Description
===========

This method contains an implementation of Algorithm 3.9 in the paper Primary
ideals and their differential equations (see https://arxiv.org/abs/2001.04700
). This method can be seen as the reverse operation of computing a set of
Noetherian operators for a primary ideal.

Let $R$ be a polynomial ring $R = K[x_1,\ldots,x_n]$ over a field $K$ of
characteristic zero. Consider the Weyl algebra $D = R<dx_1,\ldots,dx_n>$ and a
list of differential operators $L = \{L_1,\ldots,L_m\} \,\subset\, D$. Denote
by $\mathcal{E} \,\subset\, D$ the $R$-bisubmodule of $D$ that is generated by
$L_1,\ldots,L_m$. For a given prime ideal $P \,\subset\, R$, this method
computes the $P$-primary ideal given as $$ Q = \{ f \,\in\, R\, \mid\, \delta\,
\bullet\, f\, \in P, \ \forall  \delta \in \mathcal{E} \}. $$

Next, we provide several examples to show the interplay between computing a set
of Noetherian operators and then recovering the original ideal.

The first example shows an ideal that can be described with two different sets
of Noetherian operators (this example appeared in Example 7.8 of the paper
Primary ideals and their differential equations (see
https://arxiv.org/abs/2001.04700 )).

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x_1,x_2,x_3,x_4]                                                 |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : MM = matrix {{x_3,x_1,x_2},{x_1,x_2,x_4}}                               |
|                                                                             |
|o2 = | x_3 x_1 x_2 |                                                         |
|     | x_1 x_2 x_4 |                                                         |
|                                                                             |
|             2      3                                                        |
|o2 : Matrix R  <-- R                                                         |
+-----------------------------------------------------------------------------+
|i3 : P = minors(2,MM)                                                        |
|                                                                             |
|               2                           2                                 |
|o3 = ideal (- x  + x x , - x x  + x x , - x  + x x )                         |
|               1    2 3     1 2    3 4     2    1 4                          |
|                                                                             |
|o3 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i4 : M = ideal{x_1^2,x_2^2,x_3^2,x_4^2}                                      |
|                                                                             |
|             2   2   2   2                                                   |
|o4 = ideal (x , x , x , x )                                                  |
|             1   2   3   4                                                   |
|                                                                             |
|o4 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i5 : Q = joinIdeals(P,M);                                                    |
|                                                                             |
|o5 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i6 : L1 = noetherianOperators(Q) -- A set of Noetherian operators            |
|                                                                             |
|o6 = {| 1 |, | dx_1 |, | dx_2 |, | dx_1^2 |, | dx_1dx_2 |, | dx_2^2 |, |     |
|     ------------------------------------------------------------------------|
|     dx_1^3 |, | dx_1^2dx_2 |, | dx_1dx_2^2 |, | dx_2^3 |, |                 |
|     ------------------------------------------------------------------------|
|     2x_3x_4dx_1^3dx_2+5x_1x_4dx_1^2dx_2^2+2x_2x_4dx_1dx_2^3 |}              |
|                                                                             |
|o6 : List                                                                    |
+-----------------------------------------------------------------------------+
|i7 : Q1 = getIdealFromNoetherianOperators(L1, P);                            |
|                                                                             |
|o7 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i8 : Q == Q1                                                                 |
|                                                                             |
|o8 = true                                                                    |
+-----------------------------------------------------------------------------+
|i9 : L2 = noetherianOperators(M) -- Another set of Noetherian operators      |
|                                                                             |
|o9 = {| 1 |, | dx_1 |, | dx_2 |, | dx_3 |, | dx_4 |, | dx_1dx_2 |, | dx_1dx_3|
|     ------------------------------------------------------------------------|
|     |, | dx_1dx_4 |, | dx_2dx_3 |, | dx_2dx_4 |, | dx_3dx_4 |, |            |
|     ------------------------------------------------------------------------|
|     dx_1dx_2dx_3 |, | dx_1dx_2dx_4 |, | dx_1dx_3dx_4 |, | dx_2dx_3dx_4 |, | |
|     ------------------------------------------------------------------------|
|     dx_1dx_2dx_3dx_4 |}                                                     |
|                                                                             |
|o9 : List                                                                    |
+-----------------------------------------------------------------------------+
|i10 : Q2 = getIdealFromNoetherianOperators(L2, P);                           |
|                                                                             |
|o10 : Ideal of R                                                             |
+-----------------------------------------------------------------------------+
|i11 : Q == Q2                                                                |
|                                                                             |
|o11 = true                                                                   |
+-----------------------------------------------------------------------------+

The following example was given as the running example in the Introduction of
the paper Primary ideals and their differential equations (see
https://arxiv.org/abs/2001.04700 ).

+-----------------------------------------------------------------------------+
|i12 : Q = ideal(3*x_1^2*x_2^2-x_2^3*x_3-x_1^3*x_4-3*x_1*x_2*x_3*x_4+2*x_3^2*x|
|                                                                             |
|               2 2    3      3                   2 2    3           2        |
|o12 = ideal (3x x  - x x  - x x  - 3x x x x  + 2x x , 3x x x  - 3x x x x  -  |
|               1 2    2 3    1 4     1 2 3 4     3 4    1 2 4     1 2 3 4    |
|      -----------------------------------------------------------------------|
|        2   2       2 2     3       2    4         2         2   2    3      |
|      3x x x  + 3x x x  + 2x  - 2x x , 3x x  - 6x x x x  + 3x x x  + x  -    |
|        1 3 4     2 3 4     2     3 4    2 3     1 2 3 4     1 3 4    2      |
|      -----------------------------------------------------------------------|
|         2      3      4       2           2 2         2 2   5      3        |
|      x x , 4x x x  + x x  - 6x x x x  - 3x x x  + 4x x x , x  - x x x  -    |
|       3 4    1 2 3    1 4     1 2 3 4     2 3 4     1 3 4   2    1 2 4      |
|      -----------------------------------------------------------------------|
|       2   2        3     4    3              2    2 3   4      3 2     3    |
|      x x x  + x x x , x x  - x x x  - x x x x  + x x , x x  - x x  - 2x x x |
|       2 3 4    1 3 4   1 2    2 3 4    1 2 3 4    3 4   1 2    2 3     1 3 4|
|      -----------------------------------------------------------------------|
|              2     5     3           2 2     2 2         3      4           |
|      + 2x x x x , x  - 4x x x  + 3x x x  + 2x x x  - 2x x x , 3x x x  -     |
|          1 2 3 4   1     1 2 3     1 2 3     1 3 4     2 3 4    1 3 4       |
|      -----------------------------------------------------------------------|
|        2   2       2 3       3         2       2           2      3 3       |
|      6x x x x  + 3x x x  + 2x x  + 6x x x  - 6x x x  - 2x x x , 4x x  +     |
|        1 2 3 4     2 3 4     1 2     1 2 3     1 3 4     2 3 4    2 3       |
|      -----------------------------------------------------------------------|
|        3 2            3       4 2    4     2         2 2       2            |
|      4x x x  - 12x x x x  + 4x x  - x  + 6x x x  + 3x x  - 8x x x )         |
|        1 3 4      1 2 3 4     3 4    1     1 2 3     2 3     1 3 4          |
|                                                                             |
|o12 : Ideal of R                                                             |
|-----------------------------------------------------------------------------|
|_4^2,3*x_1^3*x_2*x_4-3*x_1*x_2^2*x_3*x_4-3*x_1^2*x_3*x_4^2+3*x_2*x_3^2*x_4^2+|
|-----------------------------------------------------------------------------|
|2*x_2^3-2*x_3*x_4^2,3*x_2^4*x_3-6*x_1*x_2^2*x_3*x_4+3*x_1^2*x_3*x_4^2+x_2^3-x|
|-----------------------------------------------------------------------------|
|_3*x_4^2,4*x_1*x_2^3*x_3+x_1^4*x_4-6*x_1^2*x_2*x_3*x_4-3*x_2^2*x_3^2*x_4+4*x_|
|-----------------------------------------------------------------------------|
|1*x_3^2*x_4^2,x_2^5-x_1*x_2^3*x_4-x_2^2*x_3*x_4^2+x_1*x_3*x_4^3,x_1*x_2^4-x_2|
|-----------------------------------------------------------------------------|
|^3*x_3*x_4-x_1*x_2*x_3*x_4^2+x_3^2*x_4^3,x_1^4*x_2-x_2^3*x_3^2-2*x_1^3*x_3*x_|
|-----------------------------------------------------------------------------|
|4+2*x_1*x_2*x_3^2*x_4,x_1^5-4*x_1^3*x_2*x_3+3*x_1*x_2^2*x_3^2+2*x_1^2*x_3^2*x|
|-----------------------------------------------------------------------------|
|_4-2*x_2*x_3^3*x_4,3*x_1^4*x_3*x_4-6*x_1^2*x_2*x_3^2*x_4+3*x_2^2*x_3^3*x_4+2*|
|-----------------------------------------------------------------------------|
|x_1^3*x_2+6*x_1*x_2^2*x_3-6*x_1^2*x_3*x_4-2*x_2*x_3^2*x_4,4*x_2^3*x_3^3+4*x_1|
|-----------------------------------------------------------------------------|
|^3*x_3^2*x_4-12*x_1*x_2*x_3^3*x_4+4*x_3^4*x_4^2-x_1^4+6*x_1^2*x_2*x_3+3*x_2^2|
|-----------------------------------------------------------------------------|
|*x_3^2-8*x_1*x_3^2*x_4)                                                      |
+-----------------------------------------------------------------------------+
|i13 : L = noetherianOperators(Q)                                             |
|                                                                             |
|o13 = {| 1 |, | dx_1 |, | dx_1^2-2x_2dx_2 |}                                 |
|                                                                             |
|o13 : List                                                                   |
+-----------------------------------------------------------------------------+
|i14 : Q' = getIdealFromNoetherianOperators(L, P);                            |
|                                                                             |
|o14 : Ideal of R                                                             |
+-----------------------------------------------------------------------------+
|i15 : Q == Q'                                                                |
|                                                                             |
|o15 = true                                                                   |
+-----------------------------------------------------------------------------+

The next example was given by Palamodov to show that there exist primary ideals
that cannot be described by using differential operators with constant
coefficients.

+-------------------------------------------------------+
|i16 : R = QQ[x_1, x_2, x_3]                            |
|                                                       |
|o16 = R                                                |
|                                                       |
|o16 : PolynomialRing                                   |
+-------------------------------------------------------+
|i17 : Q = ideal(x_1^2, x_2^2, x_1-x_2*x_3)             |
|                                                       |
|              2   2                                    |
|o17 = ideal (x , x , - x x  + x )                      |
|              1   2     2 3    1                       |
|                                                       |
|o17 : Ideal of R                                       |
+-------------------------------------------------------+
|i18 : L = noetherianOperators(Q)                       |
|                                                       |
|o18 = {| 1 |, | x_3dx_1+dx_2 |}                        |
|                                                       |
|o18 : List                                             |
+-------------------------------------------------------+
|i19 : Q' = getIdealFromNoetherianOperators(L,radical Q)|
|                                                       |
|                         2         2                   |
|o19 = ideal (x x  - x , x , x x , x )                  |
|              2 3    1   2   1 2   1                   |
|                                                       |
|o19 : Ideal of R                                       |
+-------------------------------------------------------+
|i20 : Q == Q'                                          |
|                                                       |
|o20 = true                                             |
+-------------------------------------------------------+

For the last example we consider an ideal defined by using the join
construction.

+-----------------------------------------------------------------------------+
|i21 : R = QQ[x_1..x_9]                                                       |
|                                                                             |
|o21 = R                                                                      |
|                                                                             |
|o21 : PolynomialRing                                                         |
+-----------------------------------------------------------------------------+
|i22 : MM = genericMatrix(R, 3, 3)                                            |
|                                                                             |
|o22 = | x_1 x_4 x_7 |                                                        |
|      | x_2 x_5 x_8 |                                                        |
|      | x_3 x_6 x_9 |                                                        |
|                                                                             |
|              3      3                                                       |
|o22 : Matrix R  <-- R                                                        |
+-----------------------------------------------------------------------------+
|i23 : P = minors(2, MM)                                                      |
|                                                                             |
|o23 = ideal (- x x  + x x , - x x  + x x , - x x  + x x , - x x  + x x , -   |
|                2 4    1 5     3 4    1 6     3 5    2 6     2 7    1 8      |
|      -----------------------------------------------------------------------|
|      x x  + x x , - x x  + x x , - x x  + x x , - x x  + x x , - x x  +     |
|       3 7    1 9     3 8    2 9     5 7    4 8     6 7    4 9     6 8       |
|      -----------------------------------------------------------------------|
|      x x )                                                                  |
|       5 9                                                                   |
|                                                                             |
|o23 : Ideal of R                                                             |
+-----------------------------------------------------------------------------+
|i24 : M = ideal(x_1^2, x_5^2, x_9^2, x_2, x_3, x_4, x_6, x_7, x_8)           |
|                                                                             |
|              2   2   2                                                      |
|o24 = ideal (x , x , x , x , x , x , x , x , x )                             |
|              1   5   9   2   3   4   6   7   8                              |
|                                                                             |
|o24 : Ideal of R                                                             |
+-----------------------------------------------------------------------------+
|i25 : Q = joinIdeals(P, M)                                                   |
|                                                                             |
|                               2 2                2 2                        |
|o25 = ideal (x x x  - x x x , x x  - 2x x x x  + x x , x x x x  - 2x x x x  +|
|              2 6 7    3 4 8   3 8     2 3 8 9    2 9   3 6 7 8     3 4 8 9  |
|      -----------------------------------------------------------------------|
|           2   2 2                2 2   2 2                2 2   2 2         |
|      x x x , x x  - 2x x x x  + x x , x x  - 2x x x x  + x x , x x  -       |
|       2 4 9   6 7     4 6 7 9    4 9   5 7     4 5 7 8    4 8   2 7         |
|      -----------------------------------------------------------------------|
|                   2 2     2                                                 |
|      2x x x x  + x x , x x x  - 2x x x x  + x x x x , x x x x  - 2x x x x  +|
|        1 2 7 8    1 8   3 5 7     3 4 5 8    2 4 6 8   2 3 4 7     1 3 4 8  |
|      -----------------------------------------------------------------------|
|       2       2 2                2 2   2 2                2 2   2   2       |
|      x x x , x x  - 2x x x x  + x x , x x  - 2x x x x  + x x , x x x  -     |
|       1 6 8   3 5     2 3 5 6    2 6   3 4     1 3 4 6    1 6   6 7 8       |
|      -----------------------------------------------------------------------|
|            2      2   2           2     2 2       2          2 2            |
|      2x x x x  - x x x  + 2x x x x , x x x  - 2x x x x  - x x x  +          |
|        4 6 8 9    5 7 9     4 5 8 9   3 6 8     2 6 8 9    3 5 9            |
|      -----------------------------------------------------------------------|
|              2          2      2 2       2          2 2   2 2           2   |
|      2x x x x , 2x x x x  - x x x  - 2x x x x  + x x x , x x x  - 2x x x x  |
|        2 5 6 9    3 5 6 8    2 6 8     3 5 8 9    2 5 9   3 7 8     2 3 7 9 |
|      -----------------------------------------------------------------------|
|                2    2   2      2          2 2     2            2 2   2   2  |
|      + 2x x x x  - x x x , 2x x x x  - x x x  - 2x x x x  + x x x , x x x  -|
|          1 2 7 9    1 8 9    5 6 7 8    4 6 8     5 6 7 9    4 5 9   3 6 7  |
|      -----------------------------------------------------------------------|
|        2                 2    2   2   2   2           2     2               |
|      2x x x x  + 2x x x x  - x x x , x x x  - 2x x x x  + 2x x x x  -       |
|        3 4 7 9     1 3 4 9    1 6 9   3 4 7     1 3 6 7     1 6 7 9         |
|      -----------------------------------------------------------------------|
|       2   2     2 2       2         2          2   2    2               2   |
|      x x x , x x x  - 2x x x x  + 2x x x x  - x x x , 2x x x x  - 2x x x x  |
|       1 4 9   2 3 7     1 3 7 8     1 3 8 9    1 2 9    2 4 5 7     1 2 5 7 |
|      -----------------------------------------------------------------------|
|         2 2      2 2     2 2      2 2           2       2              2    |
|      - x x x  + x x x , x x x  - x x x  - 2x x x x  + 2x x x x , 2x x x x  -|
|         2 4 8    1 5 8   2 4 7    1 5 7     1 2 4 8     1 4 5 8    2 3 4 5  |
|      -----------------------------------------------------------------------|
|              2    2 2      2 2     2   2    2   2       2         2         |
|      2x x x x  - x x x  + x x x , x x x  - x x x  - 2x x x x  + 2x x x x ,  |
|        1 3 4 5    2 4 6    1 5 6   2 3 4    1 3 5     1 2 4 6     1 2 5 6   |
|      -----------------------------------------------------------------------|
|       3 3       2 2       2     2    3 3     2 3     2                   2  |
|      x x  - 3x x x x  + 3x x x x  - x x , x x x  + 2x x x x x  - 4x x x x x |
|       6 8     5 6 8 9     5 6 8 9    5 9   4 6 8     5 6 7 8 9     4 5 6 8 9|
|      -----------------------------------------------------------------------|
|          3   2       2   2         3           2        2   2             2 |
|      - 2x x x  + 3x x x x , x x x x  - 2x x x x x  - x x x x  + 2x x x x x ,|
|          5 7 9     4 5 8 9   3 4 6 8     2 4 6 8 9    2 5 7 9     2 4 5 8 9 |
|      -----------------------------------------------------------------------|
|         3 2       2               2           3 2       2   2            2  |
|      x x x  + 2x x x x x  - 4x x x x x  - 2x x x  + 3x x x x , 2x x x x x  -|
|       2 6 8     3 5 6 8 9     2 5 6 8 9     3 5 9     2 5 6 9    1 3 4 6 8  |
|      -----------------------------------------------------------------------|
|       2 2 2         2        2 2 2    2            2 2 2         2          |
|      x x x  - 2x x x x x  + x x x , 2x x x x x  + x x x  + 2x x x x x  -    |
|       1 6 8     2 3 4 8 9    2 4 9    3 4 5 7 8    1 6 8     2 3 4 8 9      |
|      -----------------------------------------------------------------------|
|                                        2                     2    2 2 2     |
|      8x x x x x x  - 4x x x x x x  + 4x x x x x  + 4x x x x x  - x x x ,    |
|        1 3 4 5 8 9     1 2 4 6 8 9     1 5 6 8 9     1 2 4 5 9    1 5 9     |
|      -----------------------------------------------------------------------|
|       3             2           2                     2    2     2   3 3    |
|      x x x x  - 4x x x x x  + 2x x x x x  + 2x x x x x  - x x x x , x x  -  |
|       3 4 7 8     1 3 4 8 9     1 3 6 8 9     1 2 3 4 9    1 2 6 9   3 7    |
|      -----------------------------------------------------------------------|
|          2 2       2     2    3 3   2 2            3           2            |
|      3x x x x  + 3x x x x  - x x , x x x x  - x x x x  + 2x x x x x  -      |
|        1 3 7 9     1 3 7 9    1 9   1 5 6 7    2 3 4 8     1 2 4 6 8        |
|      -----------------------------------------------------------------------|
|        2                  2       2 3      2 3           2         2   2    |
|      2x x x x x , 2x x x x x  - 2x x x  + x x x  - 4x x x x x  + 3x x x x , |
|        1 4 5 6 8    1 2 4 5 7     1 5 7    2 4 8     1 2 4 5 8     1 4 5 8  |
|      -----------------------------------------------------------------------|
|                2     2   3    3 2         2           2   2     3 3         |
|      2x x x x x  - 2x x x  + x x x  - 4x x x x x  + 3x x x x , x x  -       |
|        1 2 3 4 5     1 3 5    2 4 6     1 2 4 5 6     1 2 5 6   2 4         |
|      -----------------------------------------------------------------------|
|          2 2       2     2    3 3   2 2 3         2 2               2       |
|      3x x x x  + 3x x x x  - x x , x x x  + 2x x x x x  - 4x x x x x x  -   |
|        1 2 4 5     1 2 4 5    1 5   1 6 8     2 3 4 8 9     1 2 4 6 8 9     |
|      -----------------------------------------------------------------------|
|            2   2    2 2   2               2   2     3     2     2           |
|      2x x x x x  - x x x x  + 4x x x x x x , x x x x  - 2x x x x x  -       |
|        1 2 5 7 9    2 4 8 9     1 2 4 5 8 9   1 3 6 8     1 2 6 8 9         |
|      -----------------------------------------------------------------------|
|       3     2       2     2   2 3 2         2                 2             |
|      x x x x  + 2x x x x x , x x x  + 2x x x x x x  - 4x x x x x x  -       |
|       2 4 7 9     1 2 4 8 9   1 6 8     2 3 4 6 8 9     1 2 4 6 8 9         |
|      -----------------------------------------------------------------------|
|              2 2    2 2   2               2         2 2    2   2 2          |
|      2x x x x x  - x x x x  + 4x x x x x x , x x x x x  - x x x x  -        |
|        1 3 4 5 9    2 4 6 9     1 2 4 5 6 9   1 2 4 6 8    1 5 6 8          |
|      -----------------------------------------------------------------------|
|       2 2          2 2          2 2   2          2 2        2   2    2   2 2|
|      x x x x x  + x x x x x  + x x x x  - x x x x x , 2x x x x x  - x x x x |
|       2 4 6 8 9    1 5 6 8 9    2 4 5 9    1 2 4 5 9    2 3 4 6 8    1 5 6 8|
|      -----------------------------------------------------------------------|
|                2         2 2                            2 2                 |
|      - 2x x x x x x  - 5x x x x x  + 4x x x x x x x  + x x x x x  +         |
|          1 3 4 5 8 9     2 4 6 8 9     1 2 4 5 6 8 9    1 5 6 8 9           |
|      -----------------------------------------------------------------------|
|        2 2   2           2 2    2 3 2   2 3         2   2            3 2    |
|      4x x x x  - 4x x x x x  + x x x , x x x x  - 2x x x x x  - x x x x  +  |
|        2 4 5 9     1 2 4 5 9    1 5 9   1 6 7 8     1 4 6 8 9    2 3 4 9    |
|      -----------------------------------------------------------------------|
|            2   2   2 2 2 2    2 2 2 2     2 2                   2           |
|      2x x x x x , x x x x  - x x x x  - 4x x x x x x  + 4x x x x x x x  +   |
|        1 2 4 6 9   2 4 6 8    1 5 6 8     2 4 5 6 8 9     1 2 4 5 6 8 9     |
|      -----------------------------------------------------------------------|
|        2 2 2 2           3 2    2 4 2                                       |
|      3x x x x  - 4x x x x x  + x x x )                                      |
|        2 4 5 9     1 2 4 5 9    1 5 9                                       |
|                                                                             |
|o25 : Ideal of R                                                             |
+-----------------------------------------------------------------------------+
|i26 : L = noetherianOperators(Q)                                             |
|                                                                             |
|o26 = {| 1 |, | dx_5 |, | x_6dx_4+x_9dx_7 |, | x_4dx_4-x_8dx_8 |, |          |
|      -----------------------------------------------------------------------|
|      x_6dx_4dx_5+x_9dx_5dx_7 |, | x_4dx_4dx_5-x_8dx_5dx_8 |, |              |
|      -----------------------------------------------------------------------|
|      x_4x_5x_9dx_4dx_7+x_4x_8x_9dx_7^2+x_5^2x_9dx_4dx_8+x_5x_8x_9dx_7dx_8-x_|
|      -----------------------------------------------------------------------|
|      5x_6dx_4 |, | x_4x_5x_9dx_4dx_5dx_7+x_4x_8x_9dx_5dx_7^2+x_5^2x_9dx_4dx_|
|      -----------------------------------------------------------------------|
|      5dx_8+x_5x_8x_9dx_5dx_7dx_8-x_5x_6dx_4dx_5 |}                          |
|                                                                             |
|o26 : List                                                                   |
+-----------------------------------------------------------------------------+
|i27 : Q' = getIdealFromNoetherianOperators(L, radical Q)                     |
|                                                                             |
|                               2 2                2 2                        |
|o27 = ideal (x x x  - x x x , x x  - 2x x x x  + x x , x x x x  - 2x x x x  +|
|              2 6 7    3 4 8   3 8     2 3 8 9    2 9   3 6 7 8     3 4 8 9  |
|      -----------------------------------------------------------------------|
|           2   2 2                2 2   2 2                2 2   2 2         |
|      x x x , x x  - 2x x x x  + x x , x x  - 2x x x x  + x x , x x  -       |
|       2 4 9   6 7     4 6 7 9    4 9   5 7     4 5 7 8    4 8   2 7         |
|      -----------------------------------------------------------------------|
|                   2 2     2                                                 |
|      2x x x x  + x x , x x x  - 2x x x x  + x x x x , x x x x  - 2x x x x  +|
|        1 2 7 8    1 8   3 5 7     3 4 5 8    2 4 6 8   2 3 4 7     1 3 4 8  |
|      -----------------------------------------------------------------------|
|       2       2 2                2 2   2 2                2 2   2   2       |
|      x x x , x x  - 2x x x x  + x x , x x  - 2x x x x  + x x , x x x  -     |
|       1 6 8   3 5     2 3 5 6    2 6   3 4     1 3 4 6    1 6   6 7 8       |
|      -----------------------------------------------------------------------|
|            2      2   2           2     2 2       2          2 2            |
|      2x x x x  - x x x  + 2x x x x , x x x  - 2x x x x  - x x x  +          |
|        4 6 8 9    5 7 9     4 5 8 9   3 6 8     2 6 8 9    3 5 9            |
|      -----------------------------------------------------------------------|
|              2          2      2 2       2          2 2   2 2           2   |
|      2x x x x , 2x x x x  - x x x  - 2x x x x  + x x x , x x x  - 2x x x x  |
|        2 5 6 9    3 5 6 8    2 6 8     3 5 8 9    2 5 9   3 7 8     2 3 7 9 |
|      -----------------------------------------------------------------------|
|                2    2   2      2          2 2     2            2 2   2   2  |
|      + 2x x x x  - x x x , 2x x x x  - x x x  - 2x x x x  + x x x , x x x  -|
|          1 2 7 9    1 8 9    5 6 7 8    4 6 8     5 6 7 9    4 5 9   3 6 7  |
|      -----------------------------------------------------------------------|
|        2                 2    2   2   2   2           2     2               |
|      2x x x x  + 2x x x x  - x x x , x x x  - 2x x x x  + 2x x x x  -       |
|        3 4 7 9     1 3 4 9    1 6 9   3 4 7     1 3 6 7     1 6 7 9         |
|      -----------------------------------------------------------------------|
|       2   2     2 2       2         2          2   2    2               2   |
|      x x x , x x x  - 2x x x x  + 2x x x x  - x x x , 2x x x x  - 2x x x x  |
|       1 4 9   2 3 7     1 3 7 8     1 3 8 9    1 2 9    2 4 5 7     1 2 5 7 |
|      -----------------------------------------------------------------------|
|         2 2      2 2     2 2      2 2           2       2              2    |
|      - x x x  + x x x , x x x  - x x x  - 2x x x x  + 2x x x x , 2x x x x  -|
|         2 4 8    1 5 8   2 4 7    1 5 7     1 2 4 8     1 4 5 8    2 3 4 5  |
|      -----------------------------------------------------------------------|
|              2    2 2      2 2     2   2    2   2       2         2         |
|      2x x x x  - x x x  + x x x , x x x  - x x x  - 2x x x x  + 2x x x x ,  |
|        1 3 4 5    2 4 6    1 5 6   2 3 4    1 3 5     1 2 4 6     1 2 5 6   |
|      -----------------------------------------------------------------------|
|       3 3       2 2       2     2    3 3    2            2 2 2         2    |
|      x x  - 3x x x x  + 3x x x x  - x x , 2x x x x x  + x x x  + 2x x x x x |
|       6 8     5 6 8 9     5 6 8 9    5 9    3 4 5 7 8    1 6 8     2 3 4 8 9|
|      -----------------------------------------------------------------------|
|                                          2                     2    2 2 2   |
|      - 8x x x x x x  - 4x x x x x x  + 4x x x x x  + 4x x x x x  - x x x ,  |
|          1 3 4 5 8 9     1 2 4 6 8 9     1 5 6 8 9     1 2 4 5 9    1 5 9   |
|      -----------------------------------------------------------------------|
|       3 3       2 2       2     2    3 3   3 3       2 2       2     2      |
|      x x  - 3x x x x  + 3x x x x  - x x , x x  - 3x x x x  + 3x x x x  -    |
|       3 7     1 3 7 9     1 3 7 9    1 9   2 4     1 2 4 5     1 2 4 5      |
|      -----------------------------------------------------------------------|
|       3 3                                                                   |
|      x x )                                                                  |
|       1 5                                                                   |
|                                                                             |
|o27 : Ideal of R                                                             |
+-----------------------------------------------------------------------------+
|i28 : Q == Q'                                                                |
|                                                                             |
|o28 = true                                                                   |
+-----------------------------------------------------------------------------+

Ways to use getIdealFromNoetherianOperators:
============================================

  * "getIdealFromNoetherianOperators(List,Ideal)"

For the programmer
==================

The object getIdealFromNoetherianOperators,
is a method function (Macaulay2Doc).

getModuleFromNoetherianOperators -- Computes a primary submodule corresponding to a list of Noetherian operators and a prime ideal
**********************************************************************************************************************************

Synopsis
========

  * Usage:
        U = getModuleFromNoetherianOperators(P, L)
  * Inputs:
      * P, an ideal (Macaulay2Doc), a prime ideal
      * L, a list (Macaulay2Doc), of differential operators,
  * Outputs:
      * U, a module (Macaulay2Doc), the primary submodule
        corresponding to L and P

Description
===========

This method contains an implementation of Algorithm 4.3 in the paper Primary
decomposition of modules: a computational differential approach (see
https://arxiv.org/abs/2104.03385 ). This method can be seen as the reverse
operation of computing a set of Noetherian operators for a primary module. For
more details, see Section 4 of the paper Primary decomposition of modules: a
computational differential approach (see https://arxiv.org/abs/2104.03385 ).

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x1,x2,x3,x4]                                                     |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : U = image matrix{{x1*x2,x2*x3,x3*x4,x4*x1}, {x1^2,x2^2,x3^2,x4^2}}      |
|                                                                             |
|o2 = image | x1x2 x2x3 x3x4 x1x4 |                                           |
|           | x1^2 x2^2 x3^2 x4^2 |                                           |
|                                                                             |
|                             2                                               |
|o2 : R-module, submodule of R                                                |
+-----------------------------------------------------------------------------+
|i3 : dpd = differentialPrimaryDecomposition U                                |
|                                                                             |
|                                                                             |
|o3 = {{ideal (x4, x2), {| 1 |}}, {ideal (x3, x1), {| 1 |}}, {ideal (x4, x3,  |
|                        | 0 |                      | 0 |                     |
|     ------------------------------------------------------------------------|
|                                                                             |
|     x2), {| x1dx2 |, | x1dx2^2 |, | x1^2dx2^3+6x1dx2dx3 |, |                |
|           |   -1  |  |  -2dx2  |  |    -3x1dx2^2-6dx3   |  |                |
|     ------------------------------------------------------------------------|
|                                                                             |
|     x1^2dx2^4+12x1dx2^2dx3 |}}, {ideal (x4, x3, x1), {| x2dx3 |, | x2dx3^2  |
|       -4x1dx2^3-24dx2dx3   |                          |   -1  |  |  -2dx3   |
|     ------------------------------------------------------------------------|
|                                                                             |
|     |, | x2^2dx3^3+6x2dx3dx4 |, | x2^2dx3^4+12x2dx3^2dx4 |}}, {ideal (x4,   |
|     |  |    -3x2dx3^2-6dx4   |  |   -4x2dx3^3-24dx3dx4   |                  |
|     ------------------------------------------------------------------------|
|                                                                             |
|     x2, x1), {| x3dx4 |, | x3dx4^2 |, | x3^2dx4^3+6x3dx1dx4 |, |            |
|               |   -1  |  |  -2dx4  |  |    -3x3dx4^2-6dx1   |  |            |
|     ------------------------------------------------------------------------|
|                                                                             |
|     x3^2dx4^4+12x3dx1dx4^2 |}}, {ideal (x3, x2, x1), {| x4dx1 |, | x4dx1^2  |
|       -4x3dx4^3-24dx1dx4   |                          |   -1  |  |  -2dx1   |
|     ------------------------------------------------------------------------|
|                                                                             |
|     |, | x4^2dx1^3+6x4dx1dx2 |, | x4^2dx1^4+12x4dx1^2dx2 |}}, {ideal (x3 -  |
|     |  |    -3x4dx1^2-6dx2   |  |   -4x4dx1^3-24dx1dx2   |                  |
|     ------------------------------------------------------------------------|
|                                                                             |
|     x4, x2 - x4, x1 - x4), {| -1 |}}, {ideal (x3 + x4, x2 - x4, x1 + x4), {||
|                             |  1 |                                         ||
|     ------------------------------------------------------------------------|
|                                        2     2                              |
|     1 |}}, {ideal (x2 + x4, x1 + x3, x3  + x4 ), {| x4 |}}, {ideal (x4, x3, |
|     1 |                                           | x3 |                    |
|     ------------------------------------------------------------------------|
|     x2, x1), {|      -2dx1dx2dx3^2-2dx1^2dx3dx4      |}}}                   |
|               | 2dx1dx2^2dx3+dx1^2dx3^2+2dx1dx3dx4^2 |                      |
|                                                                             |
|o3 : List                                                                    |
+-----------------------------------------------------------------------------+
|i4 : M = dpd / (L -> getModuleFromNoetherianOperators(first L, last L)) //   |
|                                                                             |
|o4 = image | x1x4 x3x4 x2x3 x1x2 |                                           |
|           | x4^2 x3^2 x2^2 x1^2 |                                           |
|                                                                             |
|                             2                                               |
|o4 : R-module, submodule of R                                                |
|-----------------------------------------------------------------------------|
|intersect                                                                    |
+-----------------------------------------------------------------------------+
|i5 : M == U                                                                  |
|                                                                             |
|o5 = true                                                                    |
+-----------------------------------------------------------------------------+

Ways to use getModuleFromNoetherianOperators:
=============================================

  * "getModuleFromNoetherianOperators(Ideal,List)"

For the programmer
==================

The object getModuleFromNoetherianOperators,
is a method function (Macaulay2Doc).

hilbertFunction(ZZ,DualSpace)
*****************************

Synopsis
========

  * Function: hilbertFunction (Macaulay2Doc),
  * Usage:
        k = hilbertFunction(d,S)
        T = hilbertFunction(S)
        K = hilbertFunction(D,S)
  * Inputs:
      * d, an integer, or a list (Macaulay2Doc), of arguments
      * S, an instance of the type DualSpace (NAGtypes), a
        dual space
  * Outputs:
      * T, a tally (Macaulay2Doc), of dimensions for each degree
      * k, an integer (Macaulay2Doc), dimensions
      * K, a list (Macaulay2Doc), of dimensions

Description
===========

This function counts the dimension of the dual space in each degree.  A single
degree can be specified, or a list of degrees.  If no degree is specified, then
a Tally is returned pairing each degree with its dimension.

+-----------------------------------------------------------+
|i1 : R = CC[x,y];                                          |
+-----------------------------------------------------------+
|i2 : I = (ideal {x,y})^4;                                  |
|                                                           |
|o2 : Ideal of R                                            |
+-----------------------------------------------------------+
|i3 : D = zeroDimensionalDual(origin R, I)                  |
|                                                           |
|o3 = | 1 y x .5y2 xy .5x2 .166667y3 .5xy2 .5x2y .166667x3 ||
|                                                           |
|o3 : DualSpace                                             |
+-----------------------------------------------------------+
|i4 : hilbertFunction D                                     |
|                                                           |
|o4 = Tally{0 => 1}                                         |
|           1 => 2                                          |
|           2 => 3                                          |
|           3 => 4                                          |
|                                                           |
|o4 : Tally                                                 |
+-----------------------------------------------------------+

Ways to use this method:
========================

  * "hilbertFunction(DualSpace)"
  * "hilbertFunction(List,DualSpace)"
  * hilbertFunction(ZZ,DualSpace)

isPointEmbedded(AbstractPoint,Ideal,List) -- numerically determine if the point is an embedded component of the scheme
**********************************************************************************************************************

Synopsis
========

  * Function: isPointEmbedded,
  * Usage:
        B = isPointEmbedded(P,I,C)
  * Inputs:
      * P
      * I
      * C,  witness sets representing components of Spec(I) containing P
  * Optional inputs (Macaulay2Doc),:
      * AllVisible => ..., default value false
  * Outputs:
      * B, a Boolean value (Macaulay2Doc),

Description
===========

This method runs a numerical embedded component test described in the following
paper.

R. Krone and A. Leykin, "Numerical algorithms for detecting embedded
components.", arXiv:1405.7871

See also
========

  * isPointEmbeddedInCurve
    -- numerically determine
    if the point is an embedded component of a 1-dimensional scheme

Ways to use this method:
========================

  * isPointEmbedded(AbstractPoint,Ideal,List)
    -- numerically determine
    if the point is an embedded component of the scheme

isPointEmbeddedInCurve(AbstractPoint,Ideal) -- numerically determine if the point is an embedded component of a 1-dimensional scheme
************************************************************************************************************************************

Synopsis
========

  * Function: isPointEmbeddedInCurve,
  * Usage:
        B = isPointEmbeddedInCurve(P,I)
  * Inputs:
      * P
      * I
  * Optional inputs (Macaulay2Doc),:
      * Regularity => ..., default value -1
  * Outputs:
      * B, a Boolean value (Macaulay2Doc),

Description
===========

This method runs a numerical embedded component test described in the following
paper.

R. Krone and A. Leykin, "Numerical algorithms for detecting embedded
components.", arXiv:1405.7871

The optional parameter Regularity can be used to set the local Hilbert
regularity. The default value -1 computes it using
localHilbertRegularity,

See also
========

  * isPointEmbeddedInCurve
    -- numerically determine
    if the point is an embedded component of a 1-dimensional scheme

Ways to use this method:
========================

  * isPointEmbeddedInCurve(AbstractPoint,Ideal)
    -- numerically determine
    if the point is an embedded component of a 1-dimensional scheme

joinIdeals -- Computes the join of two ideals
*********************************************

Synopsis
========

  * Usage:
        K = joinIdeals(I, J)
  * Inputs:
      * I, an ideal (Macaulay2Doc), an ideal
      * J, an ideal (Macaulay2Doc), an ideal
  * Outputs:
      * K, an ideal (Macaulay2Doc), the join of the ideals I and J

Description
===========

This method computes the join of two ideals. The join can be used to describe
very interesting types of primary ideals that include the symbolic powers of
prime ideals. For more details reader is referred to Section 7 of the paper
Primary ideals and their differential equations (see
https://arxiv.org/abs/2001.04700 ).

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x_1..x_9]                                                        |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : MM = genericMatrix(R, 3, 3)                                             |
|                                                                             |
|o2 = | x_1 x_4 x_7 |                                                         |
|     | x_2 x_5 x_8 |                                                         |
|     | x_3 x_6 x_9 |                                                         |
|                                                                             |
|             3      3                                                        |
|o2 : Matrix R  <-- R                                                         |
+-----------------------------------------------------------------------------+
|i3 : P = minors(2, MM)                                                       |
|                                                                             |
|o3 = ideal (- x x  + x x , - x x  + x x , - x x  + x x , - x x  + x x , -    |
|               2 4    1 5     3 4    1 6     3 5    2 6     2 7    1 8       |
|     ------------------------------------------------------------------------|
|     x x  + x x , - x x  + x x , - x x  + x x , - x x  + x x , - x x  + x x )|
|      3 7    1 9     3 8    2 9     5 7    4 8     6 7    4 9     6 8    5 9 |
|                                                                             |
|o3 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i4 : M = ideal(x_1^2, x_2^2, x_3^2, x_4, x_5, x_6, x_7, x_8, x_9)            |
|                                                                             |
|             2   2   2                                                       |
|o4 = ideal (x , x , x , x , x , x , x , x , x )                              |
|             1   2   3   4   5   6   7   8   9                               |
|                                                                             |
|o4 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i5 : Q = joinIdeals(P, M)                                                    |
|                                                                             |
|                                                    3 3       2 2            |
|o5 = ideal (x x  - x x , x x  - x x , x x  - x x , x x  - 3x x x x  +        |
|             6 8    5 9   6 7    4 9   5 7    4 8   3 8     2 3 8 9          |
|     ------------------------------------------------------------------------|
|       2     2    3 3   3   2       2           2     2    3   2     2 2     |
|     3x x x x  - x x , x x x  - 3x x x x x  + 3x x x x  - x x x , x x x x  - |
|       2 3 8 9    2 9   3 5 8     2 3 5 8 9     2 3 5 9    2 6 9   2 3 7 8   |
|     ------------------------------------------------------------------------|
|        2   2    2   2      2   2        2   2    2     2     2              |
|     x x x x  - x x x x  + x x x x  + x x x x  - x x x x , x x x x x  -      |
|      1 3 7 8    2 3 7 9    1 3 8 9    1 2 7 9    1 2 8 9   2 3 4 7 8        |
|     ------------------------------------------------------------------------|
|        2   2    2            2              2   2    2     2   3 2          |
|     x x x x  - x x x x x  + x x x x x  + x x x x  - x x x x , x x x  -      |
|      1 3 4 8    2 3 4 7 9    1 3 5 8 9    1 2 4 9    1 2 5 9   3 5 8        |
|     ------------------------------------------------------------------------|
|         2 2       2            3 2       2 2        2          2   2        |
|     3x x x x  + 3x x x x x  - x x x , x x x x  - x x x x x  - x x x x  +    |
|       2 3 5 9     2 3 5 6 9    2 6 9   2 3 4 8    1 3 4 5 8    2 3 4 9      |
|     ------------------------------------------------------------------------|
|      2   2        2          2           3 3       2 2       2     2    3 3 |
|     x x x x  + x x x x x  - x x x x x , x x  - 3x x x x  + 3x x x x  - x x ,|
|      1 3 5 9    1 2 4 6 9    1 2 5 6 9   3 7     1 3 7 9     1 3 7 9    1 9 |
|     ------------------------------------------------------------------------|
|      3 3       2 2       2     2    3 3   3   2       2           2     2   |
|     x x  - 3x x x x  + 3x x x x  - x x , x x x  - 3x x x x x  + 3x x x x  - |
|      2 7     1 2 7 8     1 2 7 8    1 8   3 4 7     1 3 4 7 9     1 3 4 9   |
|     ------------------------------------------------------------------------|
|      3   2   3   2       2           2     2    3   2   3 2         2 2     |
|     x x x , x x x  - 3x x x x x  + 3x x x x  - x x x , x x x  - 3x x x x  + |
|      1 6 9   2 4 7     1 2 4 7 8     1 2 4 8    1 5 8   3 4 7     1 3 4 9   |
|     ------------------------------------------------------------------------|
|       2            3 2     3 2         2 2       2            3 2     3 3   |
|     3x x x x x  - x x x , x x x  - 3x x x x  + 3x x x x x  - x x x , x x  - |
|       1 3 4 6 9    1 6 9   2 4 7     1 2 4 8     1 2 4 5 8    1 5 8   3 5   |
|     ------------------------------------------------------------------------|
|         2 2       2     2    3 3     2 2        2   2    2   2      2   2   |
|     3x x x x  + 3x x x x  - x x , x x x x  - x x x x  - x x x x  + x x x x  |
|       2 3 5 6     2 3 5 6    2 6   2 3 4 5    1 3 4 5    2 3 4 6    1 3 5 6 |
|     ------------------------------------------------------------------------|
|          2   2    2     2   3 3       2 2       2     2    3 3   3 3        |
|     + x x x x  - x x x x , x x  - 3x x x x  + 3x x x x  - x x , x x  -      |
|        1 2 4 6    1 2 5 6   3 4     1 3 4 6     1 3 4 6    1 6   2 4        |
|     ------------------------------------------------------------------------|
|         2 2       2     2    3 3       2   3    2 2 4    3   2              |
|     3x x x x  + 3x x x x  - x x , x x x x x  - x x x  + x x x x x  -        |
|       1 2 4 5     1 2 4 5    1 5   1 2 3 7 8    1 3 8    2 3 7 8 9          |
|     ------------------------------------------------------------------------|
|         2     2       2     3      4 2 2       3     2     2 2 2 2          |
|     4x x x x x x  + 3x x x x x  - x x x  + 3x x x x x  - 2x x x x ,         |
|       1 2 3 7 8 9     1 2 3 8 9    2 7 9     1 2 7 8 9     1 2 8 9          |
|     ------------------------------------------------------------------------|
|          2   3    2 2   3    3                 2     2       2       2      |
|     x x x x x  - x x x x  + x x x x x x  - 4x x x x x x  + 3x x x x x x  -  |
|      1 2 3 4 8    1 3 5 8    2 3 4 7 8 9     1 2 3 4 8 9     1 2 3 5 8 9    |
|     ------------------------------------------------------------------------|
|      4     2       3     2     2 2     2   4 2 2       3   2      2 2 2 2   |
|     x x x x  + 3x x x x x  - 2x x x x x , x x x  - 4x x x x x  - x x x x  + |
|      2 4 7 9     1 2 4 8 9     1 2 5 8 9   3 7 8     1 3 7 8 9    2 3 7 9   |
|     ------------------------------------------------------------------------|
|           2     2     2 2 2 2     2       3    2 2 4     3 2 2    2 2 3     |
|     4x x x x x x  + 3x x x x  - 4x x x x x  + x x x , x x x x  + x x x x  - |
|       1 2 3 7 8 9     1 3 8 9     1 2 3 8 9    1 2 9   1 3 7 8    2 3 7 9   |
|     ------------------------------------------------------------------------|
|       2 2   2         2   2 2     2         2     3   2 2     2 2   3       |
|     4x x x x x  - 4x x x x x  + 4x x x x x x  + 3x x x x  + 3x x x x  -     |
|       1 3 7 8 9     1 2 3 7 9     1 2 3 7 8 9     1 3 8 9     1 2 7 9       |
|     ------------------------------------------------------------------------|
|       3     3   4     2       3   2      2 2     2         2     2          |
|     4x x x x , x x x x  - 4x x x x x  - x x x x x  + 4x x x x x x  +        |
|       1 2 8 9   3 4 7 8     1 3 4 8 9    2 3 4 7 9     1 2 3 4 8 9          |
|     ------------------------------------------------------------------------|
|       2 2     2     2       3    2 2   3     3     2    2 2   2             |
|     3x x x x x  - 4x x x x x  + x x x x , x x x x x  + x x x x x  -         |
|       1 3 5 8 9     1 2 3 5 9    1 2 6 9   1 3 4 7 8    2 3 4 7 9           |
|     ------------------------------------------------------------------------|
|       2 2   2         2       2     2         2     3       2     2 2   3   |
|     4x x x x x  - 4x x x x x x  + 4x x x x x x  + 3x x x x x  + 3x x x x  - |
|       1 3 4 8 9     1 2 3 4 7 9     1 2 3 4 8 9     1 3 5 8 9     1 2 4 9   |
|     ------------------------------------------------------------------------|
|       3     3       2     2    2 2 2 2    3   2           2                 |
|     4x x x x , x x x x x x  - x x x x  + x x x x x  - 4x x x x x x x  +     |
|       1 2 5 9   1 2 3 4 5 8    1 3 5 8    2 3 4 8 9     1 2 3 4 5 8 9       |
|     ------------------------------------------------------------------------|
|       2     2        4 2 2       3     2     2 2 2 2   4 2 2       3        |
|     3x x x x x x  - x x x  + 3x x x x x  - 2x x x x , x x x  - 4x x x x x x |
|       1 2 3 5 8 9    2 4 9     1 2 4 5 9     1 2 5 9   3 4 8     1 3 4 5 8 9|
|     ------------------------------------------------------------------------|
|        2 2 2 2         2     2     2 2 2 2     2         2    2 2 2 2       |
|     - x x x x  + 4x x x x x x  + 3x x x x  - 4x x x x x x  + x x x x ,      |
|        2 3 4 9     1 2 3 4 5 9     1 3 5 9     1 2 3 5 6 9    1 2 6 9       |
|     ------------------------------------------------------------------------|
|        3 2 2    2 2 2         2 2               2   2 2     2         2     |
|     x x x x  + x x x x x  - 4x x x x x x  - 4x x x x x  + 4x x x x x x  +   |
|      1 3 4 8    2 3 4 7 9     1 3 4 5 8 9     1 2 3 4 9     1 2 3 4 5 9     |
|     ------------------------------------------------------------------------|
|       3   2 2     2 2     2     3       2       2   2      2 2 3            |
|     3x x x x  + 3x x x x x  - 4x x x x x , x x x x x x  - x x x x  +        |
|       1 3 5 9     1 2 4 6 9     1 2 5 6 9   1 2 3 4 5 8    1 3 5 8          |
|     ------------------------------------------------------------------------|
|      3   2           2     2       2     3      4 2           3             |
|     x x x x x  - 4x x x x x x  + 3x x x x x  - x x x x  + 3x x x x x x  -   |
|      2 3 4 5 9     1 2 3 4 5 9     1 2 3 5 9    2 4 6 9     1 2 4 5 6 9     |
|     ------------------------------------------------------------------------|
|       2 2 2       4 2           3   2      2 2 2             2              |
|     2x x x x x , x x x x  - 4x x x x x  - x x x x x  + 4x x x x x x x  +    |
|       1 2 5 6 9   3 4 5 8     1 3 4 5 9    2 3 4 6 9     1 2 3 4 5 6 9      |
|     ------------------------------------------------------------------------|
|       2 2 2         2       2      2 2 3       3 2        2 2 3             |
|     3x x x x x  - 4x x x x x x  + x x x x , x x x x x  + x x x x  -         |
|       1 3 5 6 9     1 2 3 5 6 9    1 2 6 9   1 3 4 5 8    2 3 4 9           |
|     ------------------------------------------------------------------------|
|       2 2   2         2   2         2                 3   2                 |
|     4x x x x x  - 4x x x x x x  + 4x x x x x x x  + 3x x x x x  +           |
|       1 3 4 5 9     1 2 3 4 6 9     1 2 3 4 5 6 9     1 3 5 6 9             |
|     ------------------------------------------------------------------------|
|       2 2   2       3     2     2 2 4    2 2 2 2       2   3                |
|     3x x x x x  - 4x x x x x , x x x  - x x x x  - 4x x x x x  +            |
|       1 2 4 6 9     1 2 5 6 9   2 3 7    1 3 7 8     1 2 3 7 9              |
|     ------------------------------------------------------------------------|
|       2     2         2 2 2 2     3       2    4 2 2   2 2   3    2 2     2 |
|     4x x x x x x  + 3x x x x  - 4x x x x x  + x x x , x x x x  - x x x x x  |
|       1 2 3 7 8 9     1 2 7 9     1 2 7 8 9    1 8 9   2 3 4 7    1 3 4 7 8 |
|     ------------------------------------------------------------------------|
|           2     2       2                 2 2     2     3       2           |
|     - 4x x x x x x  + 4x x x x x x x  + 3x x x x x  - 4x x x x x  +         |
|         1 2 3 4 7 9     1 2 3 4 7 8 9     1 2 4 7 9     1 2 4 8 9           |
|     ------------------------------------------------------------------------|
|      4     2   2 2 2 2    2 2 2 2       2   2         2     2               |
|     x x x x , x x x x  - x x x x  - 4x x x x x x  + 4x x x x x x  +         |
|      1 5 8 9   2 3 4 7    1 3 4 8     1 2 3 4 7 9     1 2 3 4 8 9           |
|     ------------------------------------------------------------------------|
|       2 2 2 2     3       2    4 2 2   2 2 3      2 2 2           2   3     |
|     3x x x x  - 4x x x x x  + x x x , x x x x  - x x x x x  - 4x x x x x  + |
|       1 2 4 9     1 2 4 5 9    1 5 9   2 3 4 7    1 3 4 5 8     1 2 3 4 9   |
|     ------------------------------------------------------------------------|
|       2     2         2 2 2         3              4 2           2   3      |
|     4x x x x x x  + 3x x x x x  - 4x x x x x x  + x x x x , x x x x x  -    |
|       1 2 3 4 5 9     1 2 4 6 9     1 2 4 5 6 9    1 5 6 9   1 2 3 4 5      |
|     ------------------------------------------------------------------------|
|      2 2 4    3   2           2     2       2     3      4 2 2       3     2|
|     x x x  + x x x x x  - 4x x x x x x  + 3x x x x x  - x x x  + 3x x x x x |
|      1 3 5    2 3 4 5 6     1 2 3 4 5 6     1 2 3 5 6    2 4 6     1 2 4 5 6|
|     ------------------------------------------------------------------------|
|         2 2 2 2   4 2 2       3   2      2 2 2 2         2     2     2 2 2 2|
|     - 2x x x x , x x x  - 4x x x x x  - x x x x  + 4x x x x x x  + 3x x x x |
|         1 2 5 6   3 4 5     1 3 4 5 6    2 3 4 6     1 2 3 4 5 6     1 3 5 6|
|     ------------------------------------------------------------------------|
|         2       3    2 2 4     3 2 2    2 2 3       2 2   2         2   2 2 |
|     - 4x x x x x  + x x x , x x x x  + x x x x  - 4x x x x x  - 4x x x x x  |
|         1 2 3 5 6    1 2 6   1 3 4 5    2 3 4 6     1 3 4 5 6     1 2 3 4 6 |
|     ------------------------------------------------------------------------|
|         2         2     3   2 2     2 2   3     3     3   2 2 4    2 2 2 2  |
|     + 4x x x x x x  + 3x x x x  + 3x x x x  - 4x x x x , x x x  - x x x x  -|
|         1 2 3 4 5 6     1 3 5 6     1 2 4 6     1 2 5 6   2 3 4    1 3 4 5  |
|     ------------------------------------------------------------------------|
|         2   3       2     2         2 2 2 2     3       2    4 2 2          |
|     4x x x x x  + 4x x x x x x  + 3x x x x  - 4x x x x x  + x x x )         |
|       1 2 3 4 6     1 2 3 4 5 6     1 2 4 6     1 2 4 5 6    1 5 6          |
|                                                                             |
|o5 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i6 : isPrimary Q                                                             |
|                                                                             |
|o6 = true                                                                    |
+-----------------------------------------------------------------------------+

Ways to use joinIdeals:
=======================

  * "joinIdeals(Ideal,Ideal)"

For the programmer
==================

The object joinIdeals, is a method function (Macaulay2Doc).

localHilbertRegularity -- regularity of the local Hilbert function of a polynomial ideal
****************************************************************************************

Synopsis
========

  * Usage:
        d = localHilbertRegularity(p,I)
  * Inputs:
      * p, a point (NAGtypes),
      * I, an ideal (Macaulay2Doc), or a one-row matrix (Macaulay2Doc),
        of generators
  * Optional inputs (Macaulay2Doc),:
      * Tolerance, => ..., default
        value null, optional argument for numerical tolerance
  * Outputs:
      * d, an integer (Macaulay2Doc),

Description
===========

The g-corners, of the ideal are computed in order to find the
Hilbert polynomial, which is compared to the Hilbert function to find the
degree of regularity, which is the degree at which the two become equal.

+----------------------------------------------------------------+
|i1 : R = CC[x,y];                                               |
+----------------------------------------------------------------+
|i2 : I = ideal{x^2,x*y}                                         |
|                                                                |
|             2                                                  |
|o2 = ideal (x , x*y)                                            |
|                                                                |
|o2 : Ideal of R                                                 |
+----------------------------------------------------------------+
|i3 : d = localHilbertRegularity(origin R, I)                    |
|-- warning: experimental computation over inexact field begun   |
|--          results not reliable (one warning given per session)|
|                                                                |
|o3 = 2                                                          |
+----------------------------------------------------------------+
|i4 : D = truncatedDual(origin R, I, 3)                          |
|                                                                |
|o4 = | 1 y x .5y2 .166667y3 |                                   |
|                                                                |
|o4 : DualSpace                                                  |
+----------------------------------------------------------------+
|i5 : L = hilbertFunction({0,1,2,3}, D)                          |
|                                                                |
|o5 = {1, 2, 1, 1}                                               |
|                                                                |
|o5 : List                                                       |
+----------------------------------------------------------------+

See also gCorners.

Ways to use localHilbertRegularity:
===================================

  * "localHilbertRegularity(AbstractPoint,Ideal)"
  * "localHilbertRegularity(AbstractPoint,Matrix)"

For the programmer
==================

The object localHilbertRegularity, is a
method function with options (Macaulay2Doc).

mapToPunctualHilbertScheme -- maps an ideal into a point in a certain punctual Hilbert scheme
*********************************************************************************************

Synopsis
========

  * Usage:
        I = mapToPunctualHilbertScheme(Q)
  * Inputs:
      * Q, an ideal (Macaulay2Doc), a primary ideal
  * Outputs:
      * I, an ideal (Macaulay2Doc), an ideal that parametrizes Q
        in a punctual Hilbert scheme

Description
===========

This method maps a P-primary ideal Q into a point in a punctual Hilbert scheme.
Let $\mathbb{K}$ be a field of characteristic zero and let $P$ be a prime ideal
of codimension $c$ in the polynomial ring $R = \mathbb{K}[x_1,\ldots,x_n]$. We
write $\mathbb{F}$ for the field of fractions of the integral domain $R/P$. To
simplify our notation, perhaps after a linear change of coordinates, we assume
that $\{ x_{c+1}, \ldots, x_n \}$ is a maximal independent set of variables
modulo $P$.

The main purpose of this method is to reduce the study of arbitrary $P$-primary
ideals in $R = \mathbb{K}[x_1,\ldots,x_n]$ to a zero-dimensional setting over
the function field $\mathbb{F}$. This reduction is made by parametrizing
$P$-primary ideals with the punctual Hilbert scheme $ {\rm Hilb}^m (
\,\mathbb{F}[[y_1,\ldots,y_c]] \,). $ This is a quasiprojective scheme over the
function field $\mathbb{F}$. Its classical points are   ideals of colength $m$
in the local ring $\mathbb{F}[[y_1,\ldots,y_c]]$.

This method maps a P-primary ideal Q into a unique point in ${\rm Hilb}^m (
\,\mathbb{F}[[y_1,\ldots,y_c]] \,)$ that corresponds with Q. This method can be
seen as an implementation of the map $\gamma$ described in Section 2 of the
paper Primary ideals and their differential equations (see
https://arxiv.org/abs/2001.04700 ).

+-----------------------------------------+
|i1 : R = QQ[x_1, x_2, x_3]               |
|                                         |
|o1 = R                                   |
|                                         |
|o1 : PolynomialRing                      |
+-----------------------------------------+
|i2 : Q = ideal(x_1^2, x_2^2, x_1-x_2*x_3)|
|                                         |
|             2   2                       |
|o2 = ideal (x , x , - x x  + x )         |
|             1   2     2 3    1          |
|                                         |
|o2 : Ideal of R                          |
+-----------------------------------------+
|i3 : mapToPunctualHilbertScheme Q        |
|                                         |
|                           2             |
|o3 = ideal (hx  - x hx , hx )            |
|              1    3  2    2             |
|                                         |
|                  /    R   \             |
|o3 : Ideal of frac|--------|[hx ..hx ]   |
|                  |(x , x )|   1    2    |
|                  \  2   1 /             |
+-----------------------------------------+

See also
========

  * Strategy => "PunctualQuot" -- strategy
    for computing Noetherian operators

Ways to use mapToPunctualHilbertScheme:
=======================================

  * "mapToPunctualHilbertScheme(Ideal)"

For the programmer
==================

The object mapToPunctualHilbertScheme, is a
method function (Macaulay2Doc).

noetherianOperators -- Noetherian operators
*******************************************

Description
===========

Let $R$ be a polynomial ring $R = K[x_1,\ldots,x_n]$ over a field $K$ of
characteristic zero. Consider the Weyl algebra $D = R<dx_1,\ldots,dx_n>$, a
prime ideal $P \subset R$ and a $P$-primary ideal. When this method is applied
we obtain a finite list of differential operators $L_1,\ldots,L_m \in D$ such
that $$ Q = \{ f \,\in\, R\, \mid\, L_i\, \bullet\, f\, \in P, \ \forall  1 \le
i \le m \}. $$ We say that $\{L_1,\ldots,L_m\}$ is a set of Noetherian
operators for the primary ideal $Q$. In the output of the algorithm we always
have that $m$ (the number of Noetherian operators) is equal to the multiplicity
of $Q$ over the prime ideal $P$.

+-----------------------------------------------------------------------------+
|i1 : R=QQ[x_1,x_2,x_3,x_4]                                                   |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : Q = ideal(x_1^2,x_1*x_2,x_1*x_3,x_1*x_4-x_3^2+x_1,x_3^2*x_4-x_2^2,x_3^2*|
|                                                                             |
|             2                 2               2      2   2             2    |
|o2 = ideal (x , x x , x x , - x  + x x  + x , x x  - x , x x  - x x  - x  +  |
|             1   1 2   1 3     3    1 4    1   3 4    2   3 4    2 3    3    |
|     ------------------------------------------------------------------------|
|     2x )                                                                    |
|       1                                                                     |
|                                                                             |
|o2 : Ideal of R                                                              |
|-----------------------------------------------------------------------------|
|x_4-x_3^2-x_2*x_3+2*x_1)                                                     |
+-----------------------------------------------------------------------------+
|i3 : isPrimary Q                                                             |
|                                                                             |
|o3 = true                                                                    |
+-----------------------------------------------------------------------------+
|i4 : noetherianOperators(Q, Strategy => "PunctualQuot")                      |
|                                                                             |
|o4 = {| 1 |, | dx_2 |, | dx_3 |, |                                           |
|     ------------------------------------------------------------------------|
|     (x_4^2+x_4)dx_2^2+(2x_4^2+2)dx_2dx_3+(x_4+1)dx_3^2+2dx_1 |}             |
|                                                                             |
|o4 : List                                                                    |
+-----------------------------------------------------------------------------+

Ways to use noetherianOperators:
================================

  * noetherianOperators(Ideal), --
    Noetherian operators of a primary ideal
  * noetherianOperators(Ideal,Ideal)
    -- Noetherian operators of a
    primary component
  * noetherianOperators(Module), --
    Noetherian operators of a primary submodule
  * noetherianOperators(Module,Ideal)
    -- Noetherian operators of a
    primary component

For the programmer
==================

The object noetherianOperators, is a method function with options (Macaulay2Doc).

noetherianOperators(Ideal) -- Noetherian operators of a primary ideal
*********************************************************************

Synopsis
========

  * Function: noetherianOperators,
  * Usage:
        noetherianOperators Q
        noetherianOperators (Q, Strategy => "PunctualQuot")
  * Inputs:
      * Q, an ideal (Macaulay2Doc), assumed to be primary
  * Outputs:
      * a list (Macaulay2Doc), of differential operators,

Description
===========

Compute a set of Noetherian operators for the primary ideal I.

+----------------------------------------------------+
|i1 : R = QQ[x,y,t];                                 |
+----------------------------------------------------+
|i2 : I = ideal(x^2, y^2-x*t);                       |
|                                                    |
|o2 : Ideal of R                                     |
+----------------------------------------------------+
|i3 : noetherianOperators I                          |
|                                                    |
|o3 = {| 1 |, | dy |, | tdy^2+2dx |, | tdy^3+6dxdy |}|
|                                                    |
|o3 : List                                           |
+----------------------------------------------------+

The optional argument Strategy can be used to choose different algorithms. Each
strategy may accept additional optional arguments, see the documentation page
for each strategy for details.

  * "PunctualQuot" (default)
  * "PunctualHilbert" (alias for "PunctualQuot")
  * "MacaulayMatrix"
  * "Hybrid"

Caveat
======

The behavior is undefined if Q is not primary. For non-primary ideals, use
noetherianOperators(Ideal,Ideal),

Ways to use this method:
========================

  * noetherianOperators(Ideal), --
    Noetherian operators of a primary ideal

noetherianOperators(Ideal,Ideal) -- Noetherian operators of a primary component
*******************************************************************************

Synopsis
========

  * Function: noetherianOperators,
  * Usage:
        noetherianOperators (I, P)
        noetherianOperators (I, P, Strategy => "MacaulayMatrix")
        noetherianOperators (I, P, Rational => false)
  * Inputs:
      * I, an ideal (Macaulay2Doc),
      * P, an ideal (Macaulay2Doc), a minimal prime of $I$
  * Outputs:
      * a list (Macaulay2Doc), of differential operators,

Description
===========

Compute a set of Noetherian operators for the $P$-primary component of $I$.

+------------------------------------------------------------------+
|i1 : R = QQ[x,y,t];                                               |
+------------------------------------------------------------------+
|i2 : I1 = ideal(x^2, y^2-x*t);                                    |
|                                                                  |
|o2 : Ideal of R                                                   |
+------------------------------------------------------------------+
|i3 : I2 = ideal((x-t)^2);                                         |
|                                                                  |
|o3 : Ideal of R                                                   |
+------------------------------------------------------------------+
|i4 : I = intersect(I1, I2);                                       |
|                                                                  |
|o4 : Ideal of R                                                   |
+------------------------------------------------------------------+
|i5 : noetherianOperators(I, radical I1)                           |
|                                                                  |
|o5 = {| 1 |, | dy |, | tdy^2+2dx |, | tdy^3+6dxdy |}              |
|                                                                  |
|o5 : List                                                         |
+------------------------------------------------------------------+
|i6 : noetherianOperators(I, radical I2) == noetherianOperators(I2)|
|                                                                  |
|o6 = true                                                         |
+------------------------------------------------------------------+

The optional argument Strategy can be used to choose different algorithms. Each
strategy may accept additional optional arguments, see the documentation page
for each strategy for details.

  * "MacaulayMatrix" (default)
  * "Hybrid"

If the prime $P$ is known to be a rational point, the optional argument
Rational can be set to true.  This may offer a speed-up in the computation.

See also
========

  * noetherianOperators -- Noetherian operators
  * noetherianOperators(Ideal), --
    Noetherian operators of a primary ideal
  * noetherianOperators(Module,Ideal)
    -- Noetherian operators of a
    primary component

Ways to use this method:
========================

  * noetherianOperators(Ideal,Ideal)
    -- Noetherian operators of a
    primary component

noetherianOperators(Module) -- Noetherian operators of a primary submodule
**************************************************************************

Synopsis
========

  * Function: noetherianOperators,
  * Usage:
        noetherianOperators U
  * Inputs:
      * U, a module (Macaulay2Doc), a primary submodule
  * Outputs:
      * a list (Macaulay2Doc), of differential operators,

Description
===========

Compute a set of Noetherian operators for the primary submodule U.

This method contains an implementation of Algorithm 4.1 in the paper Primary
decomposition of modules: a computational differential approach (see
https://arxiv.org/abs/2104.03385 ). For more details, see Section 4 of the
paper Primary decomposition of modules: a computational differential approach
(see https://arxiv.org/abs/2104.03385 ).

+--------------------------------------------------------------------+
|i1 : R = QQ[x_1,x_2,x_3]                                            |
|                                                                    |
|o1 = R                                                              |
|                                                                    |
|o1 : PolynomialRing                                                 |
+--------------------------------------------------------------------+
|i2 : U = image matrix {{x_1, x_2^2, 0}, {x_3, x_3^2, x_2^2-x_1*x_3}}|
|                                                                    |
|o2 = image | x_1 x_2^2 0            |                               |
|           | x_3 x_3^2 x_2^2-x_1x_3 |                               |
|                                                                    |
|                             2                                      |
|o2 : R-module, submodule of R                                       |
+--------------------------------------------------------------------+
|i3 : noetherianOperators U                                          |
|                                                                    |
|o3 = {| -x_3 |}                                                     |
|      |  x_1 |                                                      |
|                                                                    |
|o3 : List                                                           |
+--------------------------------------------------------------------+

Ways to use this method:
========================

  * noetherianOperators(Module), --
    Noetherian operators of a primary submodule

noetherianOperators(Module,Ideal) -- Noetherian operators of a primary component
********************************************************************************

Synopsis
========

  * Function: noetherianOperators,
  * Usage:
        noetherianOperators (U, P)
  * Inputs:
      * U, a module (Macaulay2Doc),
      * P, an ideal (Macaulay2Doc), a minimal prime of comodule U
  * Outputs:
      * a list (Macaulay2Doc), of differential operators,

Description
===========

Compute a set of Noetherian operators for the $P$-primary component of comodule
U.

+------------------------------------------+
|i1 : R = QQ[x,y,z]                        |
|                                          |
|o1 = R                                    |
|                                          |
|o1 : PolynomialRing                       |
+------------------------------------------+
|i2 : U = image matrix{{x,y,z},{y,z,x}}    |
|                                          |
|o2 = image | x y z |                      |
|           | y z x |                      |
|                                          |
|                             2            |
|o2 : R-module, submodule of R             |
+------------------------------------------+
|i3 : P = first associatedPrimes comodule U|
|                                          |
|o3 = ideal (y - z, x - z)                 |
|                                          |
|o3 : Ideal of R                           |
+------------------------------------------+
|i4 : noetherianOperators(U, P)            |
|                                          |
|o4 = {| -1 |}                             |
|      |  1 |                              |
|                                          |
|o4 : List                                 |
+------------------------------------------+

If there are no embedFded primes, running this command for all associated
primes is equivalent to running a differential primary decomposition.

+-----------------------------------------------------------------------------+
|i5 : associatedPrimes comodule U /                                           |
|         (P -> {P, noetherianOperators(U,P)})                                |
|                                                                             |
|                                                            2          2     |
|o5 = {{ideal (y - z, x - z), {| -1 |}}, {ideal (x + y + z, y  + y*z + z ), {||
|                              |  1 |                                        ||
|     ------------------------------------------------------------------------|
|     -z |}}}                                                                 |
|      y |                                                                    |
|                                                                             |
|o5 : List                                                                    |
+-----------------------------------------------------------------------------+
|i6 : netList differentialPrimaryDecomposition U                              |
|                                                                             |
|     +--------------------------------+--------+                             |
|o6 = |ideal (y - z, x - z)            |{| -1 |}|                             |
|     |                                | |  1 | |                             |
|     +--------------------------------+--------+                             |
|     |                   2          2 |        |                             |
|     |ideal (x + y + z, y  + y*z + z )|{| -z |}|                             |
|     |                                | |  y | |                             |
|     +--------------------------------+--------+                             |
+-----------------------------------------------------------------------------+

See also
========

  * noetherianOperators -- Noetherian operators
  * noetherianOperators(Ideal,Ideal)
    -- Noetherian operators of a
    primary component
  * differentialPrimaryDecomposition,
    -- compute a differential primary decomposition

Ways to use this method:
========================

  * noetherianOperators(Module,Ideal)
    -- Noetherian operators of a
    primary component

noethOpsFromComponents -- merge Noetherian operators for non-primary ideals
***************************************************************************

Synopsis
========

  * Usage:
        noethOpsFromComponents L
  * Inputs:
      * L, a list (Macaulay2Doc), of ordered pairs (P, N) where P
        is a minimal prime of I, and N is a set of Noetherian operators for the
        P-primary component of I
  * Outputs:
      * a list (Macaulay2Doc), of differential operators,

Description
===========

Let $I$ be an unmixed ideal in a polynomial ring $R = K[x_1, ..., x_n]$, with
primary decomposition $I = Q_1 \cap ... \cap Q_s$, where $Q_i$ is
$P_i$-primary. If $N_i$ is a set of Noetherian operators for $Q_i$, then one
can construct a set of differential operators $N$ for $I$ which satisfies the
Noetherian operator condition: given $f \in R$, one has $f \in I$ iff $D(f)
\in\sqrt{I}$ for all $D \in N$.

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x,y,t]                                                           |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : I = intersect(ideal((y+t)^2), ideal(x^2, y^2 - t*x))                    |
|                                                                             |
|             4      2      3          2    2 2      3   2 2     2       2 2  |
|o2 = ideal (y  - x*y t + 2y t - 2x*y*t  + y t  - x*t , x y  + 2x y*t + x t ) |
|                                                                             |
|o2 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i3 : radI = radical I                                                        |
|                                                                             |
|             2                                                               |
|o3 = ideal (y  + y*t, x*y + x*t)                                             |
|                                                                             |
|o3 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i4 : primes = associatedPrimes I                                             |
|                                                                             |
|o4 = {ideal(y + t), ideal (y, x)}                                            |
|                                                                             |
|o4 : List                                                                    |
+-----------------------------------------------------------------------------+
|i5 : L = primes / (P -> (P, noetherianOperators(I, P)))                      |
|                                                                             |
|o5 = {(ideal(y + t), {| 1 |, | dy |}), (ideal (y, x), {| 1 |, | dy |, |      |
|     ------------------------------------------------------------------------|
|     tdy^2+2dx |, | tdy^3+6dxdy |})}                                         |
|                                                                             |
|o5 : List                                                                    |
+-----------------------------------------------------------------------------+
|i6 : N = noethOpsFromComponents L                                            |
|                                                                             |
|o6 = {| 1 |, | dy |, | (yt+t2)dy^2+(2y+2t)dx |, | (yt+t2)dy^3+(6y+6t)dxdy |} |
|                                                                             |
|o6 : List                                                                    |
+-----------------------------------------------------------------------------+
|i7 : all(flatten table(N, I_*, (D, f) -> (D f) % radI == 0), identity)       |
|                                                                             |
|o7 = true                                                                    |
+-----------------------------------------------------------------------------+

Note that this construction justifies the focus of Noetherian operators on the
case that the ideal I is primary: in order to get a useful membership test for
a non-primary (but still unmixed) ideal, it suffices to compute Noetherian
operators on each primary component, and then combine them in the way given
above.

See also
========

  * noetherianOperators(Ideal,Ideal)
    -- Noetherian operators of a
    primary component

Ways to use noethOpsFromComponents:
===================================

  * "noethOpsFromComponents(List)"

For the programmer
==================

The object noethOpsFromComponents, is a
method function (Macaulay2Doc).

normalize -- rescale a differential operator
********************************************

Synopsis
========

  * Usage:
        normalize D
  * Inputs:
      * D, a differential operator,
  * Outputs:
      * a differential operator,

Description
===========

Rescales a differential operator so that the leading term of the leading
coefficient is 1.

+-------------------------------------------+
|i1 : S = diffOpRing (QQ[x,y,t]);           |
+-------------------------------------------+
|i2 : D = diffOp(3*x^3*dx^2*dt + (x+y)*dt^2)|
|                                           |
|o2 = | 3x3dx^2dt+(x+y)dt^2 |               |
|                                           |
|                1                          |
|o2 : DiffOp in S                           |
+-------------------------------------------+
|i3 : normalize D                           |
|                                           |
|o3 = | x3dx^2dt+(1/3x+1/3y)dt^2 |          |
|                                           |
|                1                          |
|o3 : DiffOp in S                           |
+-------------------------------------------+

This can be useful when computing "canonical" sets of Noetherian operators, as
a valid set of Noetherian operators stays valid even after rescaling.

+-----------------------------------------------------------------------------+
|i4 : I = ideal(x^2,y^2 - x*t);                                               |
|                                                                             |
|o4 : Ideal of QQ[x..y, t]                                                    |
+-----------------------------------------------------------------------------+
|i5 : nops = noetherianOperators(I, Strategy => "MacaulayMatrix");            |
+-----------------------------------------------------------------------------+
|i6 : nops / normalize == {diffOp 1_S, diffOp dy, diffOp(t*dy^2 + 2*dx),      |
|                                                                             |
|o6 = true                                                                    |
|-----------------------------------------------------------------------------|
|diffOp(t*dy^3 + 6*dx*dy)}                                                    |
+-----------------------------------------------------------------------------+

Ways to use normalize:
======================

  * "normalize(DiffOp)"

For the programmer
==================

The object normalize, is a method function (Macaulay2Doc).

numericalNoetherianOperators -- Noetherian operators via numerical interpolation
********************************************************************************

Synopsis
========

  * Usage:
        numericalNoetherianOperators(I)
  * Inputs:
      * I, an ideal (Macaulay2Doc), unmixed
  * Outputs:
      * a list (Macaulay2Doc), of expressions (Macaulay2Doc),

Description
===========

The method computes specialized Noetherian operators from many sampled points,
and attempts to find fitting rational functions using rational function
interpolation.

+-----------------------------------------------------------+
|i1 : R = CC[x,y,t];                                        |
+-----------------------------------------------------------+
|i2 : I = ideal(x^2, y^2 - x*t);                            |
|                                                           |
|o2 : Ideal of R                                            |
+-----------------------------------------------------------+
|i3 : numericalNoetherianOperators(I, DependentSet => {x,y})|
|                                                           |
|                 2    1      t   3                         |
|o3 = {1, 1*dy, dy  + ---*dx, -*dy  + 1*dx*dy}              |
|                     .5t     6                             |
|                                                           |
|o3 : List                                                  |
+-----------------------------------------------------------+

The behavior of the function can be adjusted using options. Currently only the
option DependentSet, is required. The following are
supported:

TrustedPoint => a point (NAGtypes). The function does not
compute specialized Noetherian operators from scratch for each point. Instead,
it computes it for a "trusted" point on the variety, and uses the obtained
Noetherian operators as a template for the rest of the computation. If
TrustedPoint is unset, the first point returned by the sampler will be used as
the trusted point.

NoetherianDegreeLimit => a non-negative integer (Macaulay2Doc).
Limits the degrees of the Noetherian operators (with respect to the variables
$dx_i$). If unset, will compute the full Noetherian operators of the "trusted"
point. Can introduce speed-ups when the maximal degree of the Noetherian
operators is known in advance.

Tolerance => a positive real number. This specifies the numerical precision
when computing the specialized Noetherian operators. The default value is 1e-6.
See Tolerance (NoetherianOperators).

Sampler => a function, taking inputs (n,I), where I is an ideal (Macaulay2Doc),
and n is an integer. The sampler function returns a list
of n points (NAGtypes), on the component of interest of I.
If unset, the default sampler uses Bertini, and assumes
that I is primary.

DependentSet => a list (Macaulay2Doc), of variables that are
algebraically dependent. See DependentSet, for details.

InterpolationTolerance => a positive real number. This specifies the numerical
precision for the interpolation routines. The default value is 1e-6. See
Tolerance (NoetherianOperators).

InterpolationDegreeLimit => a non-negative integer (Macaulay2Doc).
Limits the degree of the interpolated rational function coefficients. If no
rational functions are found within the degree limit, outputs an incomplete
differential operator

+-----------------------------------------------------------------------------+
|i1 : R = CC[x,y,t];                                                          |
+-----------------------------------------------------------------------------+
|i2 : I = ideal(x^2, y^2 - x*t);                                              |
|                                                                             |
|o2 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i3 : numericalNoetherianOperators(I, DependentSet => {x,y},                  |
|                                                                             |
|                   2   ?     ?   3                                           |
|o3 = {1, 1*dy, 1*dy  + -*dx, -*dy  + 1*dx*dy}                                |
|                       ?     ?                                               |
|                                                                             |
|o3 : List                                                                    |
|-----------------------------------------------------------------------------|
|InterpolationDegreeLimit => 0)                                               |
+-----------------------------------------------------------------------------+

See also
========

  * rationalInterpolation -- numerically
    interpolate rational functions

Ways to use numericalNoetherianOperators:
=========================================

  * "numericalNoetherianOperators(Ideal)"

For the programmer
==================

The object numericalNoetherianOperators, is
a method function with options (Macaulay2Doc).

orthogonalInSubspace -- Orthogonal of a space
*********************************************

Synopsis
========

  * Usage:
        S = orthogonalInSubspace(D, T, tol)
  * Inputs:
      * D, an instance of the type DualSpace (NAGtypes), or an
        instance of the type PolySpace (NAGtypes), a space of
        which to find the orthogonal
      * T, an instance of the type PolySpace (NAGtypes),
        ambient space
      * tol, a number (Macaulay2Doc), a positive number, the
        numerical tolerance
  * Outputs:
      * S, an instance of the type PolySpace (NAGtypes),

Description
===========

This functions computes the subspace of the polynomial space T that is
orthogonal to the dual space (or polynomial space) D.

+--------------------------------------------+
|i1 : R = CC[x,y];                           |
+--------------------------------------------+
|i2 : T = polySpace matrix{{1,x,y}};         |
+--------------------------------------------+
|i3 : D = dualSpace(matrix{{x-y}}, origin R);|
+--------------------------------------------+
|i4 : S = orthogonalInSubspace(D, T, 1e-6)   |
|                                            |
|o4 = | 1x+y 1 |                             |
|                                            |
|o4 : PolySpace                              |
+--------------------------------------------+

Ways to use orthogonalInSubspace:
=================================

  * "orthogonalInSubspace(DualSpace,PolySpace,Number)"
  * "orthogonalInSubspace(PolySpace,PolySpace,Number)"

For the programmer
==================

The object orthogonalInSubspace, is a method function (Macaulay2Doc).

pairingMatrix -- Applies dual space functionals to polynomials
**************************************************************

Synopsis
========

  * Usage:
        M = pairingMatrix(S, D)
  * Inputs:
      * S, an instance of the type PolySpace (NAGtypes),
      * D, an instance of the type DualSpace (NAGtypes),
  * Outputs:
      * M, a matrix (Macaulay2Doc), containing the values of the
        generators of D applied to the generators of S

Description
===========

The dual space represents functionals from the polynomial ring to the base
field. Given a polySpace S with n generators f_1,...,f_n and a dualSpace D with
m generators p_1,...,p_m, pairingMatrix returns a nxm matrix M over the base
field whose entries are p_j(f_i).

A dual functional is applied to a polynomial by taking the standard pairing of
their coefficient vectors.  In other words, the functional represented by the
monomial a acts on monomials in the polynomial ring as a(a) = 1 and a(b) = 0
for all other monomials b.

+----------------------------------------------+
|i1 : R = CC[x,y];                             |
+----------------------------------------------+
|i2 : S = polySpace matrix{{x+y,2*x+y^2}};     |
+----------------------------------------------+
|i3 : D = dualSpace(matrix{{1,x,y}}, origin R);|
+----------------------------------------------+
|i4 : M = pairingMatrix(S, D)                  |
|                                              |
|o4 = {-1} | 0 1 1 |                           |
|     {-2} | 0 2 0 |                           |
|                                              |
|             2      3                         |
|o4 : Matrix R  <-- R                          |
+----------------------------------------------+

The function pairingMatrix can also be called with one or both inputs a
ring element (Macaulay2Doc).  If both arguments are single
elements, the output is also a ring element rather than a matrix.

+--------------------------+
|i5 : pairingMatrix(S, 1+x)|
|                          |
|o5 = {-1} | 1 |           |
|     {-2} | 2 |           |
|                          |
|             2      1     |
|o5 : Matrix R  <-- R      |
+--------------------------+
|i6 : pairingMatrix(x, D)  |
|                          |
|o6 = {-1} | 0 1 0 |       |
|                          |
|             1      3     |
|o6 : Matrix R  <-- R      |
+--------------------------+
|i7 : pairingMatrix(x, 1+x)|
|                          |
|o7 = 1                    |
|                          |
|o7 : R                    |
+--------------------------+

Ways to use pairingMatrix:
==========================

  * "pairingMatrix(PolySpace,DualSpace)"
  * "pairingMatrix(PolySpace,PolySpace)"
  * "pairingMatrix(PolySpace,RingElement)"
  * "pairingMatrix(RingElement,DualSpace)"
  * "pairingMatrix(RingElement,RingElement)"

For the programmer
==================

The object pairingMatrix, is a method function (Macaulay2Doc).

rationalInterpolation -- numerically interpolate rational functions
*******************************************************************

Synopsis
========

  * Usage:
        (n,d) = rationalInterpolation(pts, vals, numBasis, denBasis)
        (n,d) = rationalInterpolation(pts, vals, numDenBasis)
  * Inputs:
      * pts, a list (Macaulay2Doc), of one-row matrices
        corresponding to points at which the rational function was evaluated
      * vals, a list (Macaulay2Doc), of numbers (Macaulay2Doc),
        corresponding to evaluations of the rational
        function
      * numBasis, a matrix (Macaulay2Doc), a one-row
        matrix (Macaulay2Doc), of monomials, the monomial support of the
        numerator
      * denBasis, a matrix (Macaulay2Doc), a one-row
        matrix (Macaulay2Doc), of monomials, the monomial support of the
        numerator
      * numDenBasis, a matrix (Macaulay2Doc), a one-row
        matrix (Macaulay2Doc), of monomials, the monomial supports of
        both the numerator and denominator
  * Optional inputs (Macaulay2Doc),:
      * Tolerance, => a real number (Macaulay2Doc),
        default value .000001, default value 1e-6
  * Outputs:
      * n, a ring element (Macaulay2Doc), the numerator of
        the rational function
      * d, a ring element (Macaulay2Doc), the denominator of
        the rational function

Description
===========

Given a list of points $pts = \{p_1,\dots,p_k\}$ and values $vals =
\{v_1,\dots,v_k\}$, attempts to find a rational function $f = g/h$, such that
$f(p_i) = v_i$. The polynomials $g$ and $h$ have monomial support numBasis and
denBasis respectively.

+-----------------------------------------------------------------------------+
|i1 : R = CC[x,y]                                                             |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : pts = {point{{1,0}}, point{{0,1}}, point{{1,1}}, point{{-1,-1}},        |
|                                                                             |
|o2 = {{1, 0}, {0, 1}, {1, 1}, {-1, -1}, {-1, 0}}                             |
|                                                                             |
|o2 : List                                                                    |
|-----------------------------------------------------------------------------|
|point{{-1,0}}}                                                               |
+-----------------------------------------------------------------------------+
|i3 : vals = {1, 0, 1/2, -1/2, -1}                                            |
|                                                                             |
|            1    1                                                           |
|o3 = {1, 0, -, - -, -1}                                                      |
|            2    2                                                           |
|                                                                             |
|o3 : List                                                                    |
+-----------------------------------------------------------------------------+
|i4 : numBasis = matrix{{x,y}}                                                |
|                                                                             |
|o4 = | x y |                                                                 |
|                                                                             |
|             1      2                                                        |
|o4 : Matrix R  <-- R                                                         |
+-----------------------------------------------------------------------------+
|i5 : denBasis = matrix{{x^2,y^2}}                                            |
|                                                                             |
|o5 = | x2 y2 |                                                               |
|                                                                             |
|             1      2                                                        |
|o5 : Matrix R  <-- R                                                         |
+-----------------------------------------------------------------------------+
|i6 : rationalInterpolation(pts, vals, numBasis, denBasis)                    |
|                                                                             |
|            2    2                                                           |
|o6 = (1x, 1x  + y )                                                          |
|                                                                             |
|o6 : Sequence                                                                |
+-----------------------------------------------------------------------------+

The output corresponds to the function $x / (x^2 + y^2)$. If no fitting
rational function is found, the method returns an error.

The method rationalInterpolation(List,List,Ring),
can be used to choose monomial
supports automatically.

Caveat
======

The method uses the first point to remove the rational functions whose
numerator and denominator would both evaluate to 0 on the point. Because of
this, the first entry of val should be non-zero.

See also
========

  * rationalInterpolation(List,List,Ring)
    -- numerically interpolate
    rational functions

Ways to use rationalInterpolation:
==================================

  * "rationalInterpolation(List,List,Matrix)"
  * "rationalInterpolation(List,List,Matrix,Matrix)"
  * rationalInterpolation(List,List,Ring)
    -- numerically interpolate
    rational functions

For the programmer
==================

The object rationalInterpolation, is a
method function with options (Macaulay2Doc).

rationalInterpolation(List,List,Ring) -- numerically interpolate rational functions
***********************************************************************************

Synopsis
========

  * Function: rationalInterpolation,
  * Usage:
        (n,d) = rationalInterpolation(pts, vals, R)
  * Inputs:
      * pts, a list (Macaulay2Doc), of one-row matrices
        corresponding to points at which the rational function was evaluated
      * vals, a list (Macaulay2Doc), of numbers (Macaulay2Doc),
        corresponding to evaluations of the rational
        function
      * R, a ring (Macaulay2Doc), the polynomial ring in which the
        numerator and denominator are sought
  * Optional inputs (Macaulay2Doc),:
      * Tolerance, => ..., default
        value .000001, optional argument for numerical tolerance
  * Outputs:
      * n, a ring element (Macaulay2Doc), the numerator of
        the rational function
      * d, a ring element (Macaulay2Doc), the denominator of
        the rational function

Description
===========

Given a list of points $pts = \{p_1,\dots,p_k\}$ and values $vals =
\{v_1,\dots,v_k\}$, attempts to find a rational function $f = g/h$, such that
$f(p_i) = v_i$. The method first tries to find polynomials $g,h$ of degree 0;
if this fails, it tries to find $g,h$ of degree 1 and so on. This procedure
stops when there are not enough points to compute the next degree, in which
case an error will be thrown.

+---------------------------------------------------------------------+
|i1 : R = CC[x]                                                       |
|                                                                     |
|o1 = R                                                               |
|                                                                     |
|o1 : PolynomialRing                                                  |
+---------------------------------------------------------------------+
|i2 : pts = {point{{0}},point{{1}},point{{2}}, point{{3}}, point{{4}}}|
|                                                                     |
|o2 = {{0}, {1}, {2}, {3}, {4}}                                       |
|                                                                     |
|o2 : List                                                            |
+---------------------------------------------------------------------+
|i3 : vals = {-1, 1/2, 1, 5/4, 7/5}                                   |
|                                                                     |
|          1     5  7                                                 |
|o3 = {-1, -, 1, -, -}                                                |
|          2     4  5                                                 |
|                                                                     |
|o3 : List                                                            |
+---------------------------------------------------------------------+
|i4 : rationalInterpolation(pts, vals, R)                             |
|-- warning: experimental computation over inexact field begun        |
|--          results not reliable (one warning given per session)     |
|                                                                     |
|o4 = (2x - 1, 1x + 1)                                                |
|                                                                     |
|o4 : Sequence                                                        |
+---------------------------------------------------------------------+

See also
========

  * rationalInterpolation(List,List,Matrix,Matrix)
    -- numerically interpolate rational functions

Ways to use this method:
========================

  * rationalInterpolation(List,List,Ring)
    -- numerically interpolate
    rational functions

Sampler -- optional sampler function
************************************

Description
===========

There are currently two methods that accept a user-specified sampler function

  * Strategy => "Hybrid", used with
    noetherianOperators,
  * numericalNoetherianOperators

See the respective documentation nodes for details and syntax. If the option
Sampler is unset, the default sampler will use Bertini.

See also
========

  * Strategy => "Hybrid" -- strategy for computing
    Noetherian operators
  * numericalNoetherianOperators, --
    Noetherian operators via numerical interpolation

For the programmer
==================

The object Sampler, is a symbol (Macaulay2Doc).

solvePDE(Module) -- solve linear systems of PDE with constant coefficients
**************************************************************************

Synopsis
========

  * Function: solvePDE,
  * Usage:
        solvePDE U
  * Inputs:
      * U, a module (Macaulay2Doc), a submodule of a free module,
        or a matrix (Macaulay2Doc), or an ideal (Macaulay2Doc),
  * Outputs:
      * a list (Macaulay2Doc),

Description
===========

A homogeneous system of $l$ linear partial differential equations for a
function $\psi \colon \RR^n \to \CC^k$ is encoded by a $(k \times l)$ matrix
with entries in a polynomial ring. For example, if $n=4, k=2, l=3$, the PDE
system $$ \frac{\partial^2 \psi_1}{\partial z_1 \partial z_3} +
\frac{\partial^2 \psi_2}{\partial z_1^2} = \frac{\partial^2 \psi_1}{\partial
z_1 \partial z_2} +   \frac{\partial^2 \psi_2}{\partial z_2^2} =
\frac{\partial^3 \psi_1}{\partial z_1^2 \partial z_2} + \frac{\partial^3
\psi_2}{\partial z_1^2 \partial z_4} =0 $$ is encoded by the matrix $$ M =
\begin{bmatrix} \partial_{1} \partial_{3} &  \partial_{1} \partial_{2} &
\partial_{1}^2 \partial_{2}\\ \partial_{1}^2 & \partial_{2}^2 & \partial_{1}^2
\partial_{4} \end{bmatrix}, $$ or more precisely, by the submodule $U$ of
$(K[\partial_1,\partial_2,\partial_3,\partial_4])^2$ generated by the columns
of $M$.

By the Ehrenpreis-Palamodov fundamental principle, solutions to such PDE are of
the form $$ \phi(\mathbf z) = \sum_{i = 1}^s \sum_{j=1}^{m_i} \int_{V_i}
B_{i,j}(\mathbf{x}, \mathbf{z}) e^{\mathbf{x}^t \mathbf{z}}
\mu_{i,j}(\mathbb{x}) $$ where $\mu_{i,j}$ are complex valued measures.

The function solvePDE computes the algebraic varieties $V_i$ and Noetherian
multipliers $B_{i,j}(\mathbf x, \mathbf z)$. The input is either a matrix $M$
or a module $U$, where the $\partial_i$ is replaced by $x_i$. The output is a
list of $s$ pairs. For the $i$th pair, the first entry is the prime ideal of
$V_i$. The second entry is the list $B_{i,1}, \dotsc, B_{i,m_j}$ of vectors of
polynomials in $2n$ variables, where the symbol $\mathbf{z}$ is replaced by the
symbol $\mathbf{\mathtt{d}x}$.

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x_1..x_4]                                                        |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : M = matrix{{x_1*x_3, x_1*x_2, x_1^2*x_2}, {x_1^2, x_2^2, x_1^2*x_4}}    |
|                                                                             |
|o2 = | x_1x_3 x_1x_2 x_1^2x_2 |                                              |
|     | x_1^2  x_2^2  x_1^2x_4 |                                              |
|                                                                             |
|             2      3                                                        |
|o2 : Matrix R  <-- R                                                         |
+-----------------------------------------------------------------------------+
|i3 : U = image M                                                             |
|                                                                             |
|o3 = image | x_1x_3 x_1x_2 x_1^2x_2 |                                        |
|           | x_1^2  x_2^2  x_1^2x_4 |                                        |
|                                                                             |
|                             2                                               |
|o3 : R-module, submodule of R                                                |
+-----------------------------------------------------------------------------+
|i4 : sols = solvePDE M                                                       |
|                                                                             |
|                                                                             |
|o4 = {{ideal x , {| 1 |}}, {ideal (x , x ), {| -x_1 |}}, {ideal (x , x ), {| |
|              1   | 0 |             4   2    |  x_3 |             3   2    | |
|     ------------------------------------------------------------------------|
|                                                                             |
|     1 |}}, {ideal (x , x ), {| x_2dx_1 |}}, {ideal (x , x ), {| 0 |, |   0  |
|     0 |             3   1    |    -1   |             2   1    | 1 |  | dx_1 |
|     ------------------------------------------------------------------------|
|                                           2                       2         |
|     |, |   0  |, |     0    |}}, {ideal (x  - x x , x x  - x x , x  - x x ),|
|     |  | dx_2 |  | dx_1dx_2 |             2    1 4   1 2    3 4   1    2 3  |
|     ------------------------------------------------------------------------|
|     {| -x_4 |}}}                                                            |
|      |  x_2 |                                                               |
|                                                                             |
|o4 : List                                                                    |
+-----------------------------------------------------------------------------+
|i5 : netList sols                                                            |
|                                                                             |
|     +-----------------------------------------+-----------------------------|
|o5 = |ideal x                                  |{| 1 |}                      |
|     |       1                                 | | 0 |                       |
|     +-----------------------------------------+-----------------------------|
|     |ideal (x , x )                           |{| -x_1 |}                   |
|     |        4   2                            | |  x_3 |                    |
|     +-----------------------------------------+-----------------------------|
|     |ideal (x , x )                           |{| 1 |}                      |
|     |        3   2                            | | 0 |                       |
|     +-----------------------------------------+-----------------------------|
|     |ideal (x , x )                           |{| x_2dx_1 |}                |
|     |        3   1                            | |    -1   |                 |
|     +-----------------------------------------+-----------------------------|
|     |ideal (x , x )                           |{| 0 |, |   0  |, |   0  |, ||
|     |        2   1                            | | 1 |  | dx_1 |  | dx_2 |  ||
|     +-----------------------------------------+-----------------------------|
|     |        2                       2        |                             |
|     |ideal (x  - x x , x x  - x x , x  - x x )|{| -x_4 |}                   |
|     |        2    1 4   1 2    3 4   1    2 3 | |  x_2 |                    |
|     +-----------------------------------------+-----------------------------|
|-----------------------------------------------------------------------------|
|------------+                                                                |
|            |                                                                |
|            |                                                                |
|------------+                                                                |
|            |                                                                |
|            |                                                                |
|------------+                                                                |
|            |                                                                |
|            |                                                                |
|------------+                                                                |
|            |                                                                |
|            |                                                                |
|------------+                                                                |
|     0    |}|                                                                |
| dx_1dx_2 | |                                                                |
|------------+                                                                |
|            |                                                                |
|            |                                                                |
|            |                                                                |
|------------+                                                                |
+-----------------------------------------------------------------------------+

This output reveals that the general solution to the example system above
consists of nine summands, one of which is $$ \phi(\mathbb z) =
\int_{V(x_3,x_1)} \begin{bmatrix} z_1 x_2 \\ -1 \end{bmatrix} e^{x_1 z_1 + x_2
z_2 + x_3 z_3 + x_4 z_4} \, d\mu(x_1,x_2,x_3,x_4) $$

The total number of Noetherian multipliers is equal to the arithmetic
multiplicity of the module $U$.

+--------------------------------------------+
|i6 : amult U == sum(sols / last / (l -> #l))|
|                                            |
|o6 = true                                   |
+--------------------------------------------+

Note that the output of solvePDE can be interpreted as a differential primary decomposition.

References
==========

Ait El Manssour, R., Härkönen, M., Sturmfels, B. (2021). Linear PDE with
constant coefficients (see https://arxiv.org/abs/2104.10146 )

See also
========

  * differentialPrimaryDecomposition,
    -- compute a differential primary decomposition
  * noetherianOperators(Module), --
    Noetherian operators of a primary submodule

Ways to use this method:
========================

  * "solvePDE(Ideal)"
  * "solvePDE(Matrix)"
  * solvePDE(Module) -- solve linear systems of
    PDE with constant coefficients

specializedNoetherianOperators -- Noetherian operators evaluated at a point
***************************************************************************

Synopsis
========

  * Usage:
        specializedNoetherianOperators(I, pt)
  * Inputs:
      * I, an ideal (Macaulay2Doc), unmixed
      * pt, a point (NAGtypes), or a one-row matrix (Macaulay2Doc),
  * Outputs:
      * a list (Macaulay2Doc), of differential operators,

Description
===========

Numerically computes evaluations of Noetherian operators. If the point p lies
on the variety of the minimal prime $P$, the function returns a set of
specialized Noetherian operators of the $P$-primary component of $I$. The
option DependentSet, is required when dealing with ideals
over numerical fields, or when dealing with non-primary ideals.

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x,y,t];                                                          |
+-----------------------------------------------------------------------------+
|i2 : Q1 = ideal(x^2, y^2 + x*t);                                             |
|                                                                             |
|o2 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i3 : Q2 = ideal((x+t)^2);                                                    |
|                                                                             |
|o3 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i4 : I = intersect(Q1, Q2);                                                  |
|                                                                             |
|o4 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i5 : P = radical Q1;                                                         |
|                                                                             |
|o5 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i6 : pt = point{{0,0,2}};                                                    |
+-----------------------------------------------------------------------------+
|i7 : A = specializedNoetherianOperators(I, pt, DependentSet => {x,y}) /      |
|                                                                             |
|o7 = {| 1 |, | dy |, | dy^2-dx |, | dy^3-3dxdy |}                            |
|                                                                             |
|o7 : List                                                                    |
|-----------------------------------------------------------------------------|
|normalize                                                                    |
+-----------------------------------------------------------------------------+
|i8 : B = noetherianOperators(I, P) /                                         |
|         (D -> evaluate(D,pt)) /                                             |
|         normalize                                                           |
|                                                                             |
|o8 = {| 1 |, | dy |, | dy^2-dx |, | dy^3-3dxdy |}                            |
|                                                                             |
|o8 : List                                                                    |
+-----------------------------------------------------------------------------+
|i9 : A == B                                                                  |
|                                                                             |
|o9 = true                                                                    |
+-----------------------------------------------------------------------------+

Over a non-exact field, the output will be non-exact

+--------------------------------------------------------------------------+
|i10 : S = CC[x,y,t]                                                       |
|                                                                          |
|o10 = S                                                                   |
|                                                                          |
|o10 : PolynomialRing                                                      |
+--------------------------------------------------------------------------+
|i11 : pt = point{{0,0,2.1}}                                               |
|                                                                          |
|o11 = pt                                                                  |
|                                                                          |
|o11 : Point                                                               |
+--------------------------------------------------------------------------+
|i12 : specializedNoetherianOperators(sub(I, S), pt, DependentSet => {x,y})|
|                                                                          |
|o12 = {| 1 |, | dy |, | .5dy^2-.47619dx |, | .166667dy^3-.47619dxdy |}    |
|                                                                          |
|o12 : List                                                                |
+--------------------------------------------------------------------------+

Caveat
======

It is assumed that the point lies on the variety of I

Ways to use specializedNoetherianOperators:
===========================================

  * "specializedNoetherianOperators(Ideal,AbstractPoint)"
  * "specializedNoetherianOperators(Ideal,Matrix)"

For the programmer
==================

The object specializedNoetherianOperators,
is a method function with options (Macaulay2Doc).

Strategy => "Hybrid" -- strategy for computing Noetherian operators
*******************************************************************

This strategy implements a numerical-symbolic hybrid algorithm for computing
Noetherian operators. The output is symbolic. "Hybrid" supports computing
Noetherian operators of either primary ideals (see
noetherianOperators(Ideal)), or primary
components of unmixed ideals (see noetherianOperators(Ideal,Ideal)).

The "Hybrid" strategy finds a point on the variety of the component of
interest, and computes a set of specialized Noetherian operators (see
specializedNoetherianOperators). Using this
numerical data is then used as a starting point for the symbolic computation of
Noetherian operators, which in many cases lead to significant performance
improvements over the fully symbolic methods.

The strategy accepts the following optional arguments:

Sampler => f, where f is a function taking a primary ideal and returning a
single point on the variety. The default sampler uses a combination of
bertiniSample (Bertini), and bertiniPosDimSolve (Bertini).
The user can supply a point to used by using a
dummy sampler, as in the example below:

+--------------------------------------------------------------------+
|i1 : R = QQ[x,y,t];                                                 |
+--------------------------------------------------------------------+
|i2 : I = ideal(x^2, y^2-x*t);                                       |
|                                                                    |
|o2 : Ideal of R                                                     |
+--------------------------------------------------------------------+
|i3 : p = point{{0_CC,0, 3}};                                        |
+--------------------------------------------------------------------+
|i4 : noetherianOperators(I, Strategy => "Hybrid", Sampler => I -> p)|
|-- warning: experimental computation over inexact field begun       |
|--          results not reliable (one warning given per session)    |
|                                                                    |
|o4 = {| 1 |, | dy |, | tdy^2+2dx |, | tdy^3+6dxdy |}                |
|                                                                    |
|o4 : List                                                           |
+--------------------------------------------------------------------+

Tolerance => ...: takes a positive real number. This specifies the numerical
precision when computing the specialized Noetherian operators. The default
value is 1e-6. See See Tolerance (NoetherianOperators).

DependentSet =>...: takes a list of variables. For details, see
DependentSet.

See also
========

  * Strategy => "PunctualQuot" -- strategy
    for computing Noetherian operators
  * Strategy => "MacaulayMatrix", --
    strategy for computing Noetherian operators

Strategy => "MacaulayMatrix" -- strategy for computing Noetherian operators
***************************************************************************

This strategy implements Algorithm 2 in the paper Noetherian Operators and
Primary Decomposition (see https://arxiv.org/abs/2006.13881 ), and supports
computing Noetherian operators of either primary ideals (see
noetherianOperators(Ideal)), or primary
components of unmixed ideals (see noetherianOperators(Ideal,Ideal)).

The strategy relies on computing the kernel of successively larger Macaulay
matrices. The behavior can be controlled with optional arguments:

DegreeLimit => ...: takes an integer $d$, and stops computation at degree $d$.
Note that if $d$ is set too low, this may lead to an incomplete answer. If
unset, stops computation when the dimension of the kernel stabilizes.

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x,y,z];                                                          |
+-----------------------------------------------------------------------------+
|i2 : I = (ideal(x,y,z))^3;                                                   |
|                                                                             |
|o2 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i3 : noetherianOperators(I, Strategy => "MacaulayMatrix")                    |
|                                                                             |
|o3 = {| 1 |, | dz |, | dy |, | dx |, | dz^2 |, | dydz |, | dxdz |, | dy^2 |, |
|     ------------------------------------------------------------------------|
|     | dxdy |, | dx^2 |}                                                     |
|                                                                             |
|o3 : List                                                                    |
+-----------------------------------------------------------------------------+
|i4 : noetherianOperators(I, Strategy => "MacaulayMatrix", DegreeLimit => 1)  |
|                                                                             |
|o4 = {| 1 |, | dz |, | dy |, | dx |}                                         |
|                                                                             |
|o4 : List                                                                    |
+-----------------------------------------------------------------------------+

KernelStrategy => ...: takes a string "Default" or "Gaussian". The "Default"
strategy uses the Macaulay2 builtin function kernel (Macaulay2Doc),
to compute kernels (via Grobner bases). The strategy
"Gaussian" computes kernels directly via a Gaussian reduction, and may offer
performance improvements compared to "Default".

IntegralStrategy => ...: takes a boolean value. If true, uses the Mourrain
algorithm to compute the kernel of the MacaulayMatrix, which constructs columns
of the Macaulay matrix by taking integrals of the columns in the previous step.
If false, uses the method outlined in Algorithm 1 in the paper Noetherian
Operators and Primary Decomposition (see https://arxiv.org/abs/2006.13881 ). If
unset, will choose automatically. See: B. Mourrain. Isolated points, duality
and residues. J. Pure Appl. Algebra, 117/118:469-493, 1997. Algorithms for
algebra (Eindhoven, 1996).

DependentSet => ...: takes a list of variables. For details, see
DependentSet.

See also
========

  * Strategy => "Hybrid" -- strategy for computing
    Noetherian operators
  * Strategy => "PunctualQuot" -- strategy
    for computing Noetherian operators
  * DependentSet -- option for computing Noetherian
    operators

Strategy => "PunctualQuot" -- strategy for computing Noetherian operators
*************************************************************************

This strategy implements Algorithm 4.1 in the paper Primary decomposition of
modules: a computational differential approach (see
https://arxiv.org/abs/2104.03385 ).

The following example deals with a rather non-trivial primary ideal to show the
capabilities of this strategy.

+-----------------------------------------------------------------------------+
|i1 : R = QQ[x_1,x_2,x_3,x_4]                                                 |
|                                                                             |
|o1 = R                                                                       |
|                                                                             |
|o1 : PolynomialRing                                                          |
+-----------------------------------------------------------------------------+
|i2 : k = 3                                                                   |
|                                                                             |
|o2 = 3                                                                       |
+-----------------------------------------------------------------------------+
|i3 : J = ideal((x_1^2-x_2*x_3)^k,(x_1*x_2-x_3*x_4)^k,(x_2^2-x_1*x_4)^k)      |
|                                                                             |
|             6     4         2 2 2    3 3   3 3     2 2             2 2      |
|o3 = ideal (x  - 3x x x  + 3x x x  - x x , x x  - 3x x x x  + 3x x x x  -    |
|             1     1 2 3     1 2 3    2 3   1 2     1 2 3 4     1 2 3 4      |
|     ------------------------------------------------------------------------|
|      3 3   6       4       2 2 2    3 3                                     |
|     x x , x  - 3x x x  + 3x x x  - x x )                                    |
|      3 4   2     1 2 4     1 2 4    1 4                                     |
|                                                                             |
|o3 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i4 : Q = saturate(J,ideal(x_1*x_2*x_3*x_4))                                  |
|                                                                             |
|             5      3 2         3        4 2     2     2    2 2 2      2 3   |
|o4 = ideal (x x  + x x x  - 4x x x x  - x x  + 3x x x x  + x x x  - x x x ,  |
|             2 3    1 2 4     1 2 3 4    1 4     1 2 3 4    2 3 4    1 3 4   |
|     ------------------------------------------------------------------------|
|        4      4         2 2        3 2      3   2         2 2    3 3   2 3  |
|     x x x  + x x x  - 3x x x x  - x x x  - x x x  + 4x x x x  - x x , x x x |
|      1 2 3    1 2 4     1 2 3 4    2 3 4    1 3 4     1 2 3 4    3 4   1 2 3|
|     ------------------------------------------------------------------------|
|        4 2    5       3             2 2      2 2 2      3 2   6       4     |
|     - x x  + x x  - 4x x x x  + 3x x x x  + x x x  - x x x , x  - 3x x x  + |
|        2 3    1 4     1 2 3 4     1 2 3 4    1 3 4    2 3 4   2     1 2 4   |
|     ------------------------------------------------------------------------|
|       2 2 2    3 3   3 3     2 2             2 2    3 3   6     4           |
|     3x x x  - x x , x x  - 3x x x x  + 3x x x x  - x x , x  - 3x x x  +     |
|       1 2 4    1 4   1 2     1 2 3 4     1 2 3 4    3 4   1     1 2 3       |
|     ------------------------------------------------------------------------|
|       2 2 2    3 3                                                          |
|     3x x x  - x x )                                                         |
|       1 2 3    2 3                                                          |
|                                                                             |
|o4 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i5 : isPrimary Q                                                             |
|                                                                             |
|o5 = true                                                                    |
+-----------------------------------------------------------------------------+
|i6 : elapsedTime noetherianOperators(Q, Strategy => "PunctualQuot")          |
| -- .0499782s elapsed                                                        |
|                                                                             |
|o6 = {| 1 |, | dx_1 |, | dx_2 |, | dx_1^2 |, | dx_1dx_2 |, | dx_2^2 |, |     |
|     ------------------------------------------------------------------------|
|     2x_1x_3dx_1^3+3x_2x_3dx_1^2dx_2-3x_3x_4dx_1dx_2^2-2x_1x_4dx_2^3 |}      |
|                                                                             |
|o6 : List                                                                    |
+-----------------------------------------------------------------------------+

See also
========

  * mapToPunctualHilbertScheme -- maps an
    ideal into a point in a certain punctual Hilbert scheme
  * Strategy => "MacaulayMatrix", --
    strategy for computing Noetherian operators
  * Strategy => "Hybrid" -- strategy for computing
    Noetherian operators

Tolerance (NoetherianOperators) -- optional argument for numerical tolerance
****************************************************************************

Many of the numerical operations require a tolerance value, below which numbers
are considered to be numerically zero. Tolerance should be given a non-negative
real number. Over exact fields the default value is zero, while for inexact
fields the default is 1e-6.

See also Tolerance (NumericalLinearAlgebra).

truncate(DualSpace,List,ZZ) -- truncate a polynomial space or dual space
************************************************************************

Synopsis
========

  * Function: truncate (Macaulay2Doc),
  * Usage:
        S = truncate(T, v, d)
  * Inputs:
      * T, an instance of the type DualSpace (NAGtypes), or an
        instance of the type PolySpace (NAGtypes),
      * v, a list (Macaulay2Doc), a list of the integers
        designating which variables to bound
      * d, an integer (Macaulay2Doc), the degree bound for the
        designated variables
  * Outputs:
      * S, an instance of the type DualSpace (NAGtypes),

Description
===========

This method truncates a dual space or polynomial space T, so that the total
degree of the specified variables is bounded by d.

+-------------------------------------------+
|i1 : R = CC[x,y];                          |
+-------------------------------------------+
|i2 : I = ideal {x,y};                      |
|                                           |
|o2 : Ideal of R                            |
+-------------------------------------------+
|i3 : D = zeroDimensionalDual(origin R, I^3)|
|                                           |
|o3 = | 1 y x .5y2 xy .5x2 |                |
|                                           |
|o3 : DualSpace                             |
+-------------------------------------------+
|i4 : truncate(D, {0}, 1)                   |
|                                           |
|o4 = | xy x y2 y 1 |                       |
|                                           |
|o4 : DualSpace                             |
+-------------------------------------------+

Ways to use this method:
========================

  * truncate(DualSpace,List,ZZ), --
    truncate a polynomial space or dual space
  * "truncate(PolySpace,List,ZZ)"

truncate(DualSpace,ZZ) -- truncate a polynomial space or dual space
*******************************************************************

Synopsis
========

  * Function: truncate (Macaulay2Doc),
  * Usage:
        S = truncate(T, d)
  * Inputs:
      * T, an instance of the type DualSpace (NAGtypes), or an
        instance of the type PolySpace (NAGtypes),
      * d, an integer (Macaulay2Doc), the degree bound for the
        designated variables
  * Outputs:
      * S, an instance of the type DualSpace (NAGtypes),

Description
===========

This method truncates a dual space or polynomial space T so that the total
degree of all terms does not exceed d.

+-----------------------------------------------------------------------------+
|i1 : R = CC[x,y];                                                            |
+-----------------------------------------------------------------------------+
|i2 : I = ideal {x-y};                                                        |
|                                                                             |
|o2 : Ideal of R                                                              |
+-----------------------------------------------------------------------------+
|i3 : D = truncatedDual(origin R, I, 5)                                       |
|                                                                             |
|o3 = | 1 x+y .5x2+xy+.5y2 .166667x3+.5x2y+.5xy2+.166667y3                    |
|                                                                             |
|o3 : DualSpace                                                               |
|-----------------------------------------------------------------------------|
|.0416667x4+.166667x3y+.25x2y2+.166667xy3+.0416667y4                          |
|-----------------------------------------------------------------------------|
|.00833333x5+.0416667x4y+.0833333x3y2+.0833333x2y3+.0416667xy4+.00833333y5 |  |
+-----------------------------------------------------------------------------+
|i4 : truncate(D, 3)                                                          |
|                                                                             |
|o4 = | x3+3x2y+3xy2+y3 x2+2xy+y2 x+y 1 |                                     |
|                                                                             |
|o4 : DualSpace                                                               |
+-----------------------------------------------------------------------------+

Ways to use this method:
========================

  * truncate(DualSpace,ZZ) -- truncate a
    polynomial space or dual space
  * "truncate(PolySpace,ZZ)"

truncatedDual -- truncated dual space of a polynomial ideal
***********************************************************

Synopsis
========

  * Usage:
        S = truncatedDual(p, I, d)
  * Inputs:
      * p, a point (NAGtypes), or a one-row matrix (Macaulay2Doc),
      * I, an ideal (Macaulay2Doc), or a one-row matrix (Macaulay2Doc),
        of generators
      * d, an integer (Macaulay2Doc),
  * Optional inputs (Macaulay2Doc),:
      * Tolerance, => ..., default
        value null, optional argument for numerical tolerance
  * Outputs:
      * S, an instance of the type DualSpace (NAGtypes),

Description
===========

This functions computes a basis for the local dual space of a polynomial ideal
localized at point p, truncated at degree d. Elements are expressed as elements
of the polynomial ring of the ideal although this is an abuse of notation. They
are really elements of the dual ring.

+---------------------------------+
|i1 : R = CC[x,y];                |
+---------------------------------+
|i2 : I = ideal{x^2, y*x}         |
|                                 |
|             2                   |
|o2 = ideal (x , x*y)             |
|                                 |
|o2 : Ideal of R                  |
+---------------------------------+
|i3 : truncatedDual(origin(R),I,3)|
|                                 |
|o3 = | 1 y x .5y2 .166667y3 |    |
|                                 |
|o3 : DualSpace                   |
+---------------------------------+

The functionals in the dual at a point p are expressed in coordinates centered
at p.

+-----------------------------------+
|i4 : p = point matrix{{0_CC, 1_CC}}|
|                                   |
|o4 = p                             |
|                                   |
|o4 : Point                         |
+-----------------------------------+
|i5 : truncatedDual(p,I,3)          |
|                                   |
|o5 = | 1 y .5y2 .166667y3 |        |
|                                   |
|o5 : DualSpace                     |
+-----------------------------------+

Over inexact fields, the computation accounts for the possibility of small
numerical error in the point p. The optional argument Tolerance (NoetherianOperators),
can be set to
adjust the tolerance of the numerical computations. Higher degree dual
computations generally require higher accuracy in the input and larger
tolerance value to complete correctly.

In this example, the point q is slightly away from the variety of I, but an
appropriate Tolerance value can overcome the error.

+-----------------------------------------------+
|i6 : q = point matrix{{0_CC + 1e-10, 1_CC}}    |
|                                               |
|o6 = q                                         |
|                                               |
|o6 : Point                                     |
+-----------------------------------------------+
|i7 : tol = 1e-6;                               |
+-----------------------------------------------+
|i8 : S = truncatedDual(q,I,3, Tolerance => tol)|
|                                               |
|o8 = | 1 y .5y2 .166667y3 |                    |
|                                               |
|o8 : DualSpace                                 |
+-----------------------------------------------+
|i9 : (m,c) = coefficients gens S;              |
+-----------------------------------------------+
|i10 : m*clean(tol, c)                          |
|                                               |
|o10 = | 1 y .5y2 .166667y3 |                   |
|                                               |
|              1      4                         |
|o10 : Matrix R  <-- R                          |
+-----------------------------------------------+

See also
========

  * zeroDimensionalDual -- dual space of a
    zero-dimensional polynomial ideal

Ways to use truncatedDual:
==========================

  * "truncatedDual(AbstractPoint,Ideal,ZZ)"
  * "truncatedDual(AbstractPoint,Matrix,ZZ)"
  * "truncatedDual(Matrix,Ideal,ZZ)"
  * "truncatedDual(Matrix,Matrix,ZZ)"

For the programmer
==================

The object truncatedDual, is a method function with options (Macaulay2Doc).

zeroDimensionalDual -- dual space of a zero-dimensional polynomial ideal
************************************************************************

Synopsis
========

  * Usage:
        S = zeroDimensionalDual(p, I)
  * Inputs:
      * p, a point (NAGtypes), or a one-row matrix (Macaulay2Doc),
      * I, an ideal (Macaulay2Doc), or a one-row matrix (Macaulay2Doc),
        of generators
  * Optional inputs (Macaulay2Doc),:
      * Tolerance, => ..., default
        value null, optional argument for numerical tolerance
  * Outputs:
      * S, an instance of the type DualSpace (NAGtypes),

Description
===========

This function computes a reduced basis of the dual space of a zero-dimensional
ideal.  It does not check if the ideal is zero-dimensional and if not then
termination will fail. Elements are expressed as elements of the polynomial
ring of the ideal although this is an abuse of notation. They are really
elements of the dual ring.

+---------------------------------------------------+
|i1 : R = QQ[a,b];                                  |
+---------------------------------------------------+
|i2 : I = ideal{a^3,b^3}                            |
|                                                   |
|             3   3                                 |
|o2 = ideal (a , b )                                |
|                                                   |
|o2 : Ideal of R                                    |
+---------------------------------------------------+
|i3 : D = zeroDimensionalDual(origin(R), I)         |
|                                                   |
|o3 = | 1 b a 1/2b2 ab 1/2a2 1/2ab2 1/2a2b 1/4a2b2 ||
|                                                   |
|o3 : DualSpace                                     |
+---------------------------------------------------+
|i4 : dim D                                         |
|                                                   |
|o4 = 9                                             |
+---------------------------------------------------+

The dimension of the dual space at p is the multiplicity of the solution at p.

+-------------------------------------------+
|i5 : S = CC[x,y];                          |
+-------------------------------------------+
|i6 : J = ideal{(y-2)^2,y-x^2}              |
|                                           |
|             2              2              |
|o6 = ideal (y  - 4y + 4, - x  + y)         |
|                                           |
|o6 : Ideal of S                            |
+-------------------------------------------+
|i7 : p = point matrix{{1.4142136_CC,2_CC}};|
+-------------------------------------------+
|i8 : D = zeroDimensionalDual(p, J)         |
|                                           |
|o8 = | 1 .353553x+y |                      |
|                                           |
|o8 : DualSpace                             |
+-------------------------------------------+
|i9 : dim D                                 |
|                                           |
|o9 = 2                                     |
+-------------------------------------------+

Caveat
======

The computation will not terminate if I is not locally zero-dimensional at the
chosen point.  This is not checked.

See also
========

  * truncatedDual -- truncated dual space of a polynomial
    ideal

Ways to use zeroDimensionalDual:
================================

  * "zeroDimensionalDual(AbstractPoint,Ideal)"
  * "zeroDimensionalDual(AbstractPoint,Matrix)"
  * "zeroDimensionalDual(Matrix,Ideal)"
  * "zeroDimensionalDual(Matrix,Matrix)"

For the programmer
==================

The object zeroDimensionalDual, is a method function with options (Macaulay2Doc).