qetaaux.spad

-------------------------------------------------------------------
---
--- FriCAS QEta
--- Copyright (C) 2018-2022  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}
\externaldocument{qeta}

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

\begin{document}
\title{Auxiliary functions}
\author{Ralf Hemmecke}
\date{26-Jan-2018}
\maketitle
\begin{abstract}
  The packages \qetatype{QEtaAuxiliaryPackage} contains helper
  functions.
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\tableofcontents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

First of all for debugging.
)endif
OF==>OutputForm
display77(x) ==> display((x::OF)::Formatter(Format1D))
yof x ==> x :: Symbol :: OF
dbgPrint(x,y) ==> display77([yof ":> ", yof x, y::OF]$List(OF))
tracePrint(x,y) ==> display77(hconcat([yof"-- ",yof x,yof":=",y::OF]$List(OF)))
errorPrint(x,y) ==> tracePrint(x,y)

)if LiterateDoc
%$
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
MZZ ==> Matrix ZZ -- will be a square matrix
MQQ ==> Matrix QQ -- matrixEtaOrder
SL2Z ==> MZZ -- represents SL_2(ZZ)
EXGCD ==> Record(coef1: ZZ, coef2: ZZ, generator: ZZ)
Rec ==> Record(red: SL2Z, triang: MZZ)
Pol ==> SparseUnivariatePolynomial ZZ
LZZ ==> List ZZ
LLZZ ==> List LZZ

INTF ==> IntegerNumberTheoryFunctions

asNN x ==> x pretend NN
asPP x ==> x pretend PP

)if LiterateDoc
%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Package QEtaAuxiliaryPackage: Check modularity of
  eta-quotients}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here we mainly deal with the sets $R^*(N)$ and $R^\infty(N)$ from
\cite{HemmeckeRadu_EtaRelations_2019}.

