packages/lpdo/lpdo.tex

% $Id: lpdo.tex 550 2010-02-24 19:31:24Z thomas-sturm $
\documentclass[a4paper]{article}

\usepackage{amsmath,amssymb}
\usepackage{enumerate}

\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Sym}{\operatorname{Sym}}
\newenvironment{reduce}{\begin{scriptsize}\begin{verbatim}}%
{\end{verbatim}\end{scriptsize}}

\newcommand{\txt}{\textrm}
\newenvironment{ebnf}{\begin{tabular}{rcll}}{\end{tabular}}
\newcommand{\ebnfvar}[1]{{\textit{#1}}\ \ }
\newcommand{\ebnfvara}[1]{$\left\langle \text{\textit{#1}}\right\rangle$\ \ }
\newcommand{\ebnfor}{\txt{\ $|$\ \ }}
\newcommand{\ebnfterm}[1]{\txt{\textbf{\texttt{#1}}}\ \ }
\newcommand{\ebnfoption}[1]{\txt{$\left[\;\text{#1}\;\right]$}\ \ }
\newcommand{\ebnfmult}[1]{\txt{$\left\{\;\text{#1}\;\right\}$}\ \ }
\newcommand{\ebnfattrib}[1]{{$\bigl\{$\ {#1}\ $\bigr\}$}\ \ }

\newcommand{\ebnfrule}[2]{{#1}&$\ \ \rightarrow\ \ $&{#2}&\\}
\newcommand{\ebnfrulev}[2]{{\ebnfvar{#1}}&$\ \ \rightarrow\ \ $&{#2}&\\}
\newcommand{\ebnfhalfrule}[1]{&&{\ebnfor {#1}}&\\}
\newcommand{\ebnfrulea}[3]{{#1}&$\ \ \rightarrow\ \ $&{#2}&{$\bigl\{$\ {#3}\ $\bigr\}$}\\}
\newcommand{\ebnfruleav}[3]{{\ebnfvar{#1}}&$\ \ \rightarrow\ \ $&{#2}&{$\bigl\{$\ {#3}\ $\bigr\}$}\\}
\newcommand{\ebnfhalfrulea}[2]{&&{\ebnfor {#1}}&{$\bigl\{$\ {#2}\ $\bigr\}$}\\}

\newcommand{\ebnfsinglerulev}[2]{\txt{{\ebnfvar{#1}}$\ \ \rightarrow\ \ ${#2}\\}}
\newcommand{\ebnfsinglerule}[2]{\txt{{#1}$\ \ \rightarrow\ \ ${#2}\\}}
\newcommand{\ebnfsingleruleav}[3]{\txt{{\ebnfvar{#1}}$\ \ \rightarrow\ \ ${#2}\ \ \ $\bigl\{$\ {#3}\ $\bigr\}$\\}}
\newcommand{\ebnfsinglerulea}[3]{\txt{{#1}$\ \ \rightarrow\ \ ${#2}\ \
    \ $\bigl\{$\ {#3}\ $\bigr\}$\\}}

\begin{document}
\title{LPDO: Linear Partial Differential Operators}
\author{Thomas Sturm\\
  Departamento de Matem\'aticas,
  Estad\'istica y Computaci\'on,\\
  Facultad de Ciencias,
  Universidad de Cantabria\\
  39071 Santander, Spain\\
  \texttt{sturmt@unican.es}}
\maketitle

\begin{abstract}

\end{abstract}

\section{Introduction}
Consider the field $F = \Q(x_1,\dots,x_n)$ of rational functions and a
set $\Delta = \{\partial_{x_1}$,\dots,~$\partial_{x_n}\}$ of
\emph{commuting derivations} acting on $F$. That is, for all
$\partial_{x_i}$, $\partial_{x_j} \in \Delta$ and all $f$, $g\in F$
the following properties are satisfied:
\begin{eqnarray}
\partial_{x_i}(f + g) &=& \partial_{x_i}(f) + \partial_{x_i}(g),\nonumber\\
\partial_{x_i}(f\cdot g) &=& f\cdot\partial_{x_i}(g)
  + \partial_{x_i}(f)\cdot g,\label{EQ:leibnitz}\\
\partial_{x_i}(\partial_{x_j}(f)) &=& \partial_{x_j}(\partial_{x_i}(f)).\label{EQ:comm}
\end{eqnarray}
%We call $(F,\Delta)$ a \emph{$\Delta$-field}.

Consider now the set $F[\partial_{x_1},\dots,\partial_{x_n}]$, where
the derivations are used as variables. This set forms a
non-commutative \emph{linear partial differential operator ring} with
pointwise addition, and multiplication defined as follows: For $f\in
F$ and $\partial_{x_i}$, $\partial_{x_j} \in \Delta$ we have for any
$g\in F$ that
\begin{eqnarray}
    (f \partial_{x_i})(g) &=& f \cdot \partial_{x_i}(g),\nonumber\\
    (\partial_{x_i} f)(g) &=& \partial_{x_i}(f \cdot g),\label{EQ:mult2}\\
    (\partial_{x_i} \partial_{x_j})(g) &=& \partial_{x_i}(\partial_{x_j}(g)).\label{EQ:mult3}
\end{eqnarray}
Here ``${}\cdot{}$'' denotes the multiplication in $F$. From
(\ref{EQ:mult3}) and (\ref{EQ:comm}) it follows that
$\partial_{x_i} \partial_{x_j}=\partial_{x_j} \partial_{x_i}$, and
using (\ref{EQ:mult2}) and (\ref{EQ:leibnitz}) the following
\emph{commutator} can be proved:
\begin{eqnarray*}
  \partial_{x_i} f &=& f \partial_{x_i} + \partial_{x_i}(f).
\end{eqnarray*}

A \emph{linear partial differential operator} (LPDO) of order $k$ is
an element
\begin{eqnarray*}
  D &=& \sum_{|j| \leq k} a_j \partial^j\in
  F[\partial_{x_1},\dots,\partial_{x_n}]
\end{eqnarray*}
in canonical form. Here the expression $|j| \leq k$ specifies the set
of all tuples of the form $j = (j_1,\dots,j_n) \in \N^n$ with
$\sum_{i=1}^n j_i \leq k$, and we define $\partial^j
= \partial_{x_1}^{j_1} \cdots \partial_{x_n}^{j_n}$.

A \emph{factorization} of $D$ is a non-trivial decomposition
\begin{eqnarray*}
D &=& D_1 \cdots D_r\in F[\partial_{x_1},\dots,\partial_{x_n}]
\end{eqnarray*}
into multiplicative factors, each of which is an LPDO $D_i$ of order
greater than $0$ and less than $k$. If such a factorization exists,
then $D$ is called \emph{reducible} or \emph{factorable}, else
\emph{irreducible}.

For the purpose of factorization it is helpful to temporarily consider
as regular commutative polynomials certain summands of the LPDO under
consideration. Consider a commutative polynomial ring over $F$ in new
indeterminates $y_1$, \dots,~$y_n$. Adopting the notational
conventions above, for $m\leq k$ the \emph{symbol of $D$ of order $m$}
is defined as
\begin{eqnarray*}
  \Sym_m(D) &=& \sum_{|j| = m} a_j y^j\in F[y_1,\dots,y_n].
\end{eqnarray*}
For $m=k$ we obtain as a special case the \emph{symbol $\Sym(D)$ of
  $D$}.

\section{Operators}

\subsection{\texttt{partial}}
There is a unary operator \texttt{partial($\cdot$)} denoting
$\partial$.

\begin{center}
  \begin{ebnf}
    \ebnfrule{\ebnfvara{partial-term}}{\ebnfterm{partial}\ebnfterm{(}\ebnfvara{id}\ebnfterm{)}}
  \end{ebnf}
\end{center}

\subsection{\texttt{***}}
There is a binary operator \texttt{***} for the non-commutative
multiplication involving partials $\partial_x$. All expressions
involving \texttt{***} are implicitly transformed into LPDOs, i.e.,
into the following normal form:
\begin{center}
  \begin{ebnf}
    \ebnfrule{\ebnfvara{normalized-lpdo}}{\ebnfvara{normalized-mon}\ebnfoption{\ebnfterm{+}\ebnfvara{normalized-lpdo}}}
    \ebnfrule{\ebnfvara{normalized-mon}}{\ebnfvara{F-element}\ebnfoption{\ebnfterm{***}\ebnfvara{partial-termprod}}}
%    \ebnfrule{\ebnfvara{F-element}}{\ebnfvara{standard-quotient}\ebnfor\ebnfvara{constant-operator}}
    \ebnfrule{\ebnfvara{partial-termprod}}{\ebnfvara{partial-term}\ebnfoption{\ebnfterm{***}\ebnfvara{partial-termprod}}}
  \end{ebnf}
\end{center}
The summands of the \textit{normalized-lpdo} are ordered in some
canonical way. As an example consider

\begin{footnotesize}
\begin{verbatim}
input: a()***partial(y)***b()***partial(x);

(a()*b()) *** partial(x) *** partial(y) + (a()*diff(b(),y,1)) *** partial(x)
\end{verbatim}
\end{footnotesize}
Here the \textit{F-elements} are polynomials, where the unknowns are
of the type \emph{constant-operator} denoting functions from $F$:

\begin{center}
  \begin{ebnf}
    \ebnfrule{\ebnfvara{constant-operator}}{\ebnfvara{id}\ebnfterm{(}\ebnfterm{)}}
  \end{ebnf}
\end{center}
We do not admit division of such constant operators since we cannot
exclude that such a constant operator denotes $0$.

The operator notation on the one hand emphasizes the fact that the
denoted elements are functions. On the other hand it distinguishes
\texttt{a()} from the variable \texttt{a} of a rational function,
which specifically denotes the corresponding projection. Consider e.g.

\begin{footnotesize}
\begin{verbatim}
input: (x+y)***partial(y)***(x-y)***partial(x);

  2    2
(x  - y ) *** partial(x) *** partial(y) + ( - x - y) *** partial(x)
\end{verbatim}
\end{footnotesize}
Here we use as \textit{F-elements} specific elements from $F=\Q(x,y)$.

\subsection{\texttt{diff}}
In our example with constant operators, the transformation into normal
form introduces a formal derivative operation
\texttt{diff($\cdot$,$\cdot$,$\cdot$)}, which cannot be evaluated.
Notice that we do not use the Reduce operator
\texttt{df($\cdot$,$\cdot$,$\cdot$)} here, which for technical reasons
cannot smoothly handle our constant operators.

In our second example with rational functions as \textit{F-elements},
derivative occurring with commutation can be computed such that
\texttt{diff} does not occur in the output.

\section{Shapes of F-elements}
Besides the generic computations with constant operators, we provide a
mechanism to globally fix a certain \emph{shape} for
\textit{F-elements} and to expand constant operators according to that
shape.
\subsection{\texttt{lpdoset}}
We give an example for a shape that fixes all constant operators to
denote generic bivariate affine linear functions:

\begin{footnotesize}
\begin{verbatim}
input: d := (a()+b())***partial(x1)***partial(x2)**2;

                                                2
d := (a() + b()) *** partial(x1) *** partial(x2)

input: lpdoset {!#10*x1+!#01*x2+!#00,x1,x2};

{-1}

input: d;

                                                                               2
(a00 + a01*x2 + a10*x1 + b00 + b01*x2 + b10*x1) *** partial(x1) *** partial(x2)
\end{verbatim}
\end{footnotesize}
Notice that the placeholder \texttt{\#} must be escaped with
\texttt{!}, which is a general convention for Rlisp/Reduce. Notice
that \texttt{lpdoset} returns the old shape and that \texttt{\{-1\}}
denotes the default state that there is no shape selected.

\subsection{\texttt{lpdoweyl}}
The command \texttt{lpdoweyl \{n,x1,x2,...\}} creates a shape for
generic polynomials of total degree \texttt{n} in variables
$\texttt{x1}$, $\texttt{x2}$, \dots.

\begin{footnotesize}
\begin{verbatim}
input: lpdoweyl(2,x1,x2);

                            2                                    2
{#_00_ + #_01_*x2 + #_02_*x2  + #_10_*x1 + #_11_*x1*x2 + #_20_*x1 ,x1,x2}

input: lpdoset ws;

{#10*x1 + #01*x2 + #00,x1,x2}

input: d;

                            2                                    2
(a_00_ + a_01_*x2 + a_02_*x2  + a_10_*x1 + a_11_*x1*x2 + a_20_*x1  + b_00_

                       2                                    2
  + b_01_*x2 + b_02_*x2  + b_10_*x1 + b_11_*x1*x2 + b_20_*x1 ) *** partial(x1)

                2
 *** partial(x2)
\end{verbatim}
\end{footnotesize}

\section{Commands}

\subsection{General}
\subsubsection{\texttt{lpdoord}}
The \emph{order} of an lpdo:
\begin{footnotesize}
\begin{verbatim}
input: lpdoord((a()+b())***partial(x1)***partial(x2)**2+3***partial(x1));

3
\end{verbatim}
\end{footnotesize}
%
\subsubsection{\texttt{lpdoptl}}
Returns the list of derivations (partials) occurring in its argument
LPDO $d$.

\begin{footnotesize}
\begin{verbatim}
input: lpdoptl(a()***partial(x1)***partial(x2)+partial(x4)+diff(a(),x3,1));

{partial(x1),partial(x2),partial(x4)}
\end{verbatim}
\end{footnotesize}
That is the smallest set $\{\dots,\partial_{x_i},\dots\}$ such that
$d$ is defined in $F[\dots,\partial_{x_i},\dots]$. Notice that formal
derivatives are not derivations in that sense.

\subsubsection{\texttt{lpdogp}}
Given a starting symbol $a$, a list of variables $l$, and a degree
$n$, \texttt{lpdogp($a$,$l$,$n$)} generates a generic (commutative)
polynomial of degree $n$ in variables $l$ with coefficients generated
from the starting symbol $a$:

\begin{footnotesize}
\begin{verbatim}
input: lpdogp(a,{x1,x2},2);

                           2                                    2
a_00_ + a_01_*x2 + a_02_*x2  + a_10_*x1 + a_11_*x1*x2 + a_20_*x1
\end{verbatim}
\end{footnotesize}

\subsubsection{\texttt{lpdogdp}}
Given a starting symbol $a$, a list of variables $l$, and a degree
$n$, \texttt{lpdogp($a$,$l$,$n$)} generates a generic differential
polynomial of degree $n$ in variables $l$ with coefficients generated
from the starting symbol $a$:

\begin{footnotesize}
\begin{verbatim}
input: lpdogdp(a,{x1,x2},2);

                     2                        2
a_20_ *** partial(x1)  + a_02_ *** partial(x2)

 + a_11_ *** partial(x1) *** partial(x2) + a_10_ *** partial(x1)

 + a_01_ *** partial(x2) + a_00_
\end{verbatim}
\end{footnotesize}

\subsection{Symbols}
\subsubsection{\texttt{lpdosym}}
The \emph{symbol} of an lpdo. That is the differential monomial of
highest order with the partials replaced by corresponding commutative
variables:

\begin{footnotesize}
\begin{verbatim}
input: lpdosym((a()+b())***partial(x1)***partial(x2)**2+3***partial(x1));

           2
y_x1_*y_x2_ *(a() + b())
\end{verbatim}
\end{footnotesize}
More generally, one can use a second optional arguments to specify a
the order of a different differential monomial to form the symbol of:

\begin{footnotesize}
\begin{verbatim}
input: lpdosym((a()+b())***partial(x1)***partial(x2)**2+3***partial(x1),1);

3*y_x1_
\end{verbatim}
\end{footnotesize}
Finally, a third optional argument can be used to specify an
alternative starting symbol for the commutative variable, which is
\texttt{y} by default. Altogether, the optional arguments default like
\texttt{lpdosym($\cdot$)}=\texttt{lpdosym($\cdot$,lpdoord($\cdot$),y)}.

\subsubsection{\texttt{lpdosym2dp}}
This converts a symbol obtained via \texttt{lpdosym} back into an LPDO
resulting in the corresponding differential monomial of the original
LPDO.

\begin{footnotesize}
\begin{verbatim}
input: d := a()***partial(x1)***partial(x2)+partial(x3)$

input: s := lpdosym d;

s := a()*y_x1_*y_x2_

input: lpdosym2dp s;

a() *** partial(x1) *** partial(x2)
\end{verbatim}
\end{footnotesize}
In analogy to \texttt{lpdosym} there is an optional argument for
specifying an alternative starting symbol for the commutative
variable, which is \texttt{y} by default.

\subsubsection{\texttt{lpdos}}
Given LPDOs $p$, $q$ and $m\in\N$ the function
\texttt{lpdos($p$,$q$,$m$)} computes the commutative polynomial $$S_m
= \sum_{|j| = m\atop |j| < k} \biggl(\sum_{i = 1}^n
p_i \partial_i(q_j) + p_0 q_j\biggr)y^{j}.$$ This is useful for the
factorization of LPDOs \cite{LasarukSturm:08b}.
\begin{footnotesize}
\begin{verbatim}
input: p := a()***partial(x1)+b()$

input: q := c()***partial(x1)+d()***partial(x2)$

input: lpdos(p,q,1);

a()*diff(c(),x1,1)*y_x1_ + a()*diff(d(),x1,1)*y_x2_ + b()*c()*y_x1_

 + b()*d()*y_x2_
\end{verbatim}
\end{footnotesize}


\subsection{Factorization}

\subsubsection{\texttt{lpdofactorize}}
Factorize the argument LPDO $d$. The ground field $F$ must be fixed via
\texttt{lpdoset}. The result is a list of lists
$\{\dots,(A_i,L_i),\dots\}$. $A_i$ is is genrally the identifiers
\texttt{true}, which indicates reducibility. The respective $L_i$ is a
list of two differential polynomial factors, the first of which has
order 1.

\begin{footnotesize}
\begin{verbatim}
input: bk := (partial(x)+partial(y)+(a10-a01)/2) ***
       (partial(x)-partial(y)+(a10+a01)/2);

                2             2
bk := partial(x)  - partial(y)  + a10 *** partial(x) + a01 *** partial(y)

          2      2
     - a01  + a10
 + ----------------
          4

input: lpdoset lpdoweyl(1,x,y);

{#_00_ + #_01_*y + #_10_*x,x,y}

input: lpdofactorize bk;

{{true,

                                 a01 - a10
  { - partial(x) - partial(y) + -----------,
                                     2

                                  - a01 - a10
    - partial(x) + partial(y) + --------------}}}
                                      2
\end{verbatim}
\end{footnotesize}
If the result is the empty list, then this guarantees that there is no
approximate factorization possible. In general it is possible to obtain
several sample factorizations. Note, however, that the result does not
provide a complete list of possible factorizations with a left factor of
order 1 but only at least one such sample factorization in case of
reducibility.

Furthermore, the procedure might fail due to polynomial degrees
exceeding certain bounds for the extended quantifier elimination by
virtual substitution used internally. In this case there is the
identifier \texttt{failed} returned. This must not be confused with the
empty list indicating irreducibility as described above.

Besides
\begin{enumerate}
\item the LPDO $d$,
\end{enumerate}
\texttt{lpdofactorizex} accepts several optional arguments:
\begin{enumerate}
\setcounter{enumi}{1}
\item An LPDO of order 1, which serves as a template for the left
  (linear) factor. The default is a generic linear LPDO with generic
  coefficient functions according from the ground field specified via
  \texttt{lpdoset}. The principle idea is to support the factorization
  by guessing that certain differential monomials are not present.
\item An LPDO of order $\operatorname{ord}(d)-1$, which serves as a
  template for the right factor. Similarly to the previous argument the
  default is fully generic.
\end{enumerate}

\subsubsection{\texttt{lpdofac}}
This is a low-level entry point to the factorization
\texttt{lpdofactorize}. It accepts the same arguments as
\texttt{lpdofactorize}. It generates factorization conditions as a quite
large first-order formula over the reals. This can be passed to extended
quantifier elimination. For example, consider \texttt{bk} as in the
example for \texttt{lpdofactorize} above:

{\samepage
\begin{footnotesize}
\begin{verbatim}
input: faccond := lpdofac bk$

input: rlqea faccond;

{{true,

               a01 - a10
  {p_00_00_ = -----------,
                   2

   p_00_01_ = 0, p_00_10_ = 0, p_01_00_ = -1, p_01_01_ = 0, p_01_10_ = 0,

   p_10_00_ = -1, p_10_01_ = 0, p_10_10_ = 0,

                - a01 - a10
   q_00_00_ = --------------,
                    2

   q_00_01_ = 0, q_00_10_ = 0, q_01_00_ = 1, q_01_01_ = 0, q_01_10_ = 0,

   q_10_00_ = -1, q_10_01_ = 0, q_10_10_ = 0}}}
\end{verbatim}
\end{footnotesize}
}
The result of the extended quantifier elimination provides coefficient
values for generic factor polynomials $p$ and $q$. These are
automatically interpreted and converted into differential polynomials
by \texttt{lpdofactorize}.

\subsection{Approximate Factorization}

\subsubsection{\texttt{lpdofactorizex}}
Approximately factorize the argument LPDO $d$. The ground field $F$ must
be fixed via \texttt{lpdoset}. The result is a list of lists
$\{\dots,(A_i,L_i),\dots\}$. Each $A_i$ is quantifier-free formula
possibly containing a variable $\texttt{epsilon}$, which describes the
precision of corresponding factorization $L_i$. $L_i$ is a list
containing two factors, the first of which is linear.

\begin{footnotesize}
\begin{verbatim}
input: off lpdocoeffnorm$

input: lpdoset lpdoweyl(0,x1,x2)$

input: f2 := partial(x1)***partial(x2) + 1$

input: lpdofactorizex f2;

{{epsilon - 1 >= 0,{partial(x1),partial(x2)}},

 {epsilon - 1 >= 0,{partial(x2),partial(x1)}}}
\end{verbatim}
\end{footnotesize}
If the result is the empty list, then this guarantees that there is no
approximate factorization possible. In our example we happen to obtain
two possible factorizations. Note, however, that the result in general
does not provide a complete list of factorizations with a left factor of
order 1 but only at least one such sample factorization.

Furthermore, the procedure might fail due to polynomial degrees
exceeding certain bounds for the extended quantifier elimination by
virtual substitution used internally. If this happens, the corresponding
$A_i$ will contain existential quantifiers \texttt{ex}, and $L_i$ will
be meaningless.\\
\verb!Da sollte besser ein failed kommen ...!

The first of the two subresults above has the semantics that
$\partial_{x_1}\partial_{x_2}$ is an approximate factorization of $f_2$
for all $\varepsilon\geq1$. Formally,
$||f_2-\partial_{x_1}\partial_{x_2}||\leq\varepsilon$ for all
$\varepsilon\geq1$, which is equivalent to
$||f_2-\partial_{x_1}\partial_{x_2}||\leq1$. That is, $1$ is an upper
bound for the approximation error over $\R^2$. Where there are two
possible choices for the seminorm $||\cdot||$:
\begin{enumerate}[(a)]
\item ...
\item ...
\end{enumerate}
\verb!explain switch lpdocoeffnorm ...!

Besides
\begin{enumerate}
\item the LPDO $d$,
\end{enumerate}
\texttt{lpdofactorizex} accepts several optional arguments:
\begin{enumerate}
\setcounter{enumi}{1}
\item A Boolean combination $\psi$ of equations, negated equations, and
  (possibly strict) ordering constraints. This $\psi$ describes a
  (semialgebraic) region over which to factorize approximately. The
  default is \texttt{true} specifying the entire $\R^n$. It is possible
  to choose $\psi$ parametrically. Then the parameters will in general
  occur in the conditions $A_i$ in the result.
\item[3., 4.] An LPDO of order 1, which serves as a template for the left
  (linear) factor, and an LPDO of order $\operatorname{ord}(d)-1$, which
  serves as a template for the right factor. See the documentation of
  \texttt{lpdofactorize} for defaults and details.
\item[5.] A bound $\varepsilon$ for describing the desired precision for
  approximate factorization. The default is the symbol \texttt{epsilon},
  i.e., a symbolic choice such that the optimal choice (with respect to
  parameters in $\psi$) is obtained during factorization. It is possible
  to fix $\varepsilon\in\Q$. This does, however, not considerably
  simplify the factorization process in most cases.
\end{enumerate}

\begin{footnotesize}
\begin{verbatim}
input: f3 := partial(x1) *** partial(x2) + x1$

input: psi1 := 0<=x1<=1 and 0<=x2<=1$

input: lpdofactorizex(f3,psi1,a()***partial(x1),b()***partial(x2));

{{epsilon - 1 >= 0,{partial(x1),partial(x2)}}}
\end{verbatim}
\end{footnotesize}

\subsubsection{\texttt{lpdofacx}}
This is a low-level entry point to the factorization
\texttt{lpdofactorizex}. It is analogous to \texttt{lpdofac} for
\texttt{lpdofactorize}; see the documentation there for details.

\subsubsection{\texttt{lpdohrect}}

\subsubsection{\texttt{lpdohcirc}}

\end{document}