qetatool.spad

-------------------------------------------------------------------
---
--- FriCAS QEta
--- Copyright (C) 2015-2021  Ralf Hemmecke <ralf@hemmecke.org>
---
-------------------------------------------------------------------
-- This program is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-- GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with this program.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------

)if LiterateDoc
\documentclass{article}
\usepackage{qeta}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
\title{Conversion, Evaluation, and Gr\"obner bases} %"
\author{Ralf Hemmecke}
\date{10-Dec-2015}
\maketitle
\begin{abstract}
  Auxiliary packages to help with the implementation of algorithms
  from \cite{Radu_RamanujanKolberg_2015},
  \cite{Hemmecke_DancingSambaRamanujan_2018}, and
  \cite{HemmeckeRadu_EtaRelations_2019}.


  This file contains a number of auxiliary packages, namely
  \begin{itemize}
  \item \qetatype{QAuxiliaryToos} helps to create variables (symbols) with
    given indices and converts a polynomial from
    \spadtype{Polynomial}\code{(QQ)} to \spadtype{Polynomial}\code{(ZZ)}
    by clearing denominators (where \code{ZZ==>}\spadtype{Integer} and
    \code{QQ==>}\spadtype{Fraction}\code{(ZZ)}).

  \item \qetatype{QEtaLaurentSeriesFunctions2} helps to embed elements of
    \qetatype{QEtaLaurentSeries}(C1) into \qetatype{QEtaLaurentSeries}(C2).

  \item \qetatype{Monomials} is introduced to add variable names so
    that an element of a \spadtype{DirectProduct} can be shown as a
    product of variables raised to some power. It basically turns the
    additive structure of \spadtype{DirectProduct} into a
    multiplicative structure just for output.

  \item \qetatype{PolynomialConversion} converts a polynomial from the
    generic \spadtype{Polynomial}\code{(C)} to a more specific
    polynomial ring of the form \qetatype{PolynomialRing}{(C, E)} with a
    given domain of the exponents.

  \item \qetatype{PolynomialEvaluation} evaluates the variables of a
    (multivariate) polynomial by given values in a structure that
    comes with the operations \code{+}, \code{*}, and \code{^}.

  \item \qetatype{PolynomialTool} helps to extract from a number of
    polynomials only those whose variables start with a certain
    letter.
    %
    This package helps to compute $I\cap C[E]$ where
    $I \subset C[Y,E]$ and $Y=(Y_1,\ldots,Y_n$),
    $E=(E_1,\ldots, E_n)$.

  \item \qetatype{QEtaGroebner} is a wrapper for the FriCAS
    \spadtype{GroebnerPackage}.
  \end{itemize}
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\tableofcontents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Overview}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This code implements a number of conversion and evaluation tools.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Helper macros}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let us start with a few common macros.

These two technical macros are necessary to distinguish between Rep
and \%.
)endif

rep x ==> (x@%) pretend Rep
per x ==> (x@Rep) pretend %

)if LiterateDoc
Now some abbreviations for common domains.
)endif