The function \code{matrixLigozat} computes the respective matrix $A_N$
whereas \code{modularGamma0} checks which of the conditions for an $r$
vector to correspond to a modular eta-quotient is violated.


)endif

)abbrev package QETAAUX QEtaAuxiliaryPackage
++ QEtaAuxiliaryPackage helps to do computations with
++ eta-quotients that are modular functions.
QEtaAuxiliaryPackage: with
    matrixAtkinLehner: (PP, PP) -> MZZ
      ++ matrixAtkinLehner(m, t) with t a Hall divisor
      ++ (https://en.wikipedia.org/wiki/Hall_subgroup) of m returns
      ++ the 2x2 matrix with entries [[t, -a], [m, t*b]] such that
      ++ n=m/t and t*b + n*a = 1, see
      ++ \cite{Kohnen_WeierstrassPointsAtInfinity_2004}.
    modularGamma0: (PP, List ZZ) -> ZZ
      ++ modularGamma0(nn, r) returns 0 if all conditions are
      ++ fulfilled. Otherwise it returns a positive number in the
      ++ range 1 to 4 that corresponds to the condition that is not
      ++ met. This corresponds to the conditions given for
      ++ R(N,i,j,k,l) on page 226 of \cite{Radu_RamanujanKolberg_2015}
      ++ and to the conditions \eqref{eq:sum=0},
      ++ \eqref{eq:pure-rhoinfinity}, \eqref{eq:pure-rho0}, and
      ++ \eqref{eq:productsquare} in qeta.tex.
      ++ It is equivalent to check whether there is an extension v of r
      ++ such that matrixModular(nn)*v is 0.
    modularGamma0?: (PP, List ZZ) -> Boolean
      ++ modularGamma0(nn, r) returns true iff the eta-quotient
      ++ corresponding to r is a modular function for Gamma_0(nn). It
      ++ is equivalent to zero?(modularGamma0(nn, r)).
    jacobiLowerStar: (ZZ, ZZ) -> ZZ
      ++ See, for example, Definition 2.26 in the PhD thesis of Silviu Radu.
      ++ http://www3.risc.jku.at/publications/download/risc_5453/main.pdf
    jacobiUpperStar: (ZZ, ZZ) -> ZZ
      ++ See, for example, Definition 2.26 in the PhD thesis of Silviu Radu.
      ++ http://www3.risc.jku.at/publications/download/risc_5453/main.pdf
    upsilonExponent: (ZZ, ZZ, ZZ, ZZ) -> ZZ
      ++ upsilonExponent(a,b,c,d) returns z \in {0..23} such that
      ++ exp(\pi*i*z/12) is a multiplier in the transformation of an
      ++ eta-function by the matrix [[a,b],[c,d]] from SL_2(ZZ). A
      ++ definition can be found in \cite[Lemma~2.27]{Radu_PhD_2010}
      ++ and also as $\kappa_{\gamma_{\delta,m,\lambda}}$ in equation
      ++ \ref{eq:eta_delta-m-lambda(gamma*tau)} of qeta.tex.
    upsilonExponent: MZZ -> ZZ
      ++ upsilonExponent(m) for a matrix m from SL_2(ZZ) returns
      ++ upsilonExponent(m(1,1), m(1,2), m(2,1), m(2,2)).
    numberOfGaps: (PP, List PP) -> NN
      ++ numberOfGaps(n, l) computes
      ++ reduce(+, [floor(x/n) for x in l], 0).
    splitMatrix: (SL2Z, ZZ, PP, NN) -> Rec
      ++ splitMatrix(mat, delta, m, lambda) computes
      ++ splitMatrix(a, b, c, d) where the arguments are the matrix
      ++ entries of matrix [[delta, delta*lambda], [0, m]]*mat.
      ++ Input condition: determinant(mat)>0.
    splitMatrix: (SL2Z, ZZ) -> Rec
      ++ splitMatrix(mat, delta) computes splitMatrix(mat, delta, 1, 0).
      ++ Input condition: determinant(mat)>0.
    splitMatrix: (ZZ, ZZ, ZZ, ZZ) -> Rec
      ++ splitMatrix(a, b, c, d) returns a record rec that contains
      ++ two matrices rec.red=m1 and rec.triang=m2 such that mat=m1*m2
      ++ according to \cide[Lemma~4.8]{Radu_PhD_2010}, but in a variation
      ++ as described in qeta.tex, i.e., such that 0<=m2(1,2)<det(mat)/g
      ++ where g = gcd(a, c) and mat=[[a,b],[c,d]].
      ++ Input condition: a*d-b*c>0.
    fractionalPart: QQ -> QQ
      ++ fractionalPart(x) returns x - floor(x)::Fraction(Integer).
    fractionalBernoulli1: QQ -> QQ
      ++ fractionalBernoulli1(x) computes t - 1/2 where t is the
      ++ fractional part of x, i.e., t = x - floor(x)
    fractionalBernoulli2: QQ -> QQ
      ++ fractionalBernoulli2(x) computes t^2 - t + 1/6 where t is the
      ++ fractional part of x, i.e., t = x - floor(x)
    sawTooth: QQ -> QQ
      ++ sawTooth(x) is 0 if x is an integer and
      ++ fractionalBernoulli1(x) otherwise.
    dedekindSum: (ZZ, ZZ) -> QQ
       ++ dedekindSum(p,q) returns for two relatively prime integers
       ++ p and q the Dedekind sum as defined by
       ++ dedekindSum(p,q)=sum(sawTooth(1/q)*sawTooth(p*i/q), i=1..q-1).
    generalizedDedekindSum: (PP, ZZ, ZZ, ZZ, ZZ) -> QQ
       ++ generalizedDedekindSum(nn,g,h,p,q) returns the generalized
       ++ Dedekind sum as defined in
       ++ \cite[p~673]{Yang:GeneralizedDedekindEtaFunctions} by
       ++ generalizedDedekindSum(nn,g,h,p,q) =
       ++ sum(sawTooth((g+i*nn)/(q*nn))*sawTooth((g'+p*i*nn)/(q*nn)),i=1..q-1).
       ++ Condition: one? gcd(p,q).
    kappaSchoeneberg: (PP, ZZ, ZZ, ZZ, ZZ, ZZ, ZZ) -> QQ
      ++ kappaSchoeneberg(nn, g, h, a, b, c, d) returns a value x such
      ++ that \eta_{g,h}(mat*tau) = exp(2*%pi*%i*x) *
      ++ \eta_{g',h'}(tau) as given at the top of page 673 of
      ++ \cite{Yang_GeneralizedDedekindEtaFunctions_2004} where
      ++ mat=matrix[[a,b],[c,d]].
      ++ See also \eqref{eq:kappa_g-h-N-gamma-Schoeneberg}.
    kappaSchoeneberg: (PP, ZZ, ZZ, MZZ) -> QQ
      ++ kappaSchoeneberg(nn, g, h, mat) returns a value x such
      ++ that \eta_{g,h}(mat*tau) = exp(2*%pi*%i*x) *
      ++ \eta_{g',h'}(tau) as given at the top of page 673 of
      ++ \cite{Yang_GeneralizedDedekindEtaFunctions_2004}.
      ++ See also \eqref{eq:kappa_g-h-N-gamma-Schoeneberg}.
    kappaYang: (PP, ZZ, ZZ, ZZ, ZZ, ZZ, ZZ) -> QQ
      ++ kappaYang(nn, g, h, a, b, c, d) returns a value x such that
      ++ E_{g,h}(mat*tau) = exp(2*%pi*%i*x) * E_{g',h'}(tau) as given
      ++ in \cite[Theorem~1]{Yang_GeneralizedDedekindEtaFunctions_2004}
      ++ where mat = matrix [[a,b],[c,d]].
      ++ Condition: not zero? c.
    kappaYang: (PP, ZZ, ZZ, MZZ) -> QQ
      ++ kappaYang(nn, g, h, mat) returns a value x such that
      ++ E_{g,h}(mat*tau) = exp(2*%pi*%i*x) * E_{g',h'}(tau) as given
      ++ in \cite[Theorem~1]{Yang_GeneralizedDedekindEtaFunctions_2004}.
      ++ Condition: not zero? mat(2,1).
    alphaSchoenebergContribution: (QQ, QQ) -> QQ
      ++ alphaSchoenebergContribution(a,b) returns
      ++ fractionalBernoulli1(b) if a is an integer and b is different
      ++ from zero, otherwise it returns 1.
      ++ See \eqref{eq:alpha-Schoeneberg-Contribution} in qeta.tex.
    minimizeVector: (Vector ZZ, List Vector ZZ) -> Vector ZZ
      ++ minizeVector(v, basis) computes
      ++ vps := concat [w for b in bas | (w:=v+b; dot(w,w)<dot(v,v)]
      ++ vms := concat [w for b in bas | (w:=v-b; dot(w,w)<dot(v,v)]
      ++ From these vectors it takes the one with minimal length as
      ++ the new v and iterates as long as concat(vps, vms) is not empty.
      ++ This minimal vector is eventually returned.
    primePower: (PP, PP) -> NN
      ++ primePower(d, p) returns e such that gcd(d/p^e, p)=1, i.e., the
      ++ highest power e such that p^e is a factor of d.
    verticalConcat: (Matrix QQ, Matrix QQ) -> Matrix QQ
      ++ verticalConcat(m1, m2) puts m1 on top of m2 while creating a
      ++ matrix big enough to hold all columns.
  == add

)if LiterateDoc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Related to Eta-Quotients}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

)endif

    primePower(d: PP, p: PP): NN ==
        e: NN := 0;
        x: ZZ := d
        while not one? gcd(x, p) repeat (x := (x exquo p)::ZZ; e:=e+1)
        return e

    matrixAtkinLehner(m: PP, t: PP): MZZ ==
        -- t is a Hall divisor of m.
        -- https://en.wikipedia.org/wiki/Hall_subgroup
        --assert((m exquo t) case ZZ)
        n: ZZ := (m exquo t)::ZZ -- input condition is that t divides m
        --assert(one? gcd(n, t))
        exgcd: EXGCD := extendedEuclidean(n, t)
        matrix [[t, -exgcd.coef1], [m, t*exgcd.coef2]]

    -- These are the conditions (1) to (4) of
    -- \cite{HemmeckeRadu_EtaRelations_2019} that have to be fulfilled
    -- by a vector r to lie in R^*(NN).
    modularGamma0(m: PP, r: List ZZ): ZZ ==
        divs: List ZZ := divisors(m)$INTF
        rdivs: List ZZ := reverse divs
        s: ZZ := 0
        a: ZZ := 0
        b: ZZ := 0
        c: Factored ZZ := 1
        for ri in r for d in divs for rd in rdivs repeat
            s := s + ri
            a := a + ri * d  -- sigma_\infty = 24*\rho_\infty
            b := b + ri * rd -- sigma_0 = 24*\rho_0
            c := c * factor(d)^asNN(abs ri)
        -- Check conditions 1, 2, 3, 4.
        not zero? s => 1
        not zero? positiveRemainder(a, 24) => 2
        not zero? positiveRemainder(b, 24) => 3
        for fe in factors c repeat if odd?(fe.exponent) then return 4
        return 0

    modularGamma0?(nn: PP, r: List ZZ): Boolean ==
        zero? modularGamma0(nn, r)

)if LiterateDoc
%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Related to Transformations of Eta-Quotients}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

)endif

    jacobiLowerStar(c: ZZ, d: ZZ): ZZ == -- result either 1 or -1
        --assert(one? gcd(c,d))
        --aaawer(odd? d)
        -- See \cite[Def.~2.26]{Radu_PhD_2010}.
        -- \left(\frac{c}{d}\right)_* =
        -- \begin{cases}
        -- \left(\frac{c}{\lvert d \rvert}\right) (-1)^{m_c m_d},
        --   & \text{if $c\ne0$},\\
        -- \sign{d}, & \text{otherwise}
        -- \end{cases}
        -- where $m x=\frac{\sign{x}-1}{2}$.
        zero? c => sign(d) -- note that d ~= 0 in this case
        j: ZZ := jacobi(c, abs(d))$INTF
        if c>0 or d>0 then j else -j -- c=0 or d=0 does not happen

    jacobiUpperStar(c: ZZ, d: ZZ): ZZ == -- result either 1 or -1
        --assert(one? gcd(c,d))
        --assert(odd? d)
        -- See \cite[Def.~2.26]{Radu_PhD_2010}.
        -- \left(\frac{c}{d}\right)^* =
        -- \begin{cases}
        -- \left(\frac{c}{\lvert d \rvert}\right), & \text{if $c\ne0$},\\
        -- 1, & \text{otherwise}
        -- \end{cases}
        zero? c => 1
        jacobi(c, abs(d))$INTF

    upsilonExponent(a: ZZ, b: ZZ, c: ZZ, d: ZZ): ZZ ==
        u: ZZ := c*(a+d)-b*d*(c^2-1)-3*c
        v: ZZ := 3*(d-1)*(c-1)
        r(x) ==> positiveRemainder(x, 24)
        -- Note that we exp(2*\pi*i*n/24)=exp(\pi*i)=-1 for n=12.
        odd? c => if jacobiUpperStar(d, c) > 0 then r(u) else r(u+12)
        if jacobiLowerStar(c, d) > 0 then r(u-v) else r(u-v+12)

    upsilonExponent(mat: MZZ): ZZ ==
        upsilonExponent(mat(1,1), mat(1,2), mat(2,1), mat(2,2))

    numberOfGaps(n: PP, list: List PP): NN ==
        asNN reduce(_+, [floor(x /$QQ n) for x in list]$List(ZZ), 0)

    -- Split matrix as in Lemma~\ref{thm:matrix-splitting} with
    -- z = delta*lambda.
    splitMatrix(mat: SL2Z, delta: ZZ, m: PP, lambda: NN): Rec ==
        z: ZZ := delta*lambda
        w11: ZZ := delta*mat(1,1) + mat(2,1)*z
        w12: ZZ := delta*mat(1,2) + mat(2,2)*z
        w21: ZZ := mat(2,1)*m
        w22: ZZ := mat(2,2)*m
        splitMatrix(w11, w12, w21, w22)

    -- Split matrix as in Lemma~\ref{thm:matrix-splitting} with
    -- z = 0 and m = 1.
    splitMatrix(mat: SL2Z, delta: ZZ): Rec ==
        splitMatrix(mat(1,1) * delta, mat(1,2) * delta, mat(2,1), mat(2,2))
        -- Same as splitMatrix(mat, delta, 1, 0).

    -- Split matrix according to Section~\ref{sec:eta-transformation} in
    -- qeta.tex.
    splitMatrix(a: ZZ, b: ZZ, c: ZZ, d: ZZ): Rec ==
        exgcd: EXGCD := extendedEuclidean(a, c)
        g: ZZ := exgcd.generator
        d0: ZZ :=  exgcd.coef1 -- initial value for d1
        b0: ZZ := -exgcd.coef2 -- initial value for b1
        a1: ZZ := (a exquo g) :: ZZ
        c1: ZZ := (c exquo g) :: ZZ
        a2: ZZ := g
        d2: ZZ := a1*d-b*c1
        -- Now try to find s such that 0 <= b*d1 - d*b1 - s*d2 < d2.
        b2init: ZZ := b*d0 - d*b0
        b2: ZZ := positiveRemainder(b2init, d2)
        s: ZZ := ((b2init - b2) exquo d2) :: ZZ
        b1: ZZ := b0 + s * a1
        d1: ZZ := d0 + s * c1
        m1: MZZ := matrix([[a1, b1], [c1, d1]]$LLZZ)
        m2: MZZ := matrix([[g, b2], [0, d2]]$LLZZ)
        [m1, m2]$Rec

)if LiterateDoc
%$

