% ------------------------------------------------ % This is an AMSLaTeX file for LaTeX 2e [16 pages] % ------------------------------------------------ \documentclass[a4paper]{article} \usepackage{amsmath} \usepackage{amsthm} \usepackage{amssymb} \setlength{\arraycolsep}{0pt} \addtolength{\oddsidemargin}{-1cm} \addtolength{\textwidth}{2.2cm} \addtolength{\medskipamount}{3pt} \renewcommand{\labelenumi}{(\roman{enumi})} \renewcommand{\thefootnote}{\fnsymbol{footnote}} \newtheoremstyle{thmsty}{\medskipamount}{\medskipamount}% {\it}{}{\bfseries}{.}{.5em}{} \newtheoremstyle{dfnsty}{\medskipamount}{\medskipamount}% {}{}{\bfseries}{.}{.5em}{} \newtheoremstyle{algsty}{\medskipamount}{\bigskipamount}% {\it}{}{\bfseries}{.}{.5em}% {Algorithm \mathversion{bold}\thmnote{#3}\mathversion{normal}} \theoremstyle{thmsty} \newtheorem{thm}{Theorem} \newtheorem{lem}{Lemma} \theoremstyle{dfnsty} \newtheorem{dfn}{Definition} \newtheorem*{exa}{Example} \theoremstyle{algsty} \newtheorem*{alg}{} \newcommand{\Fp}{F[x,y]} \newcommand{\Fr}{F(x,y)} \newcommand{\N}{\mathbb{N}} \newcommand{\Pc}{\mathcal{P}} \newcommand{\Uc}{\mspace{2mu}\mathcal{U}} \newcommand{\Vc}{\mathcal{V}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\bi}{bibasic\ } \newcommand{\hg}{hypergeometric\ } \newcommand{\Telpq}{Telescope$_{p,q}$} \newcommand{\alphap}{{\alpha_{\mspace{-1mu}+}}} \newcommand{\alpham}{{\alpha_{\mspace{-1mu}-}}} \newcommand{\betap}{{\beta_{\mspace{-2mu}+}}} \newcommand{\betam}{{\beta_{\mspace{-2mu}-}}} \newcommand{\bfbar}[1]{\mathbf{\bar{\mathnormal #1}}} % makes bold bar \newcommand{\bmp}[1]{#1^{\mspace{-1mu}\ast\mspace{-2mu}}} % monic part \newcommand{\fk}{(f_k)_{k \in \Z}} \newcommand{\gk}{(g_k)_{k \in \Z}} \newcommand{\gcdpq}{\gcd{}_{\!p,q}} \newcommand{\GFFpq}{\GFF_{\!p,q}} \newcommand{\pnb}{$p$\nobreakdash-\hspace{0pt}} % don't allow linebreak \newcommand{\qnb}{$q$\nobreakdash-\hspace{0pt}} \newcommand{\qfac}[2]{[#1]_{p,q}^{\underline{#2}}} \newcommand{\VMULT}{$\mathrm{\bmp{V}}$MULT} \DeclareMathOperator{\diag}{diag} \DeclareMathOperator{\GFF}{GFF} \DeclareMathOperator{\res}{Res} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} \title{A Generalization of Gosper's Algorithm\\ to Bibasic Hypergeometric Summation} \author{Axel Riese\\ \\ Research Institute for Symbolic Computation\\ Johannes Kepler University Linz\\ A--4040 Linz, Austria\\ \texttt{Axel.Riese@risc.uni-linz.ac.at}\\[.3cm]} \date{Submitted: May 9, 1996; Accepted: June 24, 1996.} \maketitle \pagestyle{myheadings} \markright{\sc the electronic journal of combinatorics 3 (1996), \#R19\hfill} \thispagestyle{empty} \begin{abstract} \noindent An algebraically motivated generalization of Gosper's algorithm to indefinite bibasic hypergeometric summation is presented. In particular, it is shown how Paule's concept of greatest factorial factorization of polynomials can be extended to the bibasic case. It turns out that most of the bibasic hypergeometric summation identities from literature can be proved and even found this way. A Mathematica implementation of the algorithm is available from the author. \end{abstract} \vspace{12pt} \textbf{AMS Subject Classification.} Primary 33D65, 68Q40; Secondary 33D20. \section{Introduction} % ------------ Recently, Paule and Strehl~\cite{PauleStrehl:SymbSumm} observed that the algorithm presented by Gosper~\cite{Gosper:DecProc} for indefinite \hg summation extends quite naturally to the \qnb \hg case by introducing a \qnb analogue of the canonical Gosper-Petkov\v{s}ek (GP) representation for rational functions. Based on the new algebraic concept of greatest factorial factorization (GFF), Paule~\cite{Paule:GFF} developed an alternative but equivalent approach to \hg telescoping. It was also shown by Paule (cf.~Paule and Riese~\cite{PauleRiese:qZeil}) that the problem of \qnb\hg telescoping can be treated along the same lines as the $q=1$ case by making use of a \qnb version of GFF. Built on these concepts, a Mathematica implementation of \qnb analogues of Gosper's as well as of Zeilberger's~\cite{Zeilberger:FastAlg} fast algorithm for definite \qnb \hg summation has been carried out by the author (cf.\ Paule and Riese~\cite{PauleRiese:qZeil}, and Riese~\cite{Riese:Dipl}). The original approach to definite \qnb \hg summation is due to Wilf and Zeilberger~\cite{WilfZeilberger:AlgPrTh}. The object of this paper is to describe how the algorithm $q$Telescope presented in \cite{PauleRiese:qZeil}, a \qnb analogue of Gosper's algorithm, generalizes to the \bi \hg case. In Section~2, the underlying theoretical background based on a \bi extension of GFF is discussed, which leads to the \bi counterpart of the algorithm $q$Telescope. In Section~3, the degree setting for solving the \bi key equation is established. Applications are given in Section~4 to illustrate the usage of the newly developed Mathematica implementation which is available by email request to \texttt{Axel.Riese@risc.uni-linz.ac.at}. \section{Theoretical Background} % ---------------------- In this section, \qnb greatest factorial factorization ($q$GFF) of polynomials, which has been introduced by Paule (cf.~Paule and Riese~\cite{PauleRiese:qZeil}) providing an algebraic explanation of \qnb \hg telescoping, is extended to the \bi \hg case. It turns out that to this end the \qnb case argumentation can be carried over almost word by word. \subsection{Bibasic Greatest Factorial Factorization} % ---------------------------------------- Let $\Z$ denote the set of all integers, and $\N$ the set of all non-negative integers. Let $p$, $q$, $x$, and $y$ be fixed indeterminates. Assume $K = L(\kappa_1, \dots, \kappa_n)$ to be the field of rational functions in a fixed number of indeterminates $\kappa_1, \dots, \kappa_n$, $n \in \N$, where $p \neq \kappa_i \neq y$ and $q \neq \kappa_i \neq x$, $1 \le i \le n$, over some computable field $L$ of characteristic $0$ and not containing $p$, $q$, $x$, and $y$. (For the sake of simplicity with regard to the implementation we will restrict ourselves to the case where $L$ is the rational number field $\mathbb{Q}$.) The transcendental extension of $K$ by the indeterminates $p$ and $q$ is denoted by $F$, i.e., $F = K(p,q)$. For $P \in \Fp$, let the \textit{\bi shift operator} $\epsilon$ be given by $(\epsilon P)(x,y) = P(q x,p y)$. The extension of this shift operator to the rational function field $\Fr$, the quotient field of the polynomial ring $\Fp$, will be also denoted by $\epsilon$. \begin{dfn} A polynomial $P \in F[x]$ (resp.\ $P \in F[y]$) is called \textit{\qnb monic} (resp.\ \textit{\pnb monic}) if $P(0) = 1$. A polynomial $P \in \Fp$ is called \textit{\bi monic} if $P(x,0) \ne 0 \ne P(0,y)$ and either $P(0,0) = 1$, or $P(0,0) = 0$ and the coefficients of $P$ are relatively prime polynomials in $F$.% \footnote[2]{In other words, $P$ is assumed to be primitive over $L[\kappa_1, \dots, \kappa_n, p, q]$ in this case, which will guarantee the uniqueness of the so-called \bi monic decomposition of a polynomial as shown below.} \end{dfn} \begin{exa} (i) The following polynomials are bibasic monic: \begin{align*} P_1(x,y) & = 1, & P_2(x,y) & = 1 - a p q x^2 y^3, & P_3(x,y) & = (1-q)^2 x^2 + p y.\\ \intertext{(ii) The following polynomials are not bibasic monic:} P_4(x,y) & = q, & P_5(x,y) & = x y - a p q x^2 y^3, & P_6(x,y) & = (1-q)^{-1} p x^2 + p y.\tag*{\qedsymbol} \end{align*} \end{exa} The properties of being \qnb monic, \pnb monic, and \bi monic are clearly invariant with respect to the \bi shift operator $\epsilon$, i.e., if $P$ is \qnb monic, \pnb monic, or \bi monic, then the same holds true for $\epsilon P$. Furthermore, the product of two \bi monic polynomials is again \bi monic. Also note that a \bi monic polynomial $P$ satisfies $\gcd(x,P) = 1 = \gcd(y,P)$. Evidently, any non-zero polynomial $P \in \Fp$ has a unique factorization, the \textit{\bi monic decomposition}, in the form \[ P = z \cdot x^\alpha \cdot y^\beta \cdot \bmp{P}, \] where $z \in F$, $\alpha, \beta \in \N$, and $\bmp{P} \in \Fp$ is \bi monic. The bibasic monic decomposition of a polynomial $P \ne 0$ can be computed easily as follows. Define $\alpha := \max\{i \in \N : x^i | P\}$, $\beta := \max\{j \in \N : y^j | P\}$, and put $\bfbar{P} := x^{-\alpha} \cdot y^{-\beta} \cdot P$. If $\bfbar{P}(0,0) \ne 0$ define $z := \bfbar{P}(0,0)$, otherwise let $l$ denote the least common multiple of all coefficient-denominators of $\bfbar{P}$, let $g$ denote the greatest common divisor of all coefficients of $l \cdot \bfbar{P}$, and define $z := g/l$. Then, for $\bmp{P} := z^{-1} \cdot \bfbar{P}$, the \bi monic decomposition of $P$ is given by $P = z \cdot x^\alpha \cdot y^\beta \cdot \bmp{P}$. \begin{exa} The bibasic monic decompositions of the polynomials $P_4$, $P_5$, and $P_6$ from the example above are given by \begin{equation*} P_4 = q \cdot x^0 \cdot y^0 \cdot 1, \quad P_5 = 1 \cdot x \cdot y \cdot (1 - a p q x y^2), \quad P_6 = \frac{p}{1-q} \cdot x^0 \cdot y^0 \cdot (x^2 + (1-q) y). \tag*{\qedsymbol} \end{equation*} \end{exa} Moreover, we assume the result of any $\gcd$ computation over $\Fp$ as being normalized in the following sense. If $P_1 = z_1 \cdot x^{\alpha_1} \cdot y^{\beta_1} \cdot \bmp{P}_1$ and $P_2 = z_2 \cdot x^{\alpha_2} \cdot y^{\beta_2} \cdot \bmp{P}_2$ are the \bi monic decompositions of $P_1, P_2 \in \Fp$, we define \[ \gcdpq(P_1, P_2) := \gcd(x^{\alpha_1}, x^{\alpha_2}) \cdot \gcd(y^{\beta_1}, y^{\beta_2}) \cdot\gcdpq(\bmp{P}_1, \bmp{P}_2), \] where the $\gcdpq$ of two \bi monic polynomials is understood to be \bi monic. The polynomial degree in $x$ and $y$ of any $P \in \Fp$ is denoted by $\deg_x(P)$ and $\deg_y(P)$, respectively. \begin{dfn} For any \bi monic polynomial $P \in \Fp$ and $k \in \N$, the \textit{$k$-th falling \bi factorial} $\qfac{P}{k}$ of $P$ is defined as \[ \qfac{P}{k} := \prod_{i=0}^{k-1} \epsilon^{-i} P. \] \end{dfn} Note that by the null convention $\prod_{i \in \emptyset} P_i := 1$ we have $\qfac{P}{0} = 1$. In general, polynomials arising in \bi \hg summation have several different representations in terms of falling \bi factorials. {}From all possibilities, we shall consider only the one taking care of maximal chains, which informally can be obtained as follows. One selects irreducible factors of $P$ in such a way that their product, say \[ P_{k,1}(x,y) \cdot P_{k,1}(q^{-1}x, p^{-1}y) \cdots P_{k,1}(q^{-k+1} x, p^{-k+1} y), \] forms a falling \bi factorial $\qfac{P_{k,1}}{k}$ of maximal length $k$. For the remaining irreducible factors of $P$ this procedure is applied again in order to find all $k$-th falling factorial divisors $\qfac{P_{k,1}}{k}, \dots, \qfac{P_{k,l}}{k}$ of that type. Then $\qfac{P_k}{k} := \qfac{P_{k,1} \cdots P_{k,l}}{k}$ forms the \bi factorial factor of $P$ of maximal length $k$. Iterating this procedure one gets a factorization of $P$ in terms of ``greatest" factorial factors. \begin{dfn} We say that $\langle P_1, \dots, P_k \rangle$, $P_i \in \Fp$, is a \textit{\bi GFF-form} of a \bi monic polynomial $P \in \Fp$, written as $\GFFpq (P) = \langle P_1, \dots, P_k \rangle$, if the following conditions hold:\\[2pt] \hspace*{.3cm}% ($\GFFpq\,$1) $P = \qfac{P_1}{1} \cdots \qfac{P_k}{k}$,\\[2pt] \hspace*{.3cm}% ($\GFFpq\,$2) each $P_i$ is \bi monic, and $k > 0$ implies $P_k \ne 1$,\\ \hspace*{.3cm}% ($\GFFpq\,$3) for $i \le j$ we have $\gcdpq(\qfac{P_i}{i},\epsilon P_j) = 1 = \gcdpq(\qfac{P_i}{i}, \epsilon^{-j} P_j)$. \end{dfn} Note that $\GFFpq (1) = \langle \rangle$. Condition ($\GFFpq\,$3) intuitively can be understood as prohibiting ``overlaps" of \bi factorials that violate length maximality. The following theorem states that, as in the \qnb \hg case, the \bi GFF-form is unique and thus provides a canonical form. \begin{thm} If $\langle P_1, \dots, P_k \rangle$ and $\langle P'_1, \dots, P'_l \rangle$ are \bi GFF-forms of a \bi monic polynomial $P \in \Fp$, then $k = l$ and $P_i = P'_i$ for all $1 \le i \le k$. \end{thm} \noindent \textit{Proof.} The corresponding result for the ordinary \hg case ($p = q = 1$) has been proved by Paule~\cite[Thm.~2.1]{Paule:GFF}. The arguments used there extend immediately to the \bi \hg case proceeding by induction on $d := \deg_x(P) + \deg_y(P)$. \qed \medskip {}From algorithmic point of view it is important to note that the \bi GFF-form can be computed in an iterative manner essentially involving only $\gcd$ computations. In \qnb \hg summation, the normalized $\gcd$ of a polynomial $P$ and its \qnb shift $\epsilon P$ plays a fundamental role, as the $\gcd$ of $P$ and its shift $E P$ does in ordinary \hg summation, where $(E P) (x) = P(x+1)$. The same is true for \bi \hg summation with respect to the \bi shift operator $\epsilon$. The mathematical and algorithmic essence lies in the following lemma. \begin{lem}[Fundamental $\boldsymbol{\GFFpq}$ Lemma] Let $P \in \Fp$ be a \bi monic polynomial with $\GFFpq (P) = \langle P_1, \dots, P_k \rangle$. Then \[ \gcdpq(P, \epsilon P) = \qfac{P_1}{0} \cdots \qfac{P_k}{k-1}. \] \end{lem} \noindent \textit{Proof.} Due to the choice of the \bi shift operator $\epsilon$, the proof of the so-called Fundamental $q$GFF Lemma (cf.~Paule and Riese~\cite[Lemma~1]{PauleRiese:qZeil}) can be carried over to the \bi \hg case completely unchanged. \qed \medskip Thus, if $\GFFpq (P) = \langle P_1, \dots, P_k \rangle$, then $\GFFpq (\gcdpq(P, \epsilon P)) = \langle P_2, \dots, P_k \rangle$. Consequently, dividing $P$ with $\GFFpq (P) = \langle P_1, \dots, P_k \rangle$ by $\epsilon^{-1} \gcdpq(P,\epsilon P)$ or $\gcdpq(P,\epsilon P)$ results in separating the product of the first, respectively last, falling \bi factorial entries, or in other words \[ \frac{P}{\epsilon^{-1} \gcdpq(P, \epsilon P)} = P_1 \cdot P_2 \cdots P_k \quad \text{and} \quad \frac{P}{\gcdpq(P,\epsilon P)} = P_1 \cdot (\epsilon^{-1} P_2) \cdots (\epsilon^{-k+1} P_k). \] \subsection{Bibasic Hypergeometric Telescoping} % ---------------------------------- A sequence $\fk$ is said to be \textit{\bi hypergeometric} (see, e.g.\ Petkov\v{s}ek, Wilf, and Zeilberger~\cite{PetkovsekWilfZeilberger:A=B}) in $p$ and $q$ over $F$, if there exists a rational function $\rho \in \Fr$ such that $f_{k+1}/f_{k} = \rho(q^k,p^k)$ for all $k$ where the quotient is well-defined. Assume we are given a \bi \hg sequence $\fk$. Then the problem of \bi \hg telescoping is to decide whether there exists a \bi \hg sequence $\gk$ such that \begin{equation} \label{hg1} g_{k+1} - g_k = f_k, \end{equation} and if so, to determine $\gk$ with the motive that for $a, b \in \Z$, $a \le b$, \[ \sum_{k=a}^b f_k = g_{b+1} - g_a, \] which solves the indefinite summation problem. For the rational function $\rho$, related to $f_{k+1}/f_k$ as above, there exists a representation $\rho(x,y) = z \cdot x^\alpha \cdot y^\beta \cdot \bmp{A}(x,y)/\bmp{B}(x,y)$ with \bi monic $\bmp{A}, \bmp{B} \in \Fp$, $z \in F$, and $\alpha, \beta \in \Z$, which we call a \textit{rational representation} of the \bi \hg sequence $\fk$. If additionally $\bmp{A}$ and $\bmp{B}$ are relatively prime, then $\rho(x,y)$ is called the \textit{reduced rational representation} of $\fk$. For $\alpha \in \Z$, let $\alpha_+ := \max(\alpha,0)$ and $\alpha_- := \max(-\alpha,0)$. It will be shown below that \bi \hg telescoping can be decided constructively as follows. \begin{alg}[\Telpq] \textsc{Input:} a \bi \hg sequence $\fk$ specified by its reduced rational representation $\rho = z \cdot x^\alpha \cdot y^\beta \cdot \bmp{A}/\bmp{B}$;\\ \textsc{Output:} a \bi \hg solution $\gk$ of \eqref{hg1}; in case such a solution does not exist, the algorithm stops. \begin{enumerate} \item Compute the \bi GP form of $\fk$, i.e., \begin{enumerate} \item determine unique \bi monic polynomials $\bmp{P}$, $\bmp{Q}$, $\bmp{R} \in \Fp$ such that \begin{equation} \label{gpmonic} \frac{\bmp{A}}{\bmp{B}} = \frac{\epsilon \bmp{P}}{\bmp{P}} \cdot \frac{\bmp{Q}}{\epsilon \bmp{R}}, \end{equation} where $\gcdpq(\bmp{P},\bmp{Q}) = 1 = \gcdpq(\bmp{P},\bmp{R})$ and $\gcdpq(\bmp{Q},\epsilon^j \bmp{R}) = 1$ for all $j \ge 1$, and \item let $a_x$, $b_x$, $a_y$, and $b_y$ denote the coefficients of the lowest occurring powers of $x$ and $y$ in $\bmp{A}(x,0)$, $\bmp{B}(x,0)$, $\bmp{A}(0,y)$, and $\bmp{B}(0,y)$, respectively. Define \[ (\gamma, \delta) := \begin{cases} (\varphi, \psi) & \begin{array}{l} \text{if $\alpha = 0 = \beta$ and $q^\varphi \cdot p^\mu \cdot b_y/a_y = z = p^\psi \cdot q^\nu \cdot b_x/a_x$}\\[-2pt] \text{for $\varphi, \psi \in \N$ and $\mu, \nu \in \Z$}, \end{array}\\ (\varphi, 0) & \text{if $\alpha = 0 \ne \beta$ and $z = q^\varphi \cdot p^\mu \cdot b_y/a_y$ for $\varphi \in \N$, $\mu \in \Z$},\\ (0, \psi) & \text{if $\alpha \ne 0 = \beta$ and $z = p^\psi \cdot q^\nu \cdot b_x/a_x$ for $\psi \in \N$, $\nu \in \Z$},\\ (0, 0) & \text{otherwise}, \end{cases} \] and put \begin{align} P & := x^\gamma \cdot y^\delta \cdot \bmp{P}, \notag\\ Q & := z \cdot q^{-\gamma} \cdot p^{-\delta} \cdot x^\alphap \cdot y^\betap \cdot \bmp{Q}, \label{PQRGP} \\ \epsilon R & := x^\alpham \cdot y^\betam \cdot \epsilon \bmp{R}, \notag \end{align} with the motive that then \[ \rho = \frac{\epsilon P}{P} \cdot \frac{Q}{\epsilon R}. \] \end{enumerate} \item Try to solve the \bi key equation \begin{equation} \label{gospereq} P = Q \cdot \epsilon Y - R \cdot Y \end{equation} for a polynomial $Y \in \Fp$. \item If such a polynomial solution $Y$ exists, then \begin{equation} \label{gkfk} g_k = \frac{R(q^k,p^k) \cdot Y(q^k,p^k)}{P(q^k,p^k)} \cdot f_k \end{equation} is a \bi \hg solution of \eqref{hg1}, otherwise no \bi \hg solution $\gk$ exists. \qed \end{enumerate} \end{alg} \pagebreak % !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! %\noindent The steps of Algorithm \Telpq\ are derived as follows. First, assume that a \bi \hg solution $\gk$ with rational representation $g_{k+1}/g_k = \sigma(q^k,p^k)$ of \eqref{hg1} exists. Then evidently we have \begin{equation} \label{hg2} g_k = \tau (q^k,p^k) \cdot f_k, \end{equation} where $\tau(x,y) = 1/(\sigma(x,y)-1) \in \Fr$. By relation \eqref{hg2}, equation~\eqref{hg1} is equivalent to \begin{equation} \label{hg3} z \cdot x^\alphap \cdot y^\betap \cdot \bmp{A} \cdot \epsilon \tau - x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot \tau = x^\alpham \cdot y^\betam \cdot \bmp{B}, \end{equation} where the reduced rational representation of $\fk$ is given by $\rho = z \cdot x^\alpha \cdot y^\beta \cdot \bmp{A}/\bmp{B}$. Vice versa, any rational solution $\tau \in \Fr$ of \eqref{hg3} gives rise to a \bi \hg solution $g_k := \tau(q^k,p^k) \cdot f_k$ of \eqref{hg1}. This means, \bi \hg telescoping is equivalent to finding a \textit{rational} solution $\tau$ of \eqref{hg3}. Any $\tau \in \Fr$ can be represented as the quotient of relatively prime polynomials in the form $\tau = \Uc/\Vc$ where $\Uc, \Vc \in \Fp$ with $\Vc = x^\varphi \cdot y^\psi \cdot \bmp{\Vc}$ the \bi monic decomposition of $\Vc$. In case such a solution $\tau$ of \eqref{hg3} exists, assume we know $\Vc$ or a multiple $V \in \Fp$ of $\Vc$. Then by clearing denominators in \[ z \cdot x^\alphap \cdot y^\betap \cdot \bmp{A} \cdot \frac{\epsilon U}{\epsilon V} - x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot \frac{U}{V} = x^\alpham \cdot y^\betam \cdot \bmp{B}, \] the problem reduces further to finding a polynomial solution $U \in \Fp$ of the resulting difference equation with polynomial coefficients, \begin{equation} \label{cap KEY} z \cdot x^\alphap \cdot y^\betap \cdot \bmp{A} \cdot V \cdot \epsilon U - x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot (\epsilon V) \cdot U = x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot V \cdot \epsilon V. \end{equation} Note that at least one polynomial solution, namely $U = \Uc \cdot V/\Vc$, exists. Furthermore, equations of that type simplify by canceling $\gcdpq$'s. For instance, in order to get more information about the denominator $\Vc$, let $\Vc_i := \epsilon^i \Vc / \gcdpq(\Vc,\epsilon \Vc)$, $i \in \{0,1\}$. Then \eqref{hg3} is equivalent to \begin{equation} \label{KEY} z \cdot x^\alphap \cdot y^\betap \cdot \bmp{A} \cdot \Vc_0 \cdot \epsilon \Uc - x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot \Vc_1 \cdot \Uc = x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot \Vc_0 \cdot \Vc_1 \cdot \gcdpq(\Vc,\epsilon \Vc). \end{equation} Now, if $\langle \Pc_1, \dots, \Pc_m \rangle$, $m \in \N$, is the \bi GFF-form of $\bmp{\Vc}$, it follows from $\gcdpq(\Uc,\Vc) = 1 = \gcdpq(\Vc_0,\Vc_1)$ and the Fundamental $\GFFpq$ Lemma that \[ \Vc_0 = (\epsilon^0 \Pc_1) \cdots (\epsilon^{-m+1} \Pc_m) \mid \bmp{B} \quad \text{and} \quad \Vc_1 = q^\varphi \cdot p^\psi \cdot (\epsilon \Pc_1) \cdots (\epsilon \Pc_m) \mid \bmp{A}. \] This observation gives rise to a simple and straightforward algorithm for computing a multiple $\bmp{V} := \qfac{P_1}{1} \cdots \qfac{P_n}{n}$ of $\bmp{\Vc}$. For instance, if $P_1 := \gcdpq(\epsilon^{-1} \bmp{A}, \bmp{B})$ then obviously $\Pc_1 | P_1$. Actually, one can iteratively extract \bi monic $\Pc_i$-multiples $P_i$ such that $\epsilon P_i | \bmp{A}$ and $\epsilon^{-i+1} P_i | \bmp{B}$ by the following algorithm. \begin{alg}[\VMULT] \textsc{Input:} relatively prime and \bi monic polynomials $\bmp{A}, \bmp{B} \in \Fp$ that constitute the \bi monic quotient of $\rho = z \cdot x^\alpha \cdot y^\beta \cdot \bmp{A}/\bmp{B} \in \Fr$;\\ \textsc{Output:} \bi monic polynomials $P_1, \dots, P_n$ such that $\bmp{V} := \qfac{P_1}{1} \cdots \qfac{P_n}{n}$ is a multiple of $\bmp{\Vc}$, the \bi monic part of the denominator $\Vc = x^\varphi \cdot y^\psi \cdot \bmp{\Vc}$ of $\tau \in \Fr$. \vspace{-3pt} \begin{enumerate} \item Compute $n = \min\{j \in \N \mid \gcdpq(\epsilon^{-1} \bmp{A}, \epsilon^{k-1} \bmp{B}) = 1$ for all integers $k > j \}$. \vspace{-5pt} \item Set $A_0 = \bmp{A}$, $B_0 = \bmp{B}$, and compute for $i$ from $1$ to $n$: \begin{align*} P_i & = \gcdpq(\epsilon^{-1} A_{i-1}, \epsilon^{i-1} B_{i-1}), \\ A_i & = A_{i-1}/\epsilon P_i,\\ B_i & = B_{i-1}/\epsilon^{-i+1} P_i. \tag*{\qed} \end{align*} \end{enumerate} \end{alg} A proof for the fact that the $P_i$ are indeed multiples of the $\Pc_i$ has been worked out for the ordinary \hg case by Paule~\cite[Lemma~5.1]{Paule:GFF}. It can be carried over to the \bi \hg world almost word by word. Hence we leave the steps of the verification to the reader. Note that in general step (i) of Algorithm \VMULT\ would be a rather time-consuming task involving resultant computations which could be solved by generalizing the univariate case (cf.~Abramov, Paule, and Petkov\v{s}ek~\cite{AbramovPaulePetkovsek:qDiffEq}) in a straightforward way, for instance, as follows. Define $R_1(v,w) := \res_x(\bmp{A}(x,y), \bmp{B}(v x, w y))$ and $R_2(v,w) := \res_y(\bmp{A}(x,y), \bmp{B}(v x, w y))$, viewed as polynomials of $v$ and $w$ over $F[y]$, respectively $F[x]$. Then $n$ is the maximal positive integer such that $R_1(q^n, p^n) \cdot R_2(q^n, p^n) = 0$ if such an integer exists, and $n = 0$ otherwise. However, in our implementation we make use of the fact that $\bmp{A}$ and $\bmp{B}$ already come in nicely factored form so that the computation of $n$ boils down to a comparison of those factors. Moreover, Algorithm \VMULT\ also delivers the constituents of the \bi monic part of the GP representation \eqref{gpmonic} as stated in the following lemma. \begin{lem} Let $n$, $A_n$, $B_n$, and the tuple $\langle P_1, \dots, P_n \rangle$ be computed as in Algorithm \textup{\VMULT}. Then for $\bmp{P} = \bmp{V}$, $\bmp{Q} = A_n$, and $\bmp{R} = \epsilon^{-1} B_n$ we have \[ \frac{\bmp{A}}{\bmp{B}} = \frac{\epsilon \bmp{P}}{\bmp{P}} \cdot \frac{\bmp{Q}} {\epsilon \bmp{R}}, \] where $\gcdpq(\bmp{P},\bmp{Q}) = 1 = \gcdpq(\bmp{P},\bmp{R})$ and $\gcdpq(\bmp{Q},\epsilon^j \bmp{R}) = 1$ for all $j \ge 1$. \end{lem} For more details on GP representations in the \qnb \hg case, see Abramov, Paule, and Petkov\v{s}ek~\cite{AbramovPaulePetkovsek:qDiffEq}, or Paule and Strehl~\cite{PauleStrehl:SymbSumm}. The results obtained there also apply in the \bi \hg case. With the multiple $\bmp{V}$ of $\bmp{\Vc}$ in hands, all what is left for solving \eqref{hg3}, and thus the \bi \hg telescoping problem \eqref{hg1}, is to determine appropriate multiplicities $\gamma$ and $\delta$ such that \[ V = x^{\gamma} \cdot y^{\delta} \cdot \bmp{V} \, \text{ is a multiple of } \, \Vc = x^\varphi \cdot y^\psi \cdot \bmp{\Vc}. \] For that we consider equation~\eqref{KEY} again in the equivalent version \begin{equation} \label{KEY2} z \cdot x^\alphap \cdot y^\betap \cdot \bmp{A} \cdot \bmp{\Vc} \cdot \epsilon \Uc - x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot q^\varphi \cdot p^\psi \cdot (\epsilon \bmp{\Vc}) \cdot \Uc = x^\alpham \cdot y^\betam \cdot \bmp{B} \cdot \bmp{\Vc} \cdot \epsilon \Vc, \end{equation} and distinguish the following cases corresponding to step (ib) of Algorithm \Telpq. \begin{enumerate} \item Assume that either $\alpha_- \ne 0$ or $\alpha_+ \ne 0$. In the first case we have $\alpha_+ = 0$ and $x^\alpham \,|\: \Uc$, hence $\varphi$ must be 0 because of $\gcdpq(\Uc,\Vc) = 1$. This means, we can choose $\gamma := 0$. In the second case we have $\alpha_- = 0$ and $x^{\min(\alpha_+,\varphi)} \,|\: \Uc$, because of $\epsilon \Vc = x^\varphi \cdot y^\psi \cdot q^\varphi \cdot p^\psi \cdot \epsilon \bmp{\Vc}$. Again $\varphi$ must be 0, and again we can choose $\gamma := 0$. Analogously, if $\beta \ne 0$ we can choose $\delta := 0$. \item Assume that $\alpha = 0$ and $\beta \ne 0$, hence $\psi = 0$ by (i). For $\varphi > 0$, evaluating equation~\eqref{KEY2} at $x = 0$ results in \begin{equation} \label{gamma1} z \cdot y^\betap \cdot \bmp{A}(0,y) \cdot \bmp{\Vc}(0,y) \cdot \Uc(0,py) - y^\betam \cdot \bmp{B}(0,y) \cdot q^\varphi \cdot \bmp{\Vc}(0,py) \cdot \Uc(0,y) = 0. \end{equation} In order to evaluate \eqref{gamma1} at $y = 0$, note that $P \in \Fp$ being \bi monic does not necessarily imply that $P(0,y) \in F[y]$ is \pnb monic. To overcome this problem, let us consider the \pnb monic decompositions of $\Uc(0,y)$ and $\bmp{\Vc}(0,y)$, say $\Uc(0,y) = u \cdot y^{\beta_u} \cdot \bfbar{U}(y)$ and $\bmp{\Vc}(0,y) = v \cdot y^{\beta_v} \cdot \bfbar{V}(y)$, respectively. Now, dividing equation~\eqref{gamma1} by $\Uc(0,y) \cdot \bmp{\Vc}(0,y) \ne 0$ leads to \begin{equation} \label{gamma2} z \cdot y^\betap \cdot \bmp{A}(0,y) \cdot p^{\beta_u} \cdot \frac{\bfbar{U}(py)}{\bfbar{U}(y)} - y^\betam \cdot \bmp{B}(0,y) \cdot q^\varphi \cdot p^{\beta_v} \cdot \frac{\bfbar{V}(py)}{\bfbar{V}(y)} = 0. \end{equation} \sloppy Additionally, let the \pnb monic decompositions of $\bmp{A}(0,y)$ and $\bmp{B}(0,y)$ be given by $\bmp{A}(0,y) = a_y \cdot y^{\beta_a} \cdot \bfbar{A}(y)$ and $\bmp{B}(0,y) = b_y \cdot y^{\beta_b} \cdot \bfbar{B}(y)$, respectively. Then the powers $y^{\beta_a+\betap}$ and $y^{\beta_b+\betam}$ must be equal, and after cancellation equation~\eqref{gamma2} at $y=0$ turns into \[ z \cdot a_y \cdot p^{\beta_u} - b_y \cdot q^\varphi \cdot p^{\beta_v} = 0. \] This means, we obtain as a condition for $\varphi > 0$ that $z = q^\varphi \cdot p^\mu \cdot b_y/a_y$ with $\mu \in \Z$. Hence, in this case we choose $\gamma := \varphi$, i.e., we set $\gamma$ to this \qnb power if $z$ has this particular form, and $\gamma := 0$ otherwise. Analogously, if $\alpha \ne 0$ and $\beta = 0$ we define $\delta := \psi > 0$, if $z = p^\psi \cdot q^\nu \cdot b_x/a_x$ with $\nu \in \Z$, and $\delta := 0$ otherwise. \item Finally, for the case $\alpha = 0 = \beta$ similar reasoning as in case (ii) leads to the conditions \begin{equation} \label{gamma3} q^\varphi \cdot p^\mu \cdot b_y/a_y = z = p^\psi \cdot q^\nu \cdot b_x/a_x, \end{equation} for $\varphi > 0$ or $\psi > 0$, and $\mu, \nu \in \Z$. Thus, if both conditions~\eqref{gamma3} are satisfied we choose $\gamma := \varphi$ and $\delta := \psi$, and otherwise $\gamma = \delta := 0$. \end{enumerate} \noindent The remaining steps of Algorithm \Telpq\ now are explained as follows. Once again, employing the GP representation for the \bi monic quotient of $\rho$, \[ \frac{\bmp{A}}{\bmp{B}} = \frac{\epsilon \bmp{P}}{\bmp{P}} \cdot \frac{\bmp{Q}}{\epsilon \bmp{R}}, \] it is easily seen that equation~\eqref{cap KEY} can be written as \begin{equation} \label{bibGo1} z \cdot q^{-\gamma} \cdot p^{-\delta} \cdot x^\alphap \cdot y^\betap \cdot \frac{\bmp{Q}}{\epsilon \bmp{R}} \cdot \epsilon U - x^\alpham \cdot y^\betam \cdot U = x^{\gamma+\alpham} \cdot y^{\delta+\betam} \cdot \bmp{P}. \end{equation} Because of relative primeness of certain polynomials, we observe that $x^\alpham | \, U$, $y^\betam | \, U$, and $\epsilon \bmp{R} \mid \epsilon U$. Hence by defining $Y$ by the relation \[ U = x^\alpham \cdot y^\betam \cdot q^{-\alpham} \cdot p^{-\beta_-} \cdot \bmp{R} \cdot Y, \] the task to solve equation~\eqref{cap KEY} for $U$ reduces to solve \begin{equation} \label{bibGo2} z \cdot q^{-\gamma} \cdot p^{-\delta} \cdot x^\alphap \cdot y^\betap \cdot \bmp{Q} \cdot \epsilon Y - x^\alpham \cdot y^\betam \cdot q^{-\alpham} \cdot p^{-\betam} \cdot \bmp{R} \cdot Y = x^{\gamma} \cdot y^{\delta} \cdot \bmp{P} \end{equation} for $Y \in \Fp$. By definition~\eqref{PQRGP} of $P$, $Q$, and $R$, equation~\eqref{bibGo2} immediately turns into the \bi key equation~\eqref{gospereq}, \[ Q \cdot \epsilon Y - R \cdot Y = P. \] Finally, from $U/V = R \cdot Y/P$, again by definition~\eqref{PQRGP}, it follows directly that \[ g_k = \frac{R(q^k,p^k) \cdot Y(q^k,p^k)}{P(q^k,p^k)} \cdot f_k \] as in~\eqref{gkfk} actually is a solution of the \bi \hg telescoping problem~\eqref{hg1}. This completes the proof of the correctness of Algorithm \Telpq. \section{Degree Setting for Solving the Bibasic Key Equation} % --------------------------------------------------- To solve the \bi key equation \begin{equation} \label{bibeq} P = Q \cdot \epsilon Y - R \cdot Y \end{equation} we first have to determine degree bounds $d_1$ and $d_2$, say, for the solution polynomial $Y \in \Fp$ with respect to $x$ and $y$, respectively, as shown in Theorem~\ref{degset} below. Then we put \[ Y(x,y) := \sum_{i=0}^{d_1} \sum_{j=0}^{d_2} y_{i,j} \cdot x^i \cdot y^j \] with undetermined $y_{i,j}$ and solve \eqref{bibeq} for the $y_{i,j}$ by equating to zero all coefficients of $x^i y^j$ in the equation \[ P - Q \cdot \epsilon Y + R \cdot Y = 0, \] which corresponds to solving a system of linear equations. \begin{thm} \label{degset} Let $l_Q^x(y)$, $l_Q^y(x)$, $l_R^x(y)$, and $l_R^y(x)$ denote the leading coefficient polynomials of $Q$ and $R$ with respect to $x$ and $y$, respectively. Let $QR^+ := Q+R$ and $QR^- := Q-R$. Then bounds for $\deg_x(Y)$ and $\deg_y(Y)$ are given by: \begin{enumerate} \item[(A$_x$)] If $\deg_x(QR^+) \neq \deg_x(QR^-)$, then \[ \deg_x(Y) \le \max\{\deg_x(P) - \max\{\deg_x(QR^+), \deg_x(QR^-)\}, 0\}. \] \item[(A$_y$)] If $\deg_y(QR^+) \neq \deg_y(QR^-)$, then \[ \deg_y(Y) \le \max\{\deg_y(P) - \max\{\deg_y(QR^+), \deg_y(QR^-)\}, 0\}. \] \item[(B$_x$)] If $\deg_x(QR^+) = \deg_x(QR^-)$, then \begin{enumerate} \item[(B1$_x$)] if $\deg_x(Q) \neq \deg_x(R)$, then \[ \deg_x(Y) = \deg_x(P) - \deg_x(QR^+), \] \item[(B2$_x$)] if $\deg_x(Q) = \deg_x(R)$, then \begin{enumerate} \item[(B2a$_x$)] if $l_R^x(y)/l_Q^x(y)$ is of the form $p^\mu \cdot q^\nu \cdot r(y)$ with $\mu, \nu \in \N$, and $r(y)$ a rational function with $r(0)=1$, then \[ \deg_x(Y) \le \max\{\deg_x(P) - \deg_x(QR^+), \nu\}, \] \item[(B2b$_x$)] otherwise \[ \deg_x(Y) = \deg_x(P) - \deg_x(QR^+). \] \end{enumerate} \end{enumerate} \item[(B$_y$)] If $\deg_y(QR^+) = \deg_y(QR^-)$, then \begin{enumerate} \item[(B1$_y$)] if $\deg_y(Q) \neq \deg_y(R)$, then \[ \deg_y(Y) = \deg_y(P) - \deg_y(QR^+), \] \item[(B2$_y$)] if $\deg_y(Q) = \deg_y(R)$, then \begin{enumerate} \item[(B2a$_y$)] if $l_R^y(x)/l_Q^y(x)$ is of the form $p^\mu \cdot q^\nu \cdot r(x)$ with $\mu, \nu \in \N$, and $r(x)$ a rational function with $r(0)=1$, then \[ \deg_y(Y) \le \max\{\deg_y(P) - \deg_y(QR^+), \mu \}, \] \item[(B2b$_y$)] otherwise \[ \deg_y(Y) = \deg_y(P) - \deg_y(QR^+). \] \end{enumerate} \end{enumerate} \end{enumerate} \end{thm} \noindent {\it Proof.} We rewrite the key equation to obtain \begin{equation} \label{qgosperrew} 2\,P = QR^+ \cdot (\epsilon Y - Y) + QR^- \cdot (\epsilon Y + Y). \end{equation} Cases (A$_x$) and (A$_y$) follow immediately. Note that it might happen that \[ \deg_x(QR^+) > \deg_x(P) \text{ and } \deg_x(QR^-) = \deg_x(P), \] and simultaneously \[ \deg_y(QR^+) > \deg_y(P) \text{ and } \deg_y(QR^-) = \deg_y(P). \] In this case, setting $\deg_x(Y) = \deg_y(Y) = 0$ could yield a solution, since $\epsilon Y - Y = 0$ then. For Case (B1$_x$) let $a := \deg_x(Q)$, $c := \deg_x(Y)$, and let $l_Y^x(y)$ denote the leading coefficient polynomial of $Y$ with respect to $x$. Assume that $\deg_x(Q) > \deg_x(R)$. Then \eqref{qgosperrew} gives \begin{align} 2\,P(x,y) & = (l_Q^x(y) \, x^a + \dots) \cdot [(l_Y^x(py) \, q^c - l_Y^x(y)) \, x^c + \dots] \notag\\ & + \,(l_Q^x(y) \, x^a + \dots) \cdot [(l_Y^x(py) \, q^c + l_Y^x(y)) \, x^c + \dots] \notag\\ & = 2\:l_Q^x(y) \: l_Y^x(py) \: q^c \: x^{a+c} + \dots.\label{B1x} \end{align} Clearly, the coefficient of $x^{a+c}$ in \eqref{B1x} will never vanish. Therefore we have \[ \deg_x(Y) = \deg_x(P) - \deg_x(Q). \] Including the case $\deg_x(Q) < \deg_x(R)$, we obtain \[ \deg_x(Y) = \deg_x(P) - \max\{\deg_x(Q), \deg_x(R)\} = \deg_x(P) - \deg_x(QR^+). \] Analogous reasoning leads to Case (B1$_y$). For Case (B2$_x$) we similarly observe that \begin{align} 2\,P(x,y) & = [(l_Q^x(y) + l_R^x(y)) \, x^a +\dots] \cdot [(l_Y^x(py) \, q^c - l_Y^x(y)) \, x^c + \dots] \notag\\ & + \,[(l_Q^x(y) - l_R^x(y)) \, x^a + \dots] \cdot [(l_Y^x(py) \, q^c + l_Y^x(y)) \, x^c + \dots] \notag\\ & = 2\:[l_Q^x(y) \: l_Y^x(py) \: q^c - l_R^x(y) \: l_Y^x(y)] \: x^{a+c} + \dots. \label{B2x} \end{align} Now we no longer have the guarantee that the coefficient of $x^{a+c}$ in \eqref{B2x} does not vanish, but it is easily seen that this happens only for \begin{equation} \label{vancond} q^c = \frac{l_R^x(y)}{l_Q^x(y)} \cdot \frac{l_Y^x(y)}{l_Y^x(py)}. \end{equation} Note that $l_Y^x(y)$ is actually \textit{not} known. However, for any non-zero polynomial $h(y) = h_0 + h_1\,y + \dots + h_d\,y^d$, the quotient $h(y)/h(py)$ is of the form $p^{-m} \cdot s(y)$, where $s(y)$ is a rational function with $s(0)=1$ and $m$ is the zero-root multiplicity of $h(y)$. Hence, the rightmost fraction in \eqref{vancond} may eliminate only positive integer powers of $p$ and a rational function of $y$ but never introduce a power of $q$. This proves Case (B2a$_x$), and after interchanging $x$ and $p$ with $y$ and $q$, respectively, also Case (B2a$_y$). On the other hand, if the coefficient of $x^{a+c}$ in \eqref{B2x} does not vanish, we obtain Case (B2b$_x$) and analogously Case (B2b$_y$). \qed \section{Applications} % ------------ In this section we shall illustrate the method of \bi \hg telescoping using the author's Mathematica implementation \texttt{qTelescope}, which is a \bi extension of a \qnb analogue of Gosper's algorithm originally described in Paule and Riese~\cite{PauleRiese:qZeil}. Let the \textit{\qnb shifted factorial} of $a \in F$ be defined as usual (see, e.g.~Gasper and Rahman~\cite{GasperRahman:BHS}) by \begin{align*} &(a;q)_k := \begin{cases} (1 - a) (1 - a q) \cdots (1 - a q^{k-1}), & \text{if $k > 0$,} \\ 1, & \text{if $k = 0$,} \\ \big[ (1 - a q^{-1}) (1 - a q^{-2}) \cdots (1 - a q^k) \big]^{-1}, & \text{if $k < 0$,} \end{cases}\\ \intertext{and} &(a;q)_\infty := \prod_{k=0}^\infty (1 - a q^k),\\ \intertext{where products of \qnb shifted factorials will be abbreviated by} &(a_1, a_2, \dots, a_n;q)_k := (a_1;q)_k \, (a_2;q)_k \cdots (a_n;q)_k. \end{align*} In the present implementation we allow as summand any \bi \hg sequence $\fk$ of the form \begin{align*} f_k & = \frac{\prod_r (C_r \, q^{(c_r i_r) k + d_r};q^{i_r})_ {a_r k + b_r}}{\prod_s (D_s \, q^{(v_s j_s) k + w_s};q^{j_s})_ {t_s k + u_s}} \cdot \frac{\prod_r (C_r' \, p^{(c'_r i'_r) k + d'_r};p^{i'_r})_ {a'_r k + b'_r}}{\prod_s (D_s' \, p^{(v'_s j'_s) k + w'_s}; p^{j'_s})_{t'_s k + u'_s}}\\ & \quad \times R(q^k,p^k) \cdot q^{\alpha \binom{k}{2}} \cdot p^{\beta \binom{k}{2}} \cdot z^k, \end{align*} \begin{tabbing}with \hspace{1.45cm} \= \\ $C_r, D_s$ \> power products in $K(p)$,\\ $C_r', D_s'$ \> power products in $K(q)$,\\ $a_r, t_s, a'_r, t'_s$ \> specific integers (i.e., integers free of any parameters),\\ $b_r, u_s, b'_r, u'_s$ \> integer parameters free of $k$, or $\pm \infty$ if $a_r \, (\text{resp.\ } t_s, a'_r, t'_s) = 0$,\\ $c_r, v_s, c'_r, v'_s$ \> specific integers,\\ $d_r, w_s, d'_r, w'_s$ \> integer parameters free of $k$,\\ $i_r, j_s, i'_r, j'_s$ \> specific non-zero integers,\\ $R$ \> a rational function in $F(q^k,p^k)$ such that the denominator factors completely\\ \> into a product of terms of the form $(1 - D\,q^{v k+w})$ and $(1 - D' p^{v' k+w'})$,\\ $\alpha, \beta$ \> specific integers, and\\ $z$ \> a rational function in $F$. \end{tabbing} For the actual computation of the GP representation let $\rho(x,y)$ denote the possibly non-reduced rational representation of the summand $f_k$. It is obvious from the input specification that $\rho$ can always be converted into the form \begin{align*} \rho(x,y) & = \frac{(\epsilon \bfbar{P})(x,y)}{\bfbar{P}(x,y)} \cdot \frac{\prod_i (1 - \Gamma_i\,x^{\gamma_i})} {\prod_j (1 - \Delta_j\,x^{\delta_j})} \cdot \frac{\prod_i (1 - \Gamma_i'\,y^{\gamma_i'})} {\prod_j (1 - \Delta_j'\,y^{\delta_j'})} \cdot x^{\bfbar{\alpha}} \cdot y^{\bfbar{\beta}} \cdot \bfbar{z} \\ & = \frac{(\epsilon \bfbar{P})(x,y)}{\bfbar{P}(x,y)} \cdot \frac{\bfbar{A}(x,y)}{\bfbar{B}(x,y)} \cdot x^{\bfbar{\alpha}} \cdot y^{\bfbar{\beta}} \cdot \bfbar{z}, \end{align*} where $\bfbar{P} \in \Fp$ is \bi monic and satisfies $\gcdpq(\bfbar{P}, \bfbar{A}) = 1 = \gcdpq(\epsilon \bfbar{P}, \bfbar{B})$; the $\Gamma_i, \Delta_j, \Gamma_i', \Delta_j'$ are power products in $F$, the $\gamma_i, \delta_j, \gamma_i', \delta_j'$ are positive integers, $\bfbar{\alpha}, \bfbar{\beta} \in \Z$, and $\bfbar{z} \in F$. Concerning Algorithm \VMULT, it is clear from above that any $\bfbar{P} \ne 1$ will actually contribute to $\qfac{P_1}{1}$ and thus can be treated separately. Due to our input restrictions --- this is the reason for admitting only power products instead of arbitrary rational functions --- it is possible to find $n$ in step (i) of Algorithm \VMULT\ simply by comparing all factors in $\bfbar{A}$ and $\bfbar{B}$ as already mentioned. Furthermore, since $\bfbar{A}$ and $\bfbar{B}$ are both products of a \qnb monic and a \pnb monic polynomial, they will never contribute to $b_x/a_x$ and $b_y/a_y$. Thus, $b_x/a_x$ and $b_y/a_y$ are in any case integer powers of $q$ and $p$, respectively, coming from $\epsilon \bfbar{P}/\bfbar{P}$. Therefore, they do not take influence on the computation of $\gamma$ and $\delta$ at all. \subsection{Bibasic Summation Formulas} % -------------------------- In 1989, Gasper~\cite{Gasper:BiSer} derived the indefinite \bi summation formula \begin{align} \sum_{k=0}^n f_k & = \sum_{k=0}^n \frac{(1-a p^k q^k)(1-b p^k q^{-k})} {(1-a)(1-b)} \; \frac{(a,b;p)_k \, (c,a/b c;q)_k} {(q,a q/b;q)_k \, (a p/c,b c p;p)_k} \; q^k \notag\\ & = \frac{(a p,b p;p)_n \, (c q,a q/b c;q)_n} {(q,a q/b;q)_n \, (a p/c,b c p;p)_n} = g_n \label{bigasper} \end{align} by showing that $g_k$ is a \bi \hg solution of the equation $f_k = g_k - g_{k-1}$, however, without revealing how to come up with $g_k$. With our implementation the job of finding $g_k$ is left to the computer. \small\begin{verbatim} In[1]:= <