PP ==> PositiveInteger
NN ==> NonNegativeInteger
ZZ ==> Integer
QQ ==> Fraction ZZ

)if LiterateDoc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package QAuxiliaryTools}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We want variables that are easy to parse when written to a file.
)endif

)abbrev package QAUXTOOL QAuxiliaryTools
++ Miscellaneous tools to deal with polynomials.
QAuxiliaryTools(): with
    indexedSymbol: (String, List ZZ) -> Symbol
        ++ indexedSymbol(s, l) returns a symbol with a name that
        ++ starts with the string s and is directly followed by the
        ++ numbers of l separated by underscores ("__").
    indexedSymbols: (String, NN) ->  List Symbol
        ++ indexedSymbols(s, n) returns indexedSymbols(s, [i for i in 1..n])
    indexedSymbols: (String, List ZZ) ->  List Symbol
        ++ indexedSymbols(s, l) returns indexedSymbols(s, [[n] for n in l]).
    indexedSymbols: (String, List List ZZ) ->  List Symbol
        ++ indexedSymbols(s, ll) returns [indexedSymbol(s,l) for l in ll].
    clearDenominator: Polynomial QQ -> Polynomial ZZ
        ++ clearDenominator(p) multiplies the polynomial p over rational numbers
        ++ with the least common multiple of all its coefficients and
        ++ returns the result as a polynomial over the integers.
  == add
    indexedSymbol(s: String, l: List ZZ): Symbol ==
        empty? l => s :: Symbol
        str: String := concat(s, convert(first l)@String)
        for i in rest l repeat
            a: String := concat("__", convert(i)@String)
            str := concat(str, a)
        str :: Symbol
    indexedSymbols(s: String, n: NN): List Symbol ==
        [indexedSymbol(s, [i]) for i in 1..n]
    indexedSymbols(s: String, l: List ZZ): List Symbol ==
        [indexedSymbol(s, [i]) for i in l]
    indexedSymbols(s: String, ll: List List ZZ): List Symbol ==
        [indexedSymbol(s, l) for l in ll]
    clearDenominator(p: Polynomial QQ): Polynomial ZZ ==
        zero? p => 0$Polynomial(ZZ)
        c: ZZ := lcm [denom x for x in coefficients p]
        mon(x) ==> monomial(numer(c*leadingCoefficient x), degree x)$Polynomial(ZZ)
        lm: List Polynomial ZZ := [mon(x) for x in monomials p]
        reduce(_+, lm)

)if LiterateDoc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package QEtaLaurentSeriesFunctions2}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Convert convert power series by mapping coefficients.
)endif

)abbrev package QETALS2 QEtaLaurentSeriesFunctions2
++ QEtaLaurentSeriesFunctions2 maps coefficients of the series.
QEtaLaurentSeriesFunctions2(C1, C2): Exports == Implementation where
  C1: CommutativeRing
  C2: CommutativeRing
  L1 ==> QEtaLaurentSeries C1
  L2 ==> QEtaLaurentSeries C2
  Exports ==> with
    map: (C1 -> C2, L1) -> L2
        ++ \spad{map(f, g(x))} applies the map f to the coefficients of
        ++ the series \spad{g(x)}.
  Implementation ==> add
    map(f: C1 -> C2, s: L1): L2 ==
        cs: Stream C1 := coefficients qetaTaylorRep s
        t2: QEtaTaylorSeries C2 := series(map(f, cs)$StreamFunctions2(C1, C2))
        laurent(degree s, t2)$L2

)if LiterateDoc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Domain Monomials: An exponent domain with variable names}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Turn various kinds of direct products into something that prints as
powerproducts with exponent vectors being the entries of the direct
product element.
)endif

)abbrev domain MONOMS Monomials
++ Monomials(dim, R, D, vl) behaves exactly like D, i.e. is an
++ AbelianMonoid, but prints its elements in a multiplicative form.
++ For example, if dim = 2, vl = [A,B] and x::Vector(NNI) = [2,3],
++ then the element x (coerced to OutputForm) looks like A^2*B^3.
Monomials(_
  dim: NN,_
  R: OrderedAbelianMonoid,
  D: DirectProductCategory(dim, R),_
  vl: List Symbol_
 ): DirectProductCategory(dim, R) == D add
    --assert(dim = # vl)
    noOneTest? := not (R has one?: R -> Boolean)
    notOne? x ==> noOneTest? or not one? x
    coerce(x: %): OutputForm ==
        vs: List Symbol := vl
        zero? x => (1$Integer)::OutputForm
        fst: Boolean := true
        k: NN := 1
        while zero?(x.k) repeat
            vs := rest vs
            k := k + 1
        --assert(k<=dim)
        o: OutputForm := (first vs)::OutputForm
        if notOne? x.k then o := o ^ ((x.k)::OutputForm)
        for i in k+1..dim for v in rest vs | not zero?(x.i) repeat
            oo: OutputForm := v::OutputForm
            if notOne? x.i then oo := oo ^ ((x.i)::OutputForm)
            o := o * oo
        return o


)if LiterateDoc
%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package PolynomialConversion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\qetatype{PolynomialConversion} converts a polynomial from the generic
\spadtype{Polynomial}\code{(C)} to a more specific polynomial ring of
the form \spadtype{PolynomialRing}\code{(C, E)} with a given domain of the
exponents.
)endif

)abbrev package POLYCONV PolynomialConversion
PolynomialConversion(_
  C: Ring, _
  E: OrderedAbelianMonoid with (
       directProduct: Vector NN -> %; members: % -> List NN),_
  syms: List Symbol): with
    coerce: Polynomial C -> PolynomialRing(C, E)
    coerce: PolynomialRing(C, E) -> Polynomial C
 == add
    R ==> PolynomialRing(C, E)
    coerce(p: Polynomial C): R ==
        r: R := 0
        while not zero? p repeat
            c: C := leadingCoefficient p
            v: Vector NN := vector degree(leadingMonomial p, syms)
            e: E := directProduct(v)$E
            p := reductum p
            r := r + monomial(c, e)
        r

    coerce(r: R): Polynomial C ==
        p: Polynomial(C) := 0
        while not zero? r repeat
            c: C := leadingCoefficient r
            e: E := leadingSupport r
            r := reductum r
            exponents: List NN := members e
            p := p + monomial(c::Polynomial(C), syms, exponents)
        p

)if LiterateDoc
%$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package PolynomialEvaluation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Evaluate a polynomial at a certain point.
)endif

)abbrev package POLYEVAL PolynomialEvaluation
PolynomialEvaluation(_
  C: Ring, _
  S: with (_
      _+: (%, %) -> %; _
      _*: (%, %) -> %; _
      _^: (%, NonNegativeInteger) -> %)_
 ): with
    eval: (Polynomial C, C -> S, List Symbol, List S) -> S
    eval: (C -> S, List Symbol, List S) -> Polynomial C -> S
 == add
    eval(p: Polynomial C, embed: C -> S, vars: List Symbol, vals: List S): S ==
        E ==> IndexedExponents Symbol
        PE ==> PolynomialCategoryLifting(E, Symbol, C, Polynomial C, S)
        map((s:Symbol):S +-> vals.position(s, vars), embed, p)$PE
    eval(embed: C -> S, vars: List Symbol, vals: List S): Polynomial C -> S ==
      (p: Polynomial C): S +-> eval(p, embed, vars, vals)

)if LiterateDoc
%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package PolynomialTool}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The following package helps to compute $I\cap C[E]$ where $I \subset
C[Y,E]$ and $Y=(Y_1,\ldots,Y_n$), $E=(E_1,\ldots, E_n)$.
)endif

)abbrev package POLYAUX PolynomialTool
PolynomialTool(R: Ring): with
    xPolynomials: (List Polynomial R, Character) -> List Polynomial R
      ++ xPolynomials(pols, c) returns all polynomials p from
      ++ pols such that variables(p) contains only variables starting
      ++ with the character c.
 == add
    -- variables not starting with the character c.
    xVariables(p: Polynomial(R), c: Character): List Symbol ==
        [x for x in variables p | (string x).1 ~= c]

    -- Is it a polynomial just in variables starting with the character c?
    xPolynomial?(p: Polynomial R, c: Character): Boolean ==
        degs: List NN := degree(p, xVariables(p, c))
        zero? reduce(_+, degs, 0)

    -- Extract all polynomials just in variables that begin with the
    -- character c.
    xPolynomials(pols: List Polynomial R, c: Character): List Polynomial R ==
        [x for x in pols | xPolynomial?(x, c)]

)if LiterateDoc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package QEtaGroebner}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A wrapper for \GB{} omputations.
)endif

)abbrev package QETAGB QEtaGroebner
++ QEtaGroebner(C, E) is a wrapper around the GrobnerPackage. It basically
++ converts given polynomials to elements in PolynomialRing(C, E), does
++ some Groebner basis computation or reduction and then converts back.
QEtaGroebner(C, E): Exports == Implementation where
  C: GcdDomain
  ExponentCat ==> OrderedAbelianMonoidSup with _
    (directProduct: Vector NN -> %; members: % -> List NN)
  E: ExponentCat
  Pol ==> Polynomial C
  LPol ==> List Pol
  R ==> PolynomialRing(C, E)
  X ==> Record(poly: R, repr: Vector R, mult: C)
  Exports ==> with
    groebner: (LPol, List Symbol) -> LPol
      ++ groebner(lpol, syms) considers the list of polynomials as
      ++ polynomials in syms, and computes a Groebner basis with
      ++ respect to the order given by E.
      ++ We assume that each of the polynomials in lpol is indeed a
      ++ polynomial in syms (and no other variables) and that E
      ++ corresponds to exactly #syms variables.
    groebnerExtend: (LPol, LPol, List Symbol) -> LPol
      ++ groebnerExtend(lpol, gb, syms) computes a Groebner basis of the
      ++ union of lpos and gb under the assumption that gb is already
      ++ a Groebner basis wrt. the order given by E.
      ++ We assume that each of the polynomials in lpol and gb is indeed a
      ++ polynomial in syms (and no other variables) and that E
      ++ corresponds to exactly #syms variables.
    normalForms: (LPol, LPol, List Symbol) -> LPol
      ++ normalForms(lpol, gb, syms) reduces each polynomial from lpol
      ++ with respect to the Groebner basis given by gb.
      ++ We assume that each of the polynomials in lpol and gb is indeed a
      ++ polynomial in syms (and no other variables) and that E
      ++ corresponds to exactly #syms variables.
    extendedNormalForm: (Pol, LPol, List Symbol, Symbol, String) -> Pol
      ++ extendedNormalForm(pol, gb, syms, f, g) returns
      ++ extendedNormalForm(pol, gb, syms, f, gsyms) for
      ++ gsyms := indexedSymbols(g, #gb).
    extendedNormalForm: (Pol, LPol, List Symbol, Symbol, List Symbol) -> Pol
      ++ extendedNormalForm(pol, gb, syms, f, gsyms) reduces pol with
      ++ respect to gb and returns that reduced polynomial together
      ++ with its relations in terms of gb. The original polynomial
      ++ is represented by the variable f and the Groebner basis
      ++ elements by the variables gsyms.
      ++ We assume that each of the polynomials pol and gb is indeed a
      ++ polynomial in syms (and no other variables) and that E
      ++ corresponds to exactly #syms variables.
    extendedNormalForms: (LPol, LPol, List Symbol, String, String) -> LPol
      ++ extendedNormalForms(lpol, gb, syms, f, g) returns
      ++ [extendedNormalForm(p, gb, syms, f, g)
      ++   for p in lpol for f in indexedSymbols(f, #lpol)]
  Implementation ==> add
    -- local
    -- toPol(x, syms, sym, gsyms) returns a polynomial that describes
    -- the representation of x.poly in terms of x.repr
    toPol(x: X, syms: List Symbol, sym: Symbol, gsyms: List Symbol): Pol ==
        import from PolynomialConversion(C, E, syms)
        p: Pol := (x.mult)*(sym::Pol) - (x.poly)::Pol
        for gsym in gsyms for r in members(x.repr) repeat
            p := p - (gsym::Pol)*(r::Pol)
        return primitivePart p

    -- exported
    groebner(lpol: LPol, syms: List Symbol): LPol ==
        import from PolynomialConversion(C, E, syms)
        rs: List R := [p::R for p in lpol]
        gb: List R := groebner(rs)$GroebnerPackage(C, E, R)
        [r::Pol for r in gb]

    groebnerExtend(lpol: LPol, gb: LPol, syms: List Symbol): LPol ==
        import from PolynomialConversion(C, E, syms)
        import from GroebnerInternalPackage(C, E, R)
        rs: List R := [p::R for p in lpol]
        rgb: List R := [p::R for p in gb]
        xgb: List R := gbasisExtend(rs, rgb, 2, 1)
        mgb: List R := minGbasis(sort((x, y) +-> degree x > degree y, xgb))
        [r::Pol for r in mgb]

    normalForms(lpol: LPol, gb: LPol, syms: List Symbol): LPol ==
        import from PolynomialConversion(C, E, syms)
        import from GroebnerInternalPackage(C, E, R)
        rs: List R := [p::R for p in lpol]
        rgb: List R := [p::R for p in gb]
        [(primitivePart redPol(x, rgb))::Pol for x in rs]

    extendedNormalForm(pol: Pol, gb: LPol, syms: List Symbol, f: Symbol, g: String): Pol ==
        gsyms: List Symbol := indexedSymbols(g, #gb)$QAuxiliaryTools
        extendedNormalForm(pol, gb, syms, f, gsyms)

    extendedNormalForm(pol: Pol, gb: LPol, syms: List Symbol, f: Symbol, gsyms: List Symbol): Pol ==
        import from PolynomialConversion(C, E, syms)
        import from ExtendedPolynomialReduction(C, E, R)
        rgb: List R := [p::R for p in gb]
        x: X := reduce(pol::R, rgb)
        toPol(x, syms, f, gsyms)

    extendedNormalForms(lpol: LPol, gb: LPol, syms: List Symbol, f: String, g: String): LPol ==
        import from PolynomialConversion(C, E, syms)
        import from ExtendedPolynomialReduction(C, E, R)
        rs: List R := [p::R for p in lpol]
        rgb: List R := [p::R for p in gb]
        lx: List X := [reduce(r, rgb) for r in rs]
        import from QAuxiliaryTools
        fsyms: List Symbol := indexedSymbols(f, #lpol)
        gsyms: List Symbol := indexedSymbols(g, #gb)
        [toPol(x, syms, fsym, gsyms) for fsym in fsyms for x in lx]


)if LiterateDoc
\end{document}
)endif