We use \cite[Lemma~2.45]{Radu_PhD_2010} to compute the cusps of
$\Gamma_0(N)$ as the set of all $\frac{a}{c}$ such that $c|N$ and
$a\in X_c$.

)endif

    -- local
    -- TODO: unused
    square(x: Record(root: PP, elem: Pol)): Vector ZZ ==
        pol: Pol := cyclotomic(x.root)$CyclotomicPolynomialPackage
        C ==> SimpleAlgebraicExtension(ZZ, Pol, pol)
        z: C := convert(x.elem)@C
        convert(z*z)@Vector(ZZ)

)if LiterateDoc
%$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Related to Generalized Eta-Quotients}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

)endif

    fractionalPart(x: QQ): QQ == -- local function
        d: ZZ := denom x
        positiveRemainder(numer x, d)/d

    fractionalBernoulli1(x: QQ): QQ == fractionalPart x - 1/2
    fractionalBernoulli2(x: QQ): QQ ==
        t: QQ := fractionalPart x
        t^2 - t + 1/6

    sawTooth(x: QQ): QQ ==
        one? denom x => 0
        fractionalBernoulli1 x
    dedekindSum(p: ZZ, q: ZZ): QQ ==
        not one? gcd(p, q) => error "dedekindSum: input not coprime"
        ds: QQ := 0
        for i in 1..abs(q)-1 repeat
            ds := ds + sawTooth(i/q)*sawTooth(p*i/q)
        return ds
    generalizedDedekindSum(nn: PP, g: ZZ, h: ZZ, a: ZZ, c: ZZ): QQ ==
        not one? gcd(a, c) => error "generalizedDedekindSum: input not coprime"
        gds: QQ := 0
        gg: ZZ := a*g + c*h
        for i in 0..abs(c)-1 repeat
            cn := c * nn
            s1: QQ := sawTooth((g+i*nn)/cn)
            s2: QQ := sawTooth((gg+a*i*nn)/cn)
            gds := gds + s1*s2
        return gds

    kappaSchoeneberg(nn: PP, g: ZZ, h: ZZ, a: ZZ, b: ZZ, c: ZZ, d: ZZ): QQ ==
        zero? c => fractionalPart(b/(2*d)*fractionalBernoulli2(g/nn))
        gprime: ZZ := a*g + c*h
        x := a/(2*c)*fractionalBernoulli2(g/nn) _
             + d/(2*c)*fractionalBernoulli2(gprime/nn)
        gds := generalizedDedekindSum(nn, g, h, a, c)
        c < 0 => fractionalPart(x + gds)
        fractionalPart(x - gds)

    kappaYang(nn: PP, g: ZZ, h: ZZ, a: ZZ, b: ZZ, c: ZZ, d: ZZ): QQ ==
        r: QQ := 1/2* ((g^2*a*b + 2*g*h*b*c + h^2*c*d)/nn^2 _
           - (g*b + h*(d-1))/nn) -- corresponds to delta/2 in Yang's paper
        u: ZZ := c*(a+d-3)+b*d*(1-c^2)
        odd? c => fractionalPart(u/12 + r)
        even? d => error "kappaYang: d is even"
        v: ZZ := a*c*(1-d^2)+d*(b-c+3)
        minusI: QQ := -1/4
        fractionalPart(minusI + v/12 + r)

    kappaSchoeneberg(nn: PP, g: ZZ, h: ZZ, mat: MZZ): QQ ==
        kappaSchoeneberg(nn, g, h, mat(1,1), mat(1,2), mat(2,1), mat(2,2))

    kappaYang(nn: PP, g: ZZ, h: ZZ, mat: MZZ): QQ ==
        kappaYang(nn, g, h, mat(1,1), mat(1,2), mat(2,1), mat(2,2))

    alphaSchoenebergContribution(a: QQ, b: QQ): QQ ==
        not one? denom a or one? denom b => 1
        (1/2)*fractionalBernoulli1(b)

    -- By adding interger multiples of elements of b we try to make
    -- the length of a vector v as small as possible in a cheap way.
    -- The emphasis here is on "cheap". It is actually an auxiliary
    -- function. It is not completely necessary to find the optimally
    -- minimal vector, since there is a second optimization step in
    -- the place where we are going to use it to compute the cofactor
    -- specification.
    minimizeVector(v: Vector ZZ, basis: List Vector ZZ): Vector ZZ ==
        empty? basis => v
        w := v; lw := lv := dot(v, v)
        repeat
            for b in basis repeat -- search for a minimal w
                t := v+b; lt := dot(t, t); if lt < lw then (w := t; lw := lt)
                t := v-b; lt := dot(t, t); if lt < lw then (w := t; lw := lt)
            lw = lv => break -- no change happened
            -- otherwise lw < lv
            v := w; lv := lw -- take this minimal w as the new v
        w

    verticalConcat(mat1: MQQ, mat2: MQQ): MQQ ==
        r1 := nrows mat1; c1 := ncols mat1
        r2 := nrows mat2; c2 := ncols mat2
        mat: MQQ := new(r1 + r2, max(c1, c2), 0)
        for i in 1..r1 repeat for j in 1..c1 repeat mat(i, j) :=  mat1(i, j)
        for i in 1..r2 repeat for j in 1..c2 repeat mat(r1+i, j) :=  mat2(i, j)
        mat













)if LiterateDoc
%$
\end{document}
)endif