# BEAM–BEAM EFFECTS IN SPACE CHARGE DOMINATED ION BEAMS ††thanks: Work
supported under Contract Number DE-AC02-98CH10886 with the auspices of the US
Department of Energy

C. Montag  A. Fedotov  Brookhaven National Laboratory  Upton  NY  USA

###### Abstract

During low-energy operations below the regular injection energy in the
Relativistic Heavy Ion Collider (RHIC), significant beam lifetime reductions
due to the beam–beam interaction in conjunction with large space charge tune
shifts have been observed. We report on dedicated experiments aimed at
understanding this phenomenon as well as preliminary simulation results, and
propose alternative working points to improve the beam lifetime in future low-
energy RHIC runs.

##  1 INTRODUCTION

One of the major physics programmes at the RHIC for the next 5–10 years is the
search for the critical point in the Quantum ChromoDynamics (QCD) phase
diagram (Fig. 1), which is expected to occur at centre-of-mass energies in the
range of $\sqrt{s_{\rm NN}}$ = 5–30 GeV/n. This requires colliding gold beams
with energies between 2.5 and 15 GeV/nucleon, which is well below the nominal
energy range of 10–100 GeV/n in the RHIC (Fig. 2). In conjunction with the
circumference of 3.8 km, this low energy results in a significant direct space
charge tune shift up to $\Delta Q_{\rm sc}=0.1,$ which is more than ten times
larger than the total beam–beam parameter $\xi_{\rm beam-beam}$ encountered
during low-energy operation [1]. Experiments with a large beam–beam parameter
comparable with the space charge tune shift have been performed using protons
and are reported elsewhere [2]. Table 1 lists space charge and beam–beam
parameters achieved in the RHIC for different Au beam energies. However, in
spite of the large difference in magnitude of these two effects, and the
similar functional dependence of the associated forces on transverse particle
coordinates, we have observed a significant deterioration of beam lifetimes
once beams are brought into collision.

To gain a better understanding of this phenomenon, we have performed a series
of beam experiments and developed a simulation code. In the following
sections, we describe our experimental observations during regular low-energy
operations with Au ions, as well as those dedicated Au beam experiments.
Furthermore, we report on first results obtained from simulations.

Figure 1: The QCD phase diagram. A lower centre-of-mass energy $\sqrt{s_{\rm NN}}$ corresponds to a higher baryon chemical potential. The critical point is expected to be in the energy range between $\sqrt{s_{\rm NN}}=5$ and $30\,{\rm GeV}$. Figure 2: An aerial view of the RHIC accelerator complex, with its two $3.8\,{\rm km}$ circumference storage rings, ‘Blue’ and ‘Yellow’. Table 1: Beam Lifetimes $\tau$ with and without Collisions at Different Energies in the RHIC, with the Corresponding Space Charge Tune Shifts $\Delta Q_{\rm sc}$ and Beam–Beam Parameters $\xi_{\rm beam-beam}$ E [GeV/n] | $\Delta Q_{\rm sc}$ | $\xi_{\rm beam-beam}$ |  $\tau$ [s]
---|---|---|---
9.8 | 0.03 | 0 | 2000
9.8 | 0.03 | 0.002 | 600
5.75 | 0.05 | 0 | 1600
5.75 | 0.05 | 0.0015 | 400
5.75 | 0.09 | 0 | 700
5.75 | 0.09 | 0.0027 | 260
3.85 | 0.11 | 0 | 70
3.85 | 0.08 | 0.003 | 70

##  2 EXPERIMENTAL OBSERVATIONS

When the RHIC was operated at a beam energy of $E=3.85\,{\rm GeV/n},$ a tune
scan was performed to maximize the beam lifetime. Starting at the regular RHIC
heavy ion working point of $(Q_{x}/Q_{y})=(28.23/29.22),$ the tunes were
lowered and the beam lifetime in collision was observed (see Fig. 3). This
resulted in a new working point of $(Q_{x}/Q_{y})=(28.17/29.16)$; during the
course of the run this was further lowered to $(28.13/29.12).$ This latest
working point was subsequently used at $E=5.75\,{\rm GeV/n}$ as well.

Figure 3: Beam intensities (middle plot) in the Blue and Yellow RHIC rings
during a tune scan at $E=5.75\,{\rm GeV/n}$ beam energy. The Yellow tunes are
shown in the top part of the plot and the Blue tunes in the bottom third. The
best Blue beam lifetime is achieved at a working point of
$(Q_{x}/Q_{y})=(28.17/29.16).$

During the course of the run, a strong effect of beam–beam interactions on the
lifetime of the space charge dominated beams was consistently observed, as
illustrated in Figs. 4–6. Figure 4 shows the intensity of individual bunches
in the Yellow RHIC ring at a beam energy of $E=3.85\,{\rm GeV/n}.$ Although
the initial intensity drops rather quickly for the first couple of minutes
while beams are not colliding, there is a sudden, sharp decrease in beam
lifetime, to roughly the same level as at the beginning of store, as soon as
the two beams begin colliding. Since the intensity of the Yellow bunch at this
time is only about half the initial value, which reduces the space charge tune
shift by the same factor of 2, this lifetime deterioration cannot simply be
explained by the total tune shift; that is, the sum of the space charge and
beam–beam tune shift. Moreover, a significant beam–beam effect is observed for
bunches with a much smaller intensity, as can be seen in Fig. 4.

Figure 4: The intensities of individual bunches in the Yellow ring, at
$E=3.85\,{\rm GeV/n}$ beam energy. Collisions start at 10:14, resulting in a
sudden decrease in the lifetime.

As a second example, we discuss the evolution of the total intensity of the
Yellow beam at $E=5.75\,{\rm GeV/n}$ while the Blue ring is being filled (Fig.
5). In this case, there is no transverse separation of the two beams during
the injection process, so the injection of each individual Blue bunch results
in a Yellow bunch starting to experience beam–beam collisions. As a result,
the total Yellow beam lifetime slowly deteriorates the more bunches undergo
collisions with the newly injected Blue bunches.

Figure 5: The total beam intensity in the two RHIC rings during injection,
without any transverse separation of the two beams. After the Yellow ring has
been filled, Blue is being injected, resulting in a gradual decrease of the
Yellow beam lifetime due to the beam–beam interaction.

Finally, we focus on the beam decay rate at the end of a Au store with
$E=5.75\,{\rm GeV/n}$ beam energy (Fig. 6). When the beam–beam force on the
Blue beam disappears due to dumping of the oncoming Yellow beam, its decay
rate improves dramatically.

Figure 6: Beam decay rates during several Au beam stores at $5.75\,{\rm
GeV/n}$ beam energy. The Blue beam decay rate improves dramatically as soon as
the Yellow beam is dumped at the end of each store (see insert). Note that the
algorithm to calculate the beam decay rate from the measured beam intensity
has a time constant of 20 s. Hence, the actual drop in the instantaneous beam
decay is even more dramatic than suggested in this picture.

Based on operational experience and the desire to improve beam lifetimes and
therefore integrated luminosities in future RHIC low-energy runs, a dedicated
beam experiment aimed at searching for a better working point was performed.
Since the spacing of non-linear resonances is largest in the vicinity of the
integer resonance, fractional tunes below $0.1$ were proposed as the most
promising candidates. These studies were performed at the regular RHIC
injection energy of $9.8\,{\rm GeV/n},$ with a space charge tune shift of
$\Delta Q_{\rm sc}=0.03$ and a beam–beam tune-shift parameter of $\xi_{\rm
beam-beam}=0.002.$

As already observed during low-energy operations, the beam lifetime
deteriorated substantially when the two beams were brought into collision at
the regular RHIC working point of $(Q_{x}/Q_{y})=(28.23/29.22)$ (Fig. 7).

Figure 7: The beam intensities in the two RHIC rings during a beam experiment
at regular injection energy ($E=9.8\,{\rm GeV/n}$), at the regular RHIC heavy
ion working point of $(Q_{x}/Q_{y})=(28.23/29.22).$ The beams are brought into
collision shortly after injection, resulting in significant decrease in the
lifetime.

However, when the experiment was repeated at a near-integer working point of
$(Q_{x}/Q_{y})=(28.08/29.09)$ in Yellow and $(Q_{x}/Q_{y})=(28.08/29.07)$ in
Blue, there was no discernable effect on the Blue lifetime, while the Yellow
lifetime still deteriorated somewhat when beams were brought into collision,
as shown in Fig. 8.

Figure 8: The beam intensities in the RHIC at $E=9.8\,{\rm GeV/n},$ at a near-
integer working point of $(Q_{x}/Q_{y})=(28.08/29.09)$ in Yellow and
$(Q_{x}/Q_{y})=(28.08/29.07)$ in Blue. When beams are brought into collision,
the lifetime of the Yellow beam suffers, while the Blue beam is unaffected.

The cause of the differing behaviour in the two rings is not yet understood.
It may be attributable to parameters such as chromaticity, coupling control,
or the different working point above the diagonal, which may have been less
than optimal in the Yellow ring during the experiment. However, this result is
very encouraging for future low-energy operations, although the space charge
tune shift during this injection energy experiment was a factor of 2–3 smaller
than at the lower energies.

##  3 SIMULATIONS

To investigate the root cause of the lifetime deterioration, we performed
tracking simulations with a space charge tune shift of $\Delta Q_{\rm
sc}=0.06$ and a beam–beam parameter in each of the two RHIC interaction points
of $\xi_{\rm beam-beam}=0.003.$ Using these parameters, tune scans as well as
a frequency map analysis at a fixed working point were applied.

###  3.1 The Model

Space charge simulations are usually very CPU-time consuming because of
frequent recalculations of the particle distribution and the associated
electromagnetic fields. In the particular problem studied here, however, we
can take advantage of the fact that the evolution of the particle distribution
is comparatively slow. This is indicated by the beam lifetime of several
minutes to tens of minutes. Typical simulations track particles only over a
number of turns that corresponds to seconds of real time, so we can safely
assume that the distribution of our test particles does not change appreciably
over the course of the simulation. This approach, which is equivalent to the
weak–strong simulation technique applied in numerical beam–beam studies,
significantly speeds up the computation. In addition, no artificial noise is
introduced into the simulation by the finite number of particles, since
recalculation of the electromagnetic fields from the actual particle
distribution is avoided. Instead, we assume that the distribution remains
Gaussian during the entire simulation process. The r.m.s. width of this
Gaussian distribution is calculated from the beam emittance and the local
$\beta$ function, including the dynamic $\beta$-beat introduced by the space
charge and beam–beam forces around the machine. The accelerator model is based
on the RHIC lattice as described in MAD. So far, no lattice non-linearities
except the chromaticity correction sextupoles and the sextupole error in the
main dipoles have been included. Particles are tracked element by element, and
space charge kicks are applied at every quadrupole around the machine. Two
beam–beam interaction points are included in IPs 6 and 8. Synchrotron
oscillations are included, and the modulation of the space charge kick due to
the resulting longitudinal position oscillations is taken into account.

###  3.2 Results

To study the emittance growth as a function of tune, we launch 1000 particles
with a Gaussian distribution in all six phase space coordinates and track them
over 20 000 turns. At the end of each turn $i$, we calculate the 4-D
transverse emittance:

$\displaystyle\epsilon(i)=\left|\begin{array}[]{cccc}\langle xx\rangle&\langle
xx^{\prime}\rangle&\langle xy\rangle&\langle xy^{\prime}\rangle\\\ \langle
x^{\prime}x\rangle&\langle x^{\prime}x^{\prime}\rangle&\langle
x^{\prime}y\rangle&\langle x^{\prime}y^{\prime}\rangle\\\ \langle
yx\rangle&\langle yx^{\prime}\rangle&\langle yy\rangle&\langle
yy^{\prime}\rangle\\\ \langle y^{\prime}x\rangle&\langle
y^{\prime}x^{\prime}\rangle&\langle y^{\prime}y\rangle&\langle
y^{\prime}y^{\prime}\rangle\\\ \end{array}\right|^{\frac{1}{4}},$ (5)

where $\langle\cdots\rangle$ indicates the average over all particles. The
emittance growth rate

$\displaystyle\tau_{\epsilon}^{-1}=\frac{1}{\epsilon}\frac{{\rm
d}\epsilon}{{\rm d}t}$ (6)

is then computed as a function of tune by a linear fit to this 4-D emittance
evolution.

For the initial tune scan, depicted in Fig. 9, we varied the tunes in steps of
$\Delta Q_{x,y}=0.01,$ with $Q_{y}=Q_{x}-0.01.$ In the absence of the
beam–beam interaction, the 4-D emittance growth rate at fractional tunes below
0.2 is significantly lower than above, which qualitatively agrees with
experimental observations. Once beam–beam interactions are added, the
emittance growth rate increases over most of the tune range; however, the
resulting growth rate is, within error bars, independent of tune.

Figure 9: The 4-D emittance growth rate
$\tau_{\epsilon}^{-}1=\frac{1}{\epsilon}\frac{{\rm d}\epsilon}{{\rm d}t}$ as a
function of tune, with the working point $(Q_{x}/Q_{y})$ chosen such that
$Q_{y}=Q_{x}-0.01.$

To determine the tune footprint and the tune diffusion frequency, map analysis
[3, 4] was applied at a fixed working point of $(Q_{x}/Q_{y})=(28.13/29.12)$
for the two cases with and without beam–beam interaction. For this purpose, we
track a single test particle over $2^{14}$ turns and apply fast Fourier
transforms to calculate the horizontal and vertical tunes $(Q_{x,1}/Q_{y,1})$
and $(Q_{x,2}/Q_{y,2})$ for the first and second $2^{13}$ turns. To increase
the tune resolution, we apply an interpolation technique [5].

The tune diffusion is measured as

$\displaystyle|\Delta Q|=\sqrt{|Q_{x,1}-Q_{x,2}|^{2}+|Q_{y,1}-Q_{y,2}|^{2}}.$
(7)

The resulting tune footprint and tune diffusion is plotted in Fig. 10. While
the tune footprint overlaps the coupling resonance $Q_{x}=Q_{y}$ with as well
as without beam–beam interaction, the presence of the beam–beam force
significantly enhances the tune diffusion around that resonance.

Figure 10: Tune footprints at a nominal working point of
$(Q_{x}/Q_{y})=(28.13/29.12),$ without (top) and with (bottom) beam–beam
interaction.

Plotting the same data in the amplitude space (Fig. 11) reveals that this
enhanced tune diffusion occurs for amplitudes $(A_{x},A_{y})$ in the region
$\sigma_{x}<A_{x}<2\sigma_{x}$ and $2\sigma_{y}<A_{y}<4\sigma_{y}.$

Figure 11: Tune diffusion in the amplitude space, as obtained from the
frequency map analysis at a nominal working point of
$(Q_{x}/Q_{y})=(28.13/29.12),$ without (top) and with (bottom) beam–beam
interaction.

As shown above, the largest tune diffusion occurs around the coupling
resonance. This behaviour suggests that it might be beneficial to increase the
tune split between the two planes, thus selecting a working point further away
from the coupling resonance. To study this hypothesis, we performed a tune
scan with $Q_{y}=Q_{x}-0.02$ and determined the 4-D emittance growth rates. As
shown in Fig. 12, the effect of the beam–beam interaction on the emittance
growth rate is significantly reduced.

Figure 12: The 4-D emittance growth rate
$\tau_{\epsilon}^{-}1=\frac{1}{\epsilon}\frac{{\rm d}\epsilon}{{\rm d}t}$ as a
function of tune, with the working point $(Q_{x}/Q_{y})$ chosen at an
increased tune split with $Q_{y}=Q_{x}-0.02.$

Figure 13: Tune footprints at a nominal working point with an increased tune
split between the two planes, $(Q_{x}/Q_{y})=(28.14/29.12),$ without (top) and
with (bottom) beam–beam interaction.

Figure 14: Tune diffusion in the amplitude space, as obtained from the
frequency map analysis at a nominal working point with an increased tune split
between the two planes, $(Q_{x}/Q_{y})=(28.14/29.12),$ without (top) and with
(bottom) beam–beam interaction.

This observation is supported by the results of a frequency map analysis for
the working point $(Q_{x}/Q_{y})=(28.14/29.12)$ (Figs. 13 and 14), which shows
little effect of the beam–beam interaction on tune diffusion.

While this result is certainly very encouraging in terms of improving the
machine performance, it is worthwhile repeating that our model so far does not
include any magnet non-linearities beyond sextupoles, which may lead to
increased tune diffusion around the associated non-linear resonances in the
presence of the beam–beam interaction.

##  4 Summary

We have studied the effects of beam–beam interactions in colliding beams with
large direct space charge parameters up to $\Delta Q_{\rm sc}=0.1$ both
experimentally and through simulations. During RHIC low-energy operations as
well as dedicated experiments, we have consistently observed a strong effect
of the beam–beam interaction on the lifetime of the stored beam, although the
associated beam–beam parameter was about an order of magnitude smaller than
the space charge tune shift.

To provide the maximum tune space between non-linear resonances, we have
operated the RHIC at a near-integer working point. In this case, we observed
no discernable lifetime reduction in the Blue ring, while Yellow still
suffered. The root cause of this difference between the two rings is still
unknown, and may be related to parameters such as chromaticity and coupling
control, or the particular working point during the experiment. However, this
result is very encouraging for future low-energy runs, although the
corresponding space charge tune shift of $\Delta Q_{\rm sc}=0.03$ was
comparably modest.

Although a quantitative comparison of our simulation results with experimental
observations in the RHIC is difficult due to the lack of effects such as
intra-beam scattering in the simulation code, the tracking model presented
here qualitatively reproduces the main experimental result, namely the strong
effect of the beam–beam interaction in the presence of a large space charge
tune shift. Based on these simulation results, an alternative working point
further away from the coupling resonance $Q_{x}=Q_{y},$ which appears to be
the main source of emittance growth in our simulations, may be beneficial.
However, it is of utmost importance to repeat these simulations after higher-
order multipole errors have been added to the tracking model, which may lead
to a further reduction of the usable tune space.

## References

  * [1] A. Fedotov et al., “Beam Lifetime and Limitations during Low-Energy RHIC Operation,” Proc. PAC’11.
  * [2] A. Fedotov et al., “Interplay of Space Charge and Beam–Beam Effects in a Collider,” Proc. HB2010 (Morschach, Switzerland, 2010), p. 634.
  * [3] J. Laskar, Physica D 67 (1993) 257.
  * [4] J. Laskar and D. Robin, Part. Accel. 54 (1996) 183.
  * [5] A. Bazzani et al., Part. Accel. 147 (1996) 147.


# Hypergeometric Solutions of Linear Difference Systems
In Memory of Marko Petkovšek and Manual Bronstein

Moulay Barkatou
Université de Limoges, XLIM
123, Av. A. Thomas
87060 Limoges cedex, France
<EMAIL_ADDRESS>Mark van Hoeij
Department of Mathematics
Florida State University
Tallahassee, FL 32306, USA
<EMAIL_ADDRESS>Supported by NSF grant 2007959. Johannes Middeke
Research Institute for
Symbolic Computation (RISC)
Johannes Kepler University
4040 Linz, Austria
<EMAIL_ADDRESS>Yi Zhou
Department of Mathematics
Florida State University
Tallahassee, FL 32306, USA
<EMAIL_ADDRESS>
###### Abstract

We extend Petkovšek’s algorithm for computing hypergeometric solutions of
scalar difference equations to the case of difference systems $\tau(Y)=MY$,
with $M\in\mathop{\rm GL}_{n}({C}(x))$, where $\tau$ is the shift operator.
Hypergeometric solutions are solutions of the form $\gamma P$ where
$P\in{C}(x)^{n}$ and $\gamma$ is a hypergeometric term over ${C}(x)$, i.e.
${\tau(\gamma)}/{\gamma}\in{C}(x)$. Our contributions concern efficient
computation of a set of candidates for ${\tau(\gamma)}/{\gamma}$ which we
write as $\lambda=c\frac{A}{B}$ with monic $A,B\in{C}[x]$, $c\in{C}^{*}$.
Factors of the denominators of $M^{-1}$ and $M$ give candidates for $A$ and
$B$, while another algorithm is needed for $c$. We use super-reduction
algorithm [5] to compute candidates for $c$, as well as other ingredients to
reduce the list of candidates for $A/B$. To further reduce the number of
candidates $A/B$, we bound the type of $A/B$ by bounding local types. Our
algorithm has been implemented in Maple and experiments show that our
implementation can handle systems of high dimension, which is useful for
factoring operators.

##  1 Introduction

Let $C$ be a subfield of $\mathbb{C}$. Let $\tau$ be the shift operator, which
sends a function $y(x)$ to $y(x+1)$. A nonzero element $y$ in a universal
extension of $C(x)$ is called _hypergeometric_ if $\tau(y)/y\in C(x)$, in
other words, $y$ is a solution of a first order operator. If $y_{1},y_{2}$ are
hypergeometric then $y_{1}y_{2}$ is hypergeometric as well. Denote by
$\operatorname{hyp}(\lambda)$ a nonzero solution of $\tau-\lambda$. The
notation is defined up to a nonzero constant.

The first algorithm to compute hypergeometric solutions of an operator was
given by Petkovšek. One writes hypergeometric solutions of
$L=a_{n}\tau^{n}+\cdots+a_{0}\tau^{0}\in C[x][\tau]$ over $C(x)$ in this
format:

$y=\operatorname{hyp}(\lambda)P$

where $\lambda=c\frac{A}{B}$ with $A,B,P\in C[x]$, $A$ and $B$ are coprime and
monic, and $c\in C-\\{0\\}$. Such $y$ can be rewritten as
$\operatorname{hyp}(\lambda^{\prime})$ where
$\lambda^{\prime}=\lambda\frac{\tau(P)}{P}$. Allowing a polynomial factor $P$
in candidate solutions $y=\operatorname{hyp}(\lambda)P$ makes it easier to
restrict $\lambda$ to a computable set. Petkovšek’s algorithm works as
follows:

  1. Step P1

Petkovšek proves that it suffices to consider $A$, $B$ where $A|a_{0}$ and
$B|\tau^{1-n}(a_{n})$. This leaves a finite set of candidates for $A/B$.

  2. Step P2

Compute candidates for $c$.

  3. Step P3

For each candidate $\lambda=cA/B$:
Construct an operator $L_{\lambda}$ whose solutions are the solutions of $L$
divided by $\operatorname{hyp}(\lambda)$. For all polynomial solutions $P$ of
$L_{\lambda}$: join $\operatorname{hyp}(\lambda)P$ to the output.

Consider a $n$-dimensional system over $C(x)$:

$\tau(Y)=MY,\qquad\text{where}\quad M\in\operatorname{GL}_{n}(C(x)).$ (sys)

The goal of this paper is to design efficient algorithms for finding its
hypergeometric solutions. A hypergeometric solution has the form
$\operatorname{hyp}(\lambda)P$ where $P$ is an $n$-dimensional column vector
over $C(x)$. If $\tilde{P}=rP$ then
$\operatorname{hyp}(\lambda)\tilde{P}=\operatorname{hyp}(\lambda\frac{\tau(r)}{r})P$
so we may assume without loss of generality that $P$ is in $C[x]^{n}$ and is
primitive (the content, the gcd of entries, is 1).

Before his sudden passing, Bronstein observed that Petkovšek’s strategy works
for systems as well [10] by replacing the bounds $a_{0}$ and
$\tau^{1-n}(a_{n})$ for $A$ and $B$ with the denominators of $M^{-1}$ and $M$.
Thus, the Bronstein-Petkovšek strategy for hypergeometric solutions is as
follows:

  1. BP1

Construct a set of candidates for $A/B$:

$\mathcal{S}:=\\{\frac{A}{B}:A,B\in C[x]\text{ are monic},\
A\mid\operatorname{denom}(M^{-1}),\ B\mid\operatorname{denom}(M)\\}.$

  2. BP2

For each candidate $A/B$ compute candidates for $c$.

  3. BP3

For each $\lambda=cA/B$: compute all polynomial solutions $P\in C[x]^{n}$ of
$\tau(P)=\lambda^{-1}MP$ and join $\operatorname{hyp}(\lambda)P$ to the
output.

Remarkably, Bronstein’s proof for Step BP1 (Theorem 2) is actually easier than
Petkovšek’s proof for Step P1, despite the fact that it is more general. To
obtain an algorithm, Bronstein still needed a way to find candidates for $c$,
which we will give in subsection 3.2.

The BP-strategy also applies to other cases (such as $q$-difference systems,
or the multi-basic case) provided that one can compute (1) candidates for $c$,
and (2) polynomial solutions.

For operators $L\in C(x)[\tau]$, the paper [16] addressed the following issues
in Petkovšek’s algorithm:

  1. (a)

The number of candidates $A/B$ can be much larger than it needs to be.

  2. (b)

As a side effect, the algorithm can produce duplicate solutions.

For systems, the same issues arise in the BP-strategy. A goal in this paper is
to address these issues. Of course one might discard duplicate solutions, or
take steps to prevent them, but that still leaves issue (a). Reducing the
number of candidates as much as possible leads to a more efficient algorithm
and eliminates issue (b) as a side effect. A key idea is that rather than
bounding $A/B$ using the denominators of $M^{-1}$ and $M$, we bound the type
of $A/B$ by bounding local types.

In [4] another algorithm for computing hypergeometric solutions of difference
systems is presented. It consists in reducing the input system to one or more
scalar equations and then use Petkovšek’s algorithm. This is similar to the
cyclic vector method, which limits the dimension for which the algorithm can
be efficient. Experiments show that our algorithm (implementation at [1]) can
handle systems of high dimension, which will be useful for the application in
factoring operators discussed in subsection 7.2.

##  2 Hypergeometric Solutions

######  Definition 1 (Hypergeometric).

Let $\mathbb{F}$ be a difference field, such as $\mathbb{F}=C(x)$ or
$C((x^{-1}))$. A non-zero element $\gamma$ in a universal extension of
$\mathbb{F}$ is called _hypergeometric_ if $\tau(\gamma)=f\gamma$ for some
$f\in\mathbb{F}^{*}$. In this case let $\operatorname{hyp}(f)$ denote
$\gamma$. The notation is unique up to a $\tau$-constant. When $\gamma$ is
hypergeometric, call the column vector $\gamma R$ (where $R\in\mathbb{F}^{n}$)
a _hypergeometric vector_.

It is not hard to verify, for $f,f_{1},f_{2}\in\mathbb{F}^{*}$:

  * (i)

$\operatorname{hyp}(f_{1})\operatorname{hyp}(f_{2})=\operatorname{hyp}(f_{1}f_{2})$
up to a constant

  * (ii)

$\operatorname{hyp}(\frac{\tau(f)}{f})=f$ up to a constant

The formula

$\operatorname{hyp}(f)=g\operatorname{hyp}(f\frac{g}{\tau(g)}),$

shows that a hypergeometric vector can have different representations. A
hypergeometric vector over $C(x)$ can always be written as $\gamma P$ where
$\gamma$ is hypergeometric and $P\in C[x]^{n}$ is primitive. This is called
the _standard representation_ of the hypergeometric vector.

##  3 Algorithm Version I

In this section a basic version of the algorithm is presented, which follows
the procedure in section 1. We first give the algorithm and then explain why
it works in the follow-up sections.

Algorithm: Version I
Input: $\tau Y=MY$, where $M\in\mathop{\rm GL}_{n}(C(x))$
Output: hypergeometric solutions of the system

  * BP1

Factor $d_{1}:=\operatorname{denom}(M)$ and
$d_{2}:=\operatorname{denom}(M^{-1})$ in $C[x]$ to compute

$\mathcal{S}:=\\{\frac{A}{B}:A,B\in C[x]\text{ are monic},A\mid d_{2},B\mid
d_{1}\\}.$

In subsection 3.1 we show that $\mathcal{S}$ contains all $\frac{A}{B}$ that
are needed.

  * BP2

Use the algorithm introduced in subsection 6.3 to compute

$G:=\\{\text{unramified generalized exponents of the system}\\}.$

Let $\mathcal{H}:=$

$\\{c\frac{A}{B}:cx^{s}(1+dx^{-1})\in
G,\frac{A}{B}\in\mathcal{S},\deg(A)-\deg(B)=s,d-\operatorname{slc}(\frac{A}{B})\in\mathbb{N}\\}.$

Generalized exponents and $\operatorname{slc}$ will be defined in subsection
3.2. An unramified generalized exponent is in the form $cx^{s}(1+dx^{-1})$,
where $(s,c,d)\in\mathbb{Z}\times{C}^{*}\times{C}$.

  * BP3

Let $\mathrm{Sols}=\emptyset$.
For each $c\frac{A}{B}\in\mathcal{H}$, solve the system
$\tau(P)=c^{-1}\frac{B}{A}MP$ for (a basis of) polynomial solutions using the
algorithm introduced in [3] (see also [2]). Add
$\operatorname{hyp}(c\frac{A}{B})P$ to $\mathrm{Sols}$ for any polynomial
solution $P$.
Return $\mathrm{Sols}$ as the output.

###  3.1 Step BP1

######  Theorem 2.

Suppose $Y=\operatorname{hyp}(f)P$ is a hypergeometric solution of Equation
sys that is in the standard form, where $f\in C(x)$. Then

$\operatorname{numer}(f)\mid\operatorname{denom}(M^{-1}),\quad\operatorname{denom}(f)\mid\operatorname{denom}(M).$

The theorem says that any hypergeometric solution $Y$ can be written as
$Y=\operatorname{hyp}(c\frac{A}{B})P$ for some primitive polynomial vector
$P$, some $\tau$-constant $c$, and some $\frac{A}{B}\in\mathcal{S}$, with
$\mathcal{S}$ from Step BP1.

The result holds not only for $C(x)$, but also for the field of fractions of a
general difference ring if it is a UFD.

###### Proof.

Write now $M=d_{1}^{-1}W$ where $d_{1}=\operatorname{denom}(M)$. Then,
substitution of this and $Y=\gamma P$ in standard form into Equation sys
allows us to rewrite the equation as

$\frac{\tau(\gamma)}{\gamma}\tau(P)=\frac{1}{d}_{1}WP,\qquad\text{or,}\qquad
d_{1}A\tau(P)=BWP$

using $\tau(\gamma)/\gamma=A/B$ and clearing denominators. Now $B$ divides the
right-hand side of the equation. On the left-hand side, we may assume $A$ and
$B$ are coprime, and hence $B$ must divide $d_{1}\tau(P)$. But $P$ and hence
$\tau(P)$ is primitive, and so $B\mid d_{1}$.

Write the inverse of $M$ as $M^{-1}={d_{2}}^{-1}W^{\prime}$ where $d_{2}\in
C[x]$ and $W^{\prime}\in C[x]^{n\times n}$. Then, multiplying Equation sys by
$M^{-1}$ and substituting again $Y=\gamma P$ in standard form, we obtain

$\frac{1}{d_{2}}W^{\prime}\tau(P)=\frac{B}{A}P\qquad\text{or}\qquad
AW^{\prime}\tau(P)=d_{2}BP.$

As before, we can derive $A\mid d_{2}$. ∎

######  Example 3.

Let $M=$

$\left(\begin{smallmatrix}0&0&\tfrac{x+1}{x}&0&0&0\\\
0&0&0&0&\tfrac{x+1}{x}&0\\\ 0&0&0&0&0&\tfrac{x+1}{x}\\\
\tfrac{(x+1)^{2}}{x(x+2)(x+3)}&0&-\tfrac{(x+1)(x^{4}+2x^{3}+x^{2}+x+4)}{x(x+2)^{2}(x+3)}&0&-\tfrac{x+1}{x(x+2)^{2}(x+3)}&0\\\
0&\tfrac{(x+1)^{2}}{x(x+2)(x+3)}&\tfrac{(x+1)^{2}}{x(x+2)^{2}(x+3)}&0&0&-\tfrac{x+1}{x(x+2)^{2}(x+3)}\\\
0&0&0&\tfrac{(x+1)^{2}}{x(x+2)(x+3)}&\tfrac{(x+1)^{2}}{x(x+2)^{2}(x+3)}&\tfrac{(x+1)(x^{4}+2x^{3}+x^{2}+x+4)}{x(x+2)^{2}(x+3)}\\\
\end{smallmatrix}\right)$

Then $\operatorname{denom}(M)=x(x+2)^{2}(x+3)$ and
$\operatorname{denom}(M^{-1})=(x+1)^{2}(x+2)$ and so

$\mathcal{S}=\\{x^{i_{0}}(x+1)^{i_{1}}(x+2)^{i_{2}}(x+3)^{i_{3}}:-1\leqslant
i_{0}\leqslant 0,0\leqslant i_{1}\leqslant 2,-2\leqslant i_{2}\leqslant
1,-1\leqslant i_{3}\leqslant 0\\},$

whose cardinality is

$2\cdot 3\cdot 4\cdot 2=48.$

This number is larger than necessary; a provably complete search can be done
with just 8 cases, see 18 which addresses issues (a) and (b) stated in section
1.

######  Remark 4.

Finding hypergeometric solutions of the system in Example 3 is related to the
factorization of the operator

$(x+2)^{2}(x+3)\tau^{4}+\tau^{3}+(x^{4}+2x^{3}+x^{2}+x+4)\tau^{2}+(x+1)\tau+(x+1)(x+2).$

It is explained in subsection 7.2 how to convert a factorization problem into
solving a system for hypergeometric solutions.

###  3.2 Step BP2: generalized exponents

The key of Step BP2 is the notion of _generalized exponents_ (see [12, Section
3.2] for a definition). We can define generalized exponents of a system in a
similar way. We start with some facts about the difference field
$\mathbb{C}((x^{-1}))$.

####  3.2.1 The difference field $K$

Denote $t=1/x$. Let $K=\mathbb{C}((t))$ and
$K_{r}=\mathbb{C}((t^{\frac{1}{r}}))$ for $r=1,2,3,\ldots$. According to [17]
the algebraic closure of $K$ is

$\overline{K}=\bigcup_{r=1}^{\infty}K_{r}.$

We will also denote $\overline{K}$ as $K_{\infty}$. The canonical $t$-adic
valuation on $K$ extends to $K_{r}$ naturally; we denote it by
$v:\overline{K}\mapsto\mathbb{Q}\cup\\{\infty\\}$. For $v$ to be a valuation,
it has to satisfy the following properties:

  * •

$v(a)=\infty$ if and only if $a=0$,

  * •

$v(ab)=v(a)+v(b)$,

  * •

$v(a+b)\geqslant\min\\{v(a),v(b)\\}$, with equality if $v(a)\neq v(b)$.

For a vector with entries in $\overline{K}$, define its valuation to be the
smallest valuation of the entries. The big O and little-o notations are used
for elements in $\overline{K}$ with respect to the valuation $v$. In
particular, when we write $f=g+o(t^{u})$ for $f,g\in\overline{K}$, it means
$v(f-g)>u$.

A general non-zero element in $K_{r}$ factors into
$ct^{s}(1+\sum_{i=1}^{\infty}a_{i}t^{\frac{i}{r}}).$ Call $c$ its _leading
coefficient_ and $ct^{s}$ the _leading term_.

The action of $\tau$ on $K$ is given by

$\tau(t)=\tau(\frac{1}{x})=\frac{1}{x+1}=\frac{t}{1+t}=t-t^{2}+t^{3}-\cdots,$

which extends to $K_{r}$ and $\overline{K}$ following

$\tau(t^{\frac{1}{r}})=t^{\frac{1}{r}}(1+t)^{-\frac{1}{r}}=t^{\frac{1}{r}}(1-\frac{1}{r}t+\frac{(-\frac{1}{r})(-\frac{1}{r}-1)}{2}t^{2}-\cdots).$

We briefly describe the universal extension of $\overline{K}$. For more
details see [15, Chapter 6]. Denote by $\overline{K}\\{\operatorname{hyp}\\}$
the algebra over $\overline{K}$ generated by hypergeometric elements. Let
$\tau$ act naturally on $\overline{K}\\{\operatorname{hyp}\\}$. Then the
polynomial ring $\overline{K}\\{\operatorname{hyp}\\}[l]$ is a universal
extension of $\overline{K}$, equipped with a $\tau$-action following the rule

$\tau(l)=l+t.$

The valuation $v$ extends to $\overline{K}[l]$ by setting
$v(\sum_{i=0}^{n}a_{i}l^{i})=\min\\{v(a_{i}):i=0,1,\ldots,n\\}$.

####  3.2.2 Generalized exponents

This section will show how generalized exponents classify solutions of
operators up to a factor of valuation 0. First consider the group

$\\{\operatorname{hyp}(f):f\in K_{r}^{*}\\}/\\{f\in K_{r}^{*}:v(f)=0\\},$ (1)

which classifies hypergeometric elements over $K_{r}$ up to a factor of
valuation 0. Applying the map $g\mapsto\frac{\tau(g)}{g}$, namely
$\operatorname{hyp}(f)\mapsto f$, the group (1) is isomorphic to

$\mathcal{G}_{r}:=K_{r}^{*}/K_{1,r},$

where $K_{1,r}=\\{\frac{\tau(f)}{f}:f\in K_{r}^{*},v(f)=0\\}$.

######  Lemma 5.

Claim that

$K_{1,r}=\\{g\in K_{r}^{*}:v(g-1)>1\\}.$

###### Proof.

A straight-forward calculation shows that

$\frac{\tau(f)}{f}=1-v(f)t+o(t).$ (2)

Obviously when $v(f)=0$ the right-hand side is $1+o(t)$.

The other direction follows from the proof given in [12, Lemma 3.2.4]

∎

Thus, $f,g\in K_{r}^{*}$ represent the same class in $\mathcal{G}_{r}$ when
$v(\frac{f}{g}-1)>1,$ which happens if and only if $v(f-g)>v(f)+1.$ For

$f=ct^{s}(1+\sum_{i=1}^{\infty}a_{i}t^{\frac{i}{r}})\in
K_{r}^{*},\quad\text{where }r<\infty,$

denote

$\operatorname{Trunc}(f)=ct^{s}(1+\sum_{i=1}^{r}a_{i}t^{\frac{i}{r}})\in
K_{r}^{*}.$

Clearly $\operatorname{Trunc}:K_{\infty}^{*}\rightarrow K_{\infty}^{*}$ is
well-defined. Then $v(f-\operatorname{Trunc}(f))>v(f)+1$, which means the
image of $f$ in $\mathcal{G}_{r}$ is represented by $\operatorname{Trunc}(f)$.
Hence $\mathcal{G}_{r}$ can be identified with the following set of
representatives

$E_{r}:=\operatorname{Trunc}(K_{r}^{*})=\\{ct^{s}(1+\sum_{i=1}^{r}a_{i}t^{\frac{i}{r}}):c\in\mathbb{C},s\in\frac{1}{r}\mathbb{Z},a_{i}\in\mathbb{C}\\}.$

The following is a short exact sequence of abelian groups

$1\rightarrow K_{1,r}\rightarrow
K_{r}^{*}\xrightarrow{\operatorname{Trunc}}E_{r}\rightarrow 1,$

where the set of representatives $E_{r}$ is turned into a group by

$g_{1}\circ g_{2}:=\operatorname{Trunc}(g_{1}g_{2}),\quad
g_{1},g_{2}\in{E}_{r}.$

Denote

$\overline{K}[l]_{\rm
hyp}=\\{\operatorname{hyp}(f)p:f\in\overline{K},p\in\overline{K}[l]\setminus\\{0\\}\\}\subseteq\overline{K}\\{\operatorname{hyp}\\}[l].$
(3)

It is a multiplicative monoid. An element $h\in\overline{K}[l]_{\rm hyp}$ can
be written uniquely in the form $h=\operatorname{hyp}(e)p,$ where $e\in
E_{\infty}$ and $p\in\overline{K}[l]$ has valuation 0. Define
$\mathrm{gen}(h)=e$, the _generalized exponent_ of $h$. Note that
$\mathrm{gen}$ commutes with multiplication.

Generalized exponents of operators are defined in [12]. An order $n$ operator
has exactly $n$ generalized exponents in $E$, counting with multiplicity. The
relevant property for us is:

######  Proposition 6.

For a non-zero operator $L\in\overline{K}[\tau]$, $e\in E_{\infty}$ is a
generalized exponent of $L$ if and only if there exists a solution
$h\in\overline{K}[l]_{\rm hyp}$ with $e=\mathrm{gen}(h)$.

######  Definition 7.

Based on 6, we call $e\in E_{\infty}$ a _generalized exponent_ of the system
$\tau Y=MY$ if this system has a solution of the form $\operatorname{hyp}(e)S$
where $S\in\overline{K}[l]^{n}$ and $v(S)=0$.

######  Remark 8 (Algorithms).

For operators over $\mathbb{C}(x)$ we can quickly compute all generalized
exponents with the program GeneralizedExponents in the LREtools package in
Maple 2021. For systems over $\mathbb{C}(x)\subseteq K$, we will give in
section 6 an algorithm for computing the unramified generalized exponents
(i.e. those in $\mathcal{G}_{1}$). This algorithm has also been implemented in
Maple.

######  Lemma 9.

Denote $e=\operatorname{Trunc}(f)$ for $f\in\overline{K}^{*}$. Due to 5,
$\operatorname{hyp}(\frac{f}{e})$ lies in $\overline{K}$ and has valuation
$0$. Hence if $\operatorname{hyp}(f)S$ is a hypergeometric solution of a
system where $v(S)=0$, then $e$ is a generalized exponent.

Now we discuss the relation between hypergeometric solutions of systems over
$\mathbb{C}(x)$ and their generalized exponents. Since $\mathbb{C}(x)\subseteq
K$, $\operatorname{Trunc}(f)\in\mathcal{G}_{1}$ for $f\in\mathbb{C}(x)^{*}$.
For a monic Laurent series

$f=t^{n}+dt^{n+1}+\cdots\in K^{*},$

denote $\operatorname{slc}(f)=d$, where $\operatorname{slc}$ stands for
_second leading coefficient_.

For a monic polynomial

$A=x^{n}+dx^{n-1}+\cdots,$

we have

$\operatorname{Trunc}(A)=x^{n}+dx^{n-1}=t^{-\deg(A)}(1+\operatorname{slc}(A)t),$

and for a non-zero rational function $c\frac{A}{B}$ where $A,B$ are monic
polynomials, straight-forward computations show that

$\operatorname{Trunc}(c\frac{A}{B})=ct^{\deg(B)-\deg(A)}(1+(\operatorname{slc}(A)-\operatorname{slc}(B))t).$

The following lemma justifies the definition of $\mathcal{H}$ in Step BP2
(section 3).

######  Lemma 10.

Consider Equation sys. Suppose $A,B\in C[x]$ are monic. If there is a solution
in the form $\operatorname{hyp}(cA/B)P$ where $P\in C(x)^{n}$ then the system
has a generalized exponent $ct^{s}(1+dt)\in E_{1}$ that satisfies

$\deg(B)-\deg(A)=s,\quad\operatorname{slc}(A)-\operatorname{slc}(B)-d=v(P)\in\mathbb{Z}.$
(4)

If we further require $P\in C[x]^{n}$ then the relations have a stronger form

$\deg(B)-\deg(A)=s,\quad-\operatorname{slc}(A)+\operatorname{slc}(B)+d=-v(P)=\deg(P)\in\mathbb{N}.$
(5)

We will refer to (4) as the _weak compatibility relations_ and (5) the _strong
compatibility relations_.

###### Proof.

Note that

$\operatorname{hyp}(c\frac{A}{B})P=\operatorname{hyp}(c\frac{A}{B}(1+t)^{-v(P)})\,t^{-v(P)}P,$

where $v(t^{-v(P)}P)=0$. By 9 the existence of such a solution implies

$\operatorname{Trunc}(c\frac{A}{B})\circ\operatorname{Trunc}((1+t)^{-v(P)})=ct^{\deg(B)-\deg(A)}(1+(\operatorname{slc}(A)-\operatorname{slc}(B)-v(P))t)$

is a generalized exponent in $E_{1}$. When $P\in C[x]^{n}$,
$-v(P)=\deg(P)\in\mathbb{N}$. ∎

######  Example 11.

Consider the same system as in 3. There are two unramified generalized
exponents:

$t^{2}(1+2t),\quad t^{-1}(1-4t).$

Therefore $c=1$. Eight $A/B$’s in $\mathcal{S}(M)$ from Step 1 (3) match
$t^{2}(1+2t)$; none matches the other generalized exponent.

###  3.3 Step BP3

In this step, for each potential $\lambda=cA/B$, compute polynomial solutions
of $\lambda\tau(P)=MP$ using the algorithm given in [3].

######  Example 12 (Continued from 11).

Our next step is to find all polynomial solutions of the system

$\tau(P)=(c\frac{A}{B})^{-1}MP$ (6)

for each $c\frac{A}{B}$. The set of hypergeometric solutions is one
dimensional, given by

$\operatorname{hyp}(\frac{x+1}{x(x+2)(x+3)})\begin{pmatrix}(x+1)^{3}(x+2)x^{2}\\\
(x+2)(x+1)(-x-1)\\\ (x+2)^{2}(x+1)^{2}\\\ -x^{4}-4x^{3}-3x^{2}+1\\\ -x-2\\\
(x+2)(x+3)\end{pmatrix}.$

Algorithm Version I computes this solution four times; different
$c\frac{A}{B}$’s can lead to the same solution. If
$\operatorname{hyp}(\lambda)P$ is a solution where $P\in C[x]^{n}$ and
$\operatorname{hyp}(\lambda^{\prime})=f\operatorname{hyp}(\lambda)$ for some
$f\in C[x]$, then Version I may rediscover $\operatorname{hyp}(\lambda)P$ by
computing $\operatorname{hyp}(\lambda^{\prime})fP$.

Computing too many candidates and duplicate solutions – issues (a),(b) in
section 1 – will be addressed in Version II.

##  4 Algorithm Version II

###  4.1 Type and local types

######  Definition 13 (Type).

For $f_{1},f_{2}\in C(x)^{*}$, say $f_{1}$ and $f_{2}$ have the same _type_ if
$\operatorname{hyp}(f_{1}){C}(x)=\operatorname{hyp}(f_{2}){C}(x)$.

Duplicates can occur in Algorithm Version I because the set $\mathcal{S}$ can
contain different $\frac{A}{B}$’s with the same type. If we select one single
$\frac{A}{B}$ for each type, then repeated solutions are avoided. This can be
done by bounding _local types_.

Notation: For a prime polynomial $p\in C[x]$, denote
$[p]:=\\{\tau^{i}(p):i\in\mathbb{Z}\\}$, the set of all polynomials that are
shift equivalent to $p$.

It is easy to detect if two monic irreducible polynomials are shift
equivalent, because the $i$ in $f=\tau^{i}(g)$ can be computed by dividing
$\mathrm{slc}(f)-\mathrm{slc}(g)$ by the degree, and if $i$ is an integer,
then one checks if $f$ is actually equal to $\tau^{i}(g)$.

For a prime polynomial $p\in C[x]$, let
$v_{p}:C(x)\rightarrow\mathbb{Z}\cup\\{\infty\\}$ be the $p$-valuation.

######  Definition 14 (Local Type).

For a non-zero element $a\in C(x)\setminus\\{0\\}$ and a prime polynomial
$p\in C[x]$, let

$g_{p}(a)=\sum_{k\in\mathbb{Z}}v_{\tau^{k}(p)}(a).$

In other words, $g_{p}$ is the sum of valuations with respect to all prime
polynomials that are shift equivalent to $p$. We call $g_{p}(a)$ the _local
type_ of $a$ at $[p]$.

The following theorem shows the relation between types and local types.

######  Theorem 15 ([16, Theorem 1]).

Suppose $\lambda_{1}=c_{1}\frac{A_{1}}{B_{1}},\
\lambda_{2}=c_{2}\frac{A_{2}}{B_{2}}\in C(x)^{*}$ where
$A_{1},B_{1},A_{2},B_{2}\in C[x]$ are monic and $c_{1},c_{2}\in C$. Then
$\lambda_{1}$ and $\lambda_{2}$ have the same type if and only if

  * •

$c_{1}=c_{2}$, and

  * •

$g_{p}(\lambda_{1})=g_{p}(\lambda_{2})$ for any prime $p\in C[x]$.

Let now $Y=\operatorname{hyp}(cA/B)P$ be once more a hypergeometric solution
to Equation sys where $c\in C^{*}$ and $A,B\in C[x]$ are monic. Then Theorem 2
yields $A\mid\operatorname{denom}(M^{-1})$ and $B\mid\operatorname{denom}(M)$.
The first statement implies

$0\leqslant g_{p}(A)\leqslant g_{p}(\operatorname{denom}(M^{-1}))$

for every prime polynomial $p\in C[x]$ while the second statement yields

$0\leqslant g_{p}(B)\leqslant g_{p}(\operatorname{denom}(M)).$

Now we restate Theorem 2 in terms of local types.

######  Lemma 16.

For a hypergeometric solution $Y=\operatorname{hyp}(\lambda)P$ of Equation
sys, we have

$-g_{p}(\operatorname{denom}(M))\leqslant g_{p}(\lambda)\leqslant
g_{p}(\operatorname{denom}(M^{-1}))$

for every prime $p\in C[x]$.

If the number of combinations in the algorithm below is high, then it can be
worthwhile to spend more CPU time to obtain a sharper bound:

$-g_{p}(\operatorname{denom}(\tau(M)M))\leqslant 2g_{p}(\lambda)\leqslant
g_{p}(\operatorname{denom}((\tau(M)M)^{-1})).$ (7)

To see why this holds, note that $Y$ is also a solution of
$\tau^{2}(Y)=\tau(MY)=\tau(M)\tau(Y)=\tau(M)MY$, and $\tau^{2}({\rm
hyp}(\lambda))=\tau(\lambda)\lambda{\rm hyp}(\lambda)$, and
$g_{p}(\tau(\lambda)\lambda)=2g_{p}(\lambda)$. One can derive even sharper
bounds from denominators of even larger products ($\tau^{2}(M)\tau(M)M$, etc.)
for diminishing returns.

###  4.2 Algorithm Version II

Theorem 15 and 16 lead to a more efficient algorithm for computing
hypergeometric solutions.

Algorithm: Version II
Input: $\tau Y=MY$, where $M\in\mathop{\rm GL}_{n}(C(x))$
Output: hypergeometric solutions.

  * BP1

Factor $\operatorname{denom}(M)$ and $\operatorname{denom}(M^{-1})$. Say

$\operatorname{denom}(M)=a_{1}\prod_{i=1}^{m}p_{i}^{e_{i}},\quad\operatorname{denom}(M^{-1})=a_{2}\prod_{i=1}^{m}p_{i}^{e_{i}},$

where $a_{1},a_{2}\in C-\\{0\\}$ and $p_{i}$ are monic irreducible
polynomials.
Sort $p_{1},\ldots,p_{m}$ according to their shift equivalence classes.
Suppose $[p_{1}],[p_{2}],\ldots,[p_{l}]$ are all the shift equivalence
classes.
For $i=1,2,\ldots,l$, calculate $g_{p_{i}}(\operatorname{denom}(M))$ and
$g_{p_{i}}(\operatorname{denom}(M^{-1}))$.
Return
$\mathcal{S}_{2}:=\\{\prod_{i=1}^{l}p_{i}^{f_{i}}:-g_{p_{i}}(\operatorname{denom}(M))\leqslant
f_{i}\leqslant g_{p_{i}}(\operatorname{denom}(M^{-1}))\\}$.

  * BP2

Use the algorithm introduced in subsection 6.3 to compute

$G:=\\{\text{unramified generalized exponents of the system}\\}.$

Let $\mathcal{H}_{2}:=$

$\\{c\frac{A}{B}:cx^{s}(1+dx^{-1})\in
G,\frac{A}{B}\in\mathcal{S}_{2},\deg(A)-\deg(B)=s,d-\operatorname{slc}(\frac{A}{B})\in\mathbb{Z}\\}.$

  * BP3

Let $\mathrm{Sols}=\emptyset$.
For each $c\frac{A}{B}\in\mathcal{H}_{2}$, solve the system
$\tau(P)=c^{-1}\frac{B}{A}MP$ for (a basis of) _rational_ solutions (e.g. [3,
2]). Add $\operatorname{hyp}(c\frac{A}{B})P$ to $\mathrm{Sols}$ for any
rational solution $P$.
Return $\mathrm{Sols}$ as the output.

######  Remark 17.

Steps BP2 and Step BP3 are slightly different from the same steps in Version
I. The reason is, if $\operatorname{hyp}(cA/B)P$ is a hypergeometric solution
in standard form, then Version II only computes the type of $cA/B$, which
determines $\operatorname{hyp}(cA/B)$ up to a factor in ${C}(x)$. As a result,
$P$ is no longer guaranteed to have polynomial entries. Hence in Step BP2 the
weak compatibility relations (4) are applied instead of the strong and in Step
BP3 we compute rational solutions.

######  Example 18.

Let $M$ be the same matrix as in 3. Every irreducible factor of
$\operatorname{denom}(M)=x(x+2)^{2}(x+3)$ and
$\operatorname{denom}(M^{-1})=(x+1)^{2}(x+2)$ is shift equivalent to $x$, and

$g_{x}(\operatorname{denom}(M))=4,\quad
g_{x}(\operatorname{denom}(M^{-1}))=3.$

Therefore,

$\mathcal{S}_{2}(M)=\\{x^{i}:-4\leqslant i\leqslant 3\\}.$

The cardinality of $\mathcal{S}_{2}(M)$ is 8, a significant improvement over
Version I (3). In 11 we calculated the unramified generalized exponents:

$t^{2}(1+2t),\quad t^{-1}(1-4t).$

The matching elements in $\mathcal{S}_{2}(M)$ are $x^{-2},x$ and so
$\mathcal{H}_{2}(M)=\\{x^{-2},x\\}$.

Version II has two advantages over Version I:

  * •

More efficient: there are fewer candidates for $\frac{A}{B}$.

  * •

No duplicate solutions.

as well as one drawback: we now have to compute rational solutions $P$ instead
of polynomial solutions (17). But this drawback is easily overcome, the next
section shows how one can ensure that $P$ will be in $C[x]^{n}$.

##  5 Algorithm Version III

Algorithm Version II greatly reduced the number of $c\frac{A}{B}$’s by
avoiding repeated types. This section is devoted to the technical result that
we can further improve the algorithm by choosing $A/B$ in a way such that in
Step BP3 we only need to search for polynomial solutions instead of rational
ones. The idea is to make sure $A/B$ satisfies the conditions in Theorem 2 and
at the same time minimize $\operatorname{slc}(A)-\operatorname{slc}(B)$. A few
notations are needed before stating the main result of this section.

Suppose $r_{1},r_{2}\in C(x)$ are of the same type. Say $r_{1}\preccurlyeq
r_{2}$ if $\frac{\operatorname{hyp}(r_{2})}{\operatorname{hyp}(r_{1})}\in
C[x]$. It is easy to see that $\preccurlyeq$ is a partial order (on the set of
all non-zero rational functions or those that are of the same type). Theorem 2
states that every hypergeometric solution can be written as
$\operatorname{hyp}(c\frac{A}{B})P$ where $P\in C[x]^{n}$. If
$\operatorname{hyp}(c^{\prime}\frac{A^{\prime}}{B^{\prime}})P^{\prime}=\operatorname{hyp}(c\frac{A}{B})P$
and $c^{\prime}\frac{A^{\prime}}{B^{\prime}}\preccurlyeq c\frac{A}{B}$, then a
consequence is $P^{\prime}\in C[x]^{n}$. Thus, to achieve the goal, we want
$\frac{A}{B}$ to be as small as possible.

Recall that

$\mathcal{S}(M)=\\{\frac{A}{B}:A,B\in C[x]\text{ monic,
}A\mid\operatorname{denom}(M^{-1}),B\mid\operatorname{denom}(M)\\}$

contains all potential $\frac{A}{B}$ in Algorithm Version I (section 3).
Consider the partition $\mathcal{S}(M)=\bigcup_{i}\mathcal{T}_{i}(M)$ of
$\mathcal{S}(M)$, where each $\mathcal{T}_{i}(M)$ consists of all elements of
a particular type.

######  Theorem 19.

Suppose $s_{i}\in\mathcal{T}_{i}(M)$ has the smallest slc in
$\mathcal{T}_{i}(M)$. Then $s_{i}$ is the $\preccurlyeq$-smallest element in
$\mathcal{T}_{i}(M)$.

###### Proof.

All elements in $\mathcal{T}_{i}(M)$ have the same local type, so the local
types of $\mathcal{T}_{i}(M)$ is well-defined. Suppose $\mathcal{T}_{i}(M)$
has a non-zero local type $e_{1},e_{2},\ldots,e_{N}$ at
$[p_{1}],[p_{2}],\ldots,[p_{N}]$, respectively, where
$p_{1},p_{2},\ldots,p_{N}$ are mutually shift non-equivalent. We first
consider the case $N=1$.

Assume $s_{i}$ is not the $\preccurlyeq$-smallest element in
$\mathcal{T}_{i}(M)$. Then there exists $r\in\mathcal{T}_{i}(M)$ such that
$\operatorname{hyp}(\frac{r}{s_{i}})\notin C[x]$. Therefore there exists an
irreducible polynomial $p\in C[x]$ such that

$\operatorname{hyp}(\frac{r}{s_{i}})=\frac{f}{pg},$

where $f,g\in C[x]$ and $\gcd(f,pg)=1$. Then
$s_{i}\frac{p}{\tau(p)}\in\mathcal{T}_{i}(M)$ and

$\mathrm{slc}(s_{i}\frac{p}{\tau(p)})=\mathrm{slc}(s_{i})+\mathrm{slc}(p)-\mathrm{slc}(\tau(p))=\mathrm{slc}(s_{i})-\deg(p),$

which contradicts to the assumption that $s_{i}$ has the smallest
$\mathrm{slc}$ among elements in $\mathcal{T}_{i}(M)$.

The $\preccurlyeq$-smallest element is unique because the ordering is anti-
symmetric. The proof also works for $N>1$, treating each shift-equivalent
class separately. ∎

######  Remark 20.

If $ct^{s}(1+dt)$ is the unramified generalized exponent with the largest $d$
that is compatible with $\operatorname{hyp}(c\frac{A}{B})P$ where $P\in
C[x]^{n}$, then $d-\operatorname{slc}(A/B)$ is a degree bound for $P$ due to
the relation $d=\operatorname{slc}(A/B)+\deg(P)$. With a degree bound the
problem of finding polynomial solutions reduces to solving a system of linear
equations.

######  Example 21.

Let $M$ be the same as in 3. The unramified generalized exponents are:

$t^{2}(1+2t),\quad t^{-1}(1-4t).$

There are eight types of elements in $\mathcal{S}(M)$, two of them matching
the generalized exponents, as was stated in 18. The $\preccurlyeq$-smallest
elements in $\mathcal{S}(M)$ of these two types are

$s_{1}=\frac{x+1}{(x+2)^{2}(x+3)}\quad s_{2}=\frac{(x+1)^{2}}{x+3}.$

The set of candidates for $c\frac{A}{B}$ is $\\{s_{1},s_{2}\\}$. The degree
bounds from Remark 20 are $8$ for $s_{1}$ and $-3$ for $s_{2}$. So we discard
$s_{2}$ and compute polynomial solutions for $\tau(P)=\frac{1}{s_{1}}MP$,
finding

$\begin{pmatrix}x^{3}(x+2)(x+1)^{4}\\\ -x(x+2)(x+1)^{3}\\\
x(x+2)^{2}(x+1)^{3}\\\ -x(x+1)(x^{4}+4x^{3}+3x^{2}-1)\\\ -x(x+2)(x+1)\\\
x(x+3)(x+2)(x+1))\end{pmatrix}.$

##  6 Computing the unramified generalized exponents of a system

The difference operator $\delta=\tau-1$ is given by
$\delta(f)=\tau(f)-f=f(x+1)-f(x).$

###  6.1 Local Hypergeometric Solutions

The difference field $({C}{(\\!(}x^{-1}{)\\!)},\tau)$ is an extension of
$({C}(x),\tau)$, so any system of first order linear difference equations,
$\tau(Y)=MY$ with $M\in\mathop{\rm GL}_{n}({C}(x))$, can be viewed as a
difference system over ${C}{(\\!(}x^{-1}{)\\!)}$. Moreover, if it has a
hypergeometric solution over ${C}(x)$

$\operatorname{hyp}(c\frac{A}{B})\,P$

where $c\in{C}^{*}$, $A,B\in{C}[x]$ and $P\in{C}[x]^{n}$ then it has a formal
solution of the form

$Y=\Gamma(x)^{s}c^{x}x^{d}F$ (8)

where $s\in\mathbb{Z}$, $c\in{C}^{*}$, $d\in{{C}}$ and
$F\in{{C}}[[x^{-1}]]^{n}$ with $\mathop{v}(F)=0$:

$F=\left(F_{0}+x^{-1}F_{1}+\dots\right)\,\quad F_{0}\neq 0.$

More precisely, one has the relations

$s=\deg A-\deg B\in\mathbb{Z},$
$d-\operatorname{slc}(A)+\operatorname{slc}(B)=\deg P\in\mathbb{N}.$

Note that such a formal solution is in fact a hypergeometric solution over
${C}{(\\!(}x^{-1}{)\\!)}$. Indeed, it is a solution of the form $\gamma F$
where $\gamma=\Gamma(x)^{s}c^{x}x^{d}$ is hypergeometric over
${C}{(\\!(}x^{-1}{)\\!)}$, since
$\frac{\tau(\gamma)}{\gamma}=cx^{s}(1+\frac{1}{x})^{d}=cx^{s}(1+\frac{d}{x}+\cdots)\in{{C}}{(\\!(}x^{-1}{)\\!)}$.
For this reason a formal solution of the form (8) will be called a local
hypergeometric solution (lhs, in short) of our system. Thus, a hypergeometric
solution over $C(x)$ is also a local hypergeometric solution. The set $G$ of
unramified generalized exponents of the system can be identified with the set
of all triplets $(s,c,d)\in\mathbb{Z}\times{C}^{*}\times{C}$ involved in
possible lhs (which we still denote by $G$).

We will explain in this section how to determine this set $G$, directly from
the system (i.e. without using cyclic vectors).

Note that as in the scalar case ($L(y)=0$ with $L\in{C}[x][\tau]$), the
possible values of $s$ are among the integer slopes (if any) of the
$\tau-$Newton polygon of $L$ [7]. To each slope $s$ corresponds a Newton
polynomial $P_{s}(\lambda)\in{C}[\lambda]$ whose nonzero roots (if any) in
${C}$ are candidates for the $c$’s. For each pair $(s,c)$ the corresponding
$d$’s are obtained as roots in ${C}$ of the so-called indicial equation (or
Newton polynomial) associated to the slope $0$ (if any) of the $\delta-$Newton
polygon of the operator $L_{s,c}$ obtained by performing in the equation
$L(y)=0$ the change of variable $y\mapsto(\Gamma(x)^{s}c^{x})\cdot y$.

Let us note that in the operator case the local data (Newton polygon and
polynomials) can be easily obtained from the coefficients of the operator (see
[7, 8], [13]). However, in the system case the local data cannot be obtained
directly from $M$ nor from its inverse. Our strategy consists in reducing the
given system to a suitable gauge-equivalent one from which we can easily
obtain the integer slopes of the Newton polygon and their corresponding Newton
polynomials. Note that in the systems case, we do not compute the Newton
polygon, we only need its integer slopes and those are determined in another
way, using the properties of a super-reduced form of the system.

Before explaining our method we begin by setting up some notation. Any
difference system

$\tau(Y)=MY,\quad\textit{where }M\in\mathrm{GL}_{n}({C}(x))$ (9)

can be rewritten as

$\displaystyle\delta(Y)=\widehat{M}Y,\quad\textit{where
}\widehat{M}:=M-I_{n}\in{C}(x)^{n\times n}.$ (10)

Performing the gauge transformation $Y=TZ$ where $T\in\mathop{\rm
GL}_{n}({C}(x))$ yields the following _gauge-equivalent_ system

$\delta(Z)=\widehat{N}Z\quad\;\;\widehat{N}=T_{\delta}[\widehat{M}]:=\tau(T^{-1})(\widehat{M}T-\delta(T)).$

If $Y=\Gamma(x)^{s}c^{x}x^{d}F$ with $F_{0}\neq 0$ is a lhs of the system
$\tau(Y)=MY$ then for all $T\in\mathop{\rm GL}_{n}({C}(x))$, $Z=TY$ is lhs of
the equivalent system $\tau(Z)=T_{\delta}[M]Z$. Indeed,

$Z=TY=\Gamma(x)^{s}c^{x}x^{d}TF=\Gamma(x)^{s}c^{x}x^{d+\nu}\widetilde{F}$

where $\widetilde{F}\in{{C}}[[x^{-1}]]^{n}$ with $\widetilde{F}_{0}\neq 0$ and
where $\nu\in\mathbb{Z}$ with $\nu\geqslant\mathop{v}(T)$. Thus, the pairs
$(s,c)$ are invariants under gauge transformations, while the $d$’s may change
into $d+\nu$, for some $\nu\in\mathbb{Z}$.

We are interested in the computation of all the triplets
$(s,c,d)\in\mathbb{Z}\times{{C}}^{*}\times{{C}}$ that appear in a lhs of
Equation sys. In general, those triplets cannot be obtained directly from the
matrix system $M$. However, we will show that

  * -

first, all the possible pairs $(s,c)$ can be obtained from a super-irreducible
form, in the sense of [5, 6], of the matrix $\widehat{M}$;

  * -

then, once a pair $(s,c)$ is obtained, the corresponding indices (the $d$’s)
can be computed from a $0$-simple form ([9]) of the system $\delta
Y=(x^{-s}c^{-1}M-I_{n})Y$.

###  6.2 Computing the pairs $(s,c)$

For a lhs $Y=\Gamma(x)^{s}c^{x}x^{d}F$ we will refer to $s$ as the slope, $c$
as the characteristic root and $d$ as the index of $Y$. The following facts
are easy to prove

  1. (a)

The slope $s$ of any lhs of the system $\tau(Y)=MY$ satisfies the following
inequalities:

$\mathop{v}(M^{-1})\leqslant s\leqslant-\mathop{v}(M).$

  2. (b)

If $Y$ is a lhs of slope $s$ then for any integer $p$,
$\widetilde{Y}=\Gamma(x)^{p}Y$ is a lhs, with slope $\tilde{s}:=s+p$, of the
difference system $\tau(\widetilde{Y})=x^{p}M\widetilde{Y}$. Hence, if we
choose $p>-\mathop{v}(M^{-1})$ then we can make
$\tilde{s}\geqslant\mathop{v}(x^{-p}M^{-1})>0$. Thus, even if it means to
replace $M$ by $x^{p}M$ with $p=-\mathop{v}(M^{-1})+1$, one can suppose that
$\mathop{v}(M^{-1})\geqslant 1$ and in this case $\mathop{v}(M)\leqslant-1$
and the slopes of the local hypergeometric solutions of (9) are positive.

In the sequel, we assume that $\mathop{v}(M^{-1})\geqslant 1$ so that
$\mathop{v}(M)\leqslant-1$ and the slopes of the local hypergeometric
solutions of (9) are positive. It follows that the matrix
$\widehat{M}=M-I_{n}$ has the same valuation as M.

Put $q:=-\mathop{v}(\widehat{M})\geqslant 1$ and let $1\leqslant s\leqslant q$
be an integer. We will look for a necessary and sufficient condition for $s$
to be a slope of some lhs of (9) that can be read easily from the first few
coefficients in the expansion of $\widehat{M}$.

For $1\leqslant i\leq n$, let
$\nu^{(s)}_{i}=\min{(\mathop{v}(\widehat{M}_{i,.}),-s)}$ where
$\widehat{M}_{i,.}$ denotes the $i$th row of the matrix $\widehat{M}$. Put

$D^{(s)}:=\operatorname{diag}(x^{\nu^{(s)}_{1}},\ldots,x^{\nu^{(s)}_{n}})\;\;\hbox{
and }\;\;N^{(s)}:=D^{(s)}\widehat{M}.$

It is clear that $\mathop{v}{(N^{(s)})}\geqslant 0$ and that
$\mathop{v}{(D^{(s)})}\geqslant-s$ (in fact, the entries of the matrix
$x^{s}D$ are polynomials in $x^{-1}$). It then makes sense to set

$D^{(s)}=x^{-s}\left(D^{(s)}_{0}+\mathop{\rm O}(x^{-1})\right)\;\;\hbox{ and
}\;\;N^{(s)}=N^{(s)}_{0}+\mathop{\rm O}(x^{-1}).$

Our system (9) can be rewritten as

$\displaystyle D^{(s)}\delta(Y)-N^{(s)}Y=0.$

Let us now replace $Y$ by the expression in (8). One has

$\delta(Y)=\tau\left(\Gamma(x)^{s}c^{x}x^{d}\right)\delta(F)+\delta\left(\Gamma(x)^{s}c^{x}x^{d}\right)\,F.$

Now, we have $\delta(F)=F(x+1)-F(x)=\mathop{\rm O}(x^{-2})$ and

$\tau\left(\Gamma(x)^{s}c^{x}x^{d}\right)=cx^{s}\left(1+x^{-1}\right)^{d}\Gamma(x)^{s}c^{x}x^{d}\\\
=cx^{s}\left(1+dx^{-1}+\mathop{\rm O}(x^{-2})\right)\Gamma(x)^{s}c^{x}x^{d}.$

Since $s>0$, one has

$\delta\left(\Gamma(x)^{s}c^{x}x^{d}\right)=x^{s}\left(c+\mathop{\rm
O}(x^{-1})\right)\Gamma(x)^{s}c^{x}x^{d},$

hence

$\delta\left(\Gamma(x)^{s}c^{x}x^{d}F\right)=x^{s}\Gamma(x)^{s}c^{x}x^{d}\left(cF_{0}+\mathop{\rm
O}(x^{-1})\right).$

Now

$D^{(s)}\delta(Y)-N^{(s)}Y\\\
=\left(\Gamma(x)^{s}c^{x}x^{d}\right)\left(x^{s}D^{(s)}\left(cF_{0}+\mathop{\rm
O}(x^{-1})\right)-N^{(s)}\left(F_{0}+\mathop{\rm O}(x^{-1})\right)\right).$

Thus

$\displaystyle
D^{(s)}\delta(Y)-N^{(s)}Y=\Gamma(x)^{s}c^{x}x^{d}\left((cD^{(s)}_{0}-N^{(s)}_{0})F_{0}+\mathop{\rm
O}(x^{-1})\right).$

If a nonzero formal solution of the form (8) exists then one must have

$(cD^{(s)}_{0}-N^{(s)}_{0})F_{0}=0,$

and therefore

$\det{(cD^{(s)}_{0}-N^{(s)}_{0})}=0,$ (11)

since $F_{0}\neq 0$.

Thus, if the system (10) has a solution of the form (8) then $c$ must
annihilate the polynomial $E_{s}(\lambda):=\det{(N^{(s)}_{0}-\lambda
D^{(s)}_{0})}$. However, it may happen that this polynomial vanishes
identically in $\lambda$, in which case it gives no useful information. This
motivates the following definition

######  Definition 1.

System (10) is said to be $s$-simple if the polynomial
$E_{s}(\lambda)=\det{(N^{(s)}_{0}-\lambda D^{(s)}_{0})}\not\equiv 0$. In this
case the polynomial $E_{s}(\lambda)$ will be called the $s$-characteristic
polynomial of (10).

It is important to notice that any system can be reduced to an $s$-simple
gauge-equivalent one. An algorithm to perform such a reduction is presented in
[9] where it is shown that a super-irreducible system is necessarily
$s$-simple for all $s\in\\{1,\dots,q\\}$. Hence, it is sufficient to compute a
super-irreducible form of (10) to get all the $E_{s}(\lambda)$ for
$s=1,\dots,q$. This can be done using one of the super-reduction algorithms
presented in [6, 5].

Now, for $1\leqslant s\leqslant q$, let $C_{s}$ be the set of nonzero roots of
$E_{s}(\lambda)$ in ${C}$. If $C_{s}=\emptyset$ then there is no lhs of slope
$s$ and this value of $s$ can be discarded; otherwise the set $\\{s\\}\times
C_{s}$ gives candidates for the pairs $(s,c)$ which occur in the possible
hypergeometric solutions with slope $s$. In summary, the following algorithm
produces a set of candidates for the pairs $(s,c)$ occurring in the possible
hypergeometric solutions of the input system.

Algorithm: sc-Pairs
Input: $M\in\mathop{\rm GL}_{n}({C}(x))$
Output: $SC$ the set of pairs $(s,c)$

  * 1.

If $\mathop{v}{(M^{-1})}\leqslant 0$ then set $p:=-\mathop{v}{(M^{-1})}+1$ and
replace the matrix $M$ by $x^{p}M$.

  * 2.

Let $\widehat{M}:=M-I_{n}$ and $q:=-\mathop{v}{(M)}$.

  * 3.

Apply the super-reduction algorithm to the matrix $\widehat{M}$. This produces
the list of the polynomials $E_{s}(\lambda)$ for $s=1,\dots,q$.

  * 4.

$SC:=\emptyset$.

  * 5.

For $s=1,\dots,q$ do

    * –

Compute the set $C_{s}:=\\{c\in{C}\setminus\\{0\\}\;|\;E_{s}(c)=0\\}$.

    * –

If $C_{s}\neq\emptyset$ then $SC:=SC\cup\\{s-p\\}\times C_{s}$.

###  6.3 Computing the exponents $d$’s associated to a pair $(s,c)$

Let a system (9) be given. We assume that we have computed the set $SC$ of the
$(s,c)$-pairs belonging to lhs of (9) and that $SC$ is not empty. Let
$(s,c)\in SC$. We will now explain how to determine the set of the indices
$d\in{C}$ such that the system possesses a lhs of the form (8).
If $Y$ is such a solution then $\Gamma(x)^{-s}c^{-x}Y=x^{d}F$ is a formal
solution of the system $\tau Z=M_{s,c}Z$ where $M_{s,c}:=c^{-1}x^{-s}M$ and
hence $d$ (modulo $\mathbb{Z}$) must satisfy an algebraic equation which can
be computed from a so called simple form of the matrix $\widehat{M}_{s,c}$.

Here again the system $\delta(Z)=\widehat{M}_{s,c}Z$ is rewritten as

$\displaystyle xD\delta(Z)-NZ=0,$ (12)

with

$D:=\operatorname{diag}(x^{\nu_{1}},\ldots,x^{\nu_{n}})\;\;\hbox{ and
}\;\;N:=xD\widehat{M}_{s,c}$

where $\nu_{i}=\min{(\mathop{v}(\widehat{M}_{s,c,i})-1,0)}$,
$\widehat{M}_{s,c,i}$ being the $i$th row of the matrix $\widehat{M}_{s,c}$.
By construction, one has $\mathop{v}{(N)}\geqslant 0$ and
$\mathop{v}{(D)}\geqslant 0$ so that

$D=\left(D_{0}+\mathop{\rm O}(x^{-1})\right)\;\;\hbox{ and
}\;\;N=N_{0}+\mathop{\rm O}(x^{-1}).$

Plugging the expansions of $D,N$, $Z=x^{d}(F_{0}+\mathop{\rm O}(x^{-1}))$ and
$x\delta(Z)=x^{d}(dF_{0}+\mathop{\rm O}(x^{-1}))$ in the Equation 12, one gets

$\displaystyle x^{d}\left((dD_{0}-N_{0})F_{0}+\mathop{\rm
O}(x^{-1})\right)=0.$

which implies

$(dD_{0}-N_{0})F_{0}=0$

and, since $F_{0}\neq 0$,

$\displaystyle\det\left(dD_{0}-N_{0}\right)=0.$ (13)

Thus the index $d$ must be a root of the polynomial
$E_{0}(\lambda):=\det\left(\lambda D_{0}-N_{0}\right)$.
In order to get a finite set of candidates for the d’s, it suffices to make
sure that the system $\delta(Z)=\widehat{M}_{s,c}Z$ is such that
$E_{0}(\lambda)$ is not identically zero (in which case we say that it is
$0$-simple). As mentioned before, it is always possible to reduce a given
difference system to a gauge-equivalent system which is simple [9].

The existing (and implemented) algorithms for computing simple forms [9]
produce a polynomial $E_{0}(\lambda)$ and a gauge transformation $T$ which is
polynomial in $x^{-1}$ and such that its inverse $T^{-1}$ is polynomial in
$x$.
Since $Tx^{d}F=x^{\tilde{d}}(\widetilde{F}_{0}+\mathop{\rm O}(x^{-1}))$, with
$\widetilde{F}_{0}\neq 0$ and $E_{0}(\tilde{d})=0$, it follows that
$d\in\tilde{d}+\mathop{v}(T)+\mathbb{N}$. Hence we only know $d$ modulo
$\mathbb{N}$ and this is sufficient for what we want to do!

Algorithm: Algorithm Unramified Generalized Exponents
Input: $M\in\mathop{\rm GL}_{n}({C}(x))$
Output: the set $G$ of the triplets $(s,c,d)$ involved in the lhs of $\tau
Y=MY$

  * 1.

$SC:=\textit{sc-Pairs}(M)$.

  * 2.

$G:=\emptyset$.

  * 3.

For $(s,c)\in SC$ do

    * –

Let $\widehat{M}_{s,c}:=c^{-1}x^{-s}M-I_{n}$.

    * –

Compute a $0$-simple form of $\widehat{M}_{s,c}$. This produces a polynomial
$E_{0}(\lambda)$ and a gauge transformation $T$ which is polynomial in $1/x$.

    * –

Compute the set $R_{sc}:=\\{d\in{C}\;|\;E_{0}(d)=0\\}$.

    * –

If $R_{sc}\neq\emptyset$ then
$G:=G\cup\left(\\{(s,c)\\}\times(R_{sc}+\mathop{v}(T))\right)$.

##  7 Application: Beke-Bronstein Algorithm

The Beke-Bronstein algorithm ([14] with improvements in [11]) for factoring
differential operators applies to polynomials and difference operators as
well. In the recurrence case, this leads to a recurrence system, and computing
hypergeometric solutions. In the following we reformulate the Beke-Bronstein
approach in the language of exterior algebra first for the polynomial case and
then the difference case (it works almost the same for both cases).

###  7.1 Polynomial case

Let $f\in\mathbb{Q}[y]$ be the polynomial that is to be factored. Assume
$y\nmid f$ throughout this section. If $g$ is a factor of $f$ then so is $ag$
where $a\in\mathbb{Q}-\\{0\\}$. To identify these factors we use
$\mathbb{P}(\mathbb{Q}[x])=(\mathbb{Q}[x]-\\{0\\})/(\mathbb{Q}-\\{0\\})$. If
$g\in\mathbb{Q}[x]-\\{0\\}$ then let $[g]$ denote its image in
$\mathbb{P}(\mathbb{Q}[x])$. Call $[g]$ a _projective factor_ of $f$ if $g\mid
f$. Let $\mathrm{Fact}_{m}(f)=\\{[g]:g\mid f,\deg(g)=m\\}$ be the set of
degree-$m$ projective factors.

Let $M=\mathbb{Q}[y]/(f)$. Then $M$ is a $\mathbb{Q}[y]$-module. For $g\mid f$
with $\deg(g)=m$, let $\mu(g)=g\wedge yg\wedge\cdots\wedge
y^{n-m-1}g\in\bigwedge^{n-m}M$ where $n$ is the degree of $f$. Suppose
$g=\sum_{i=0}^{m}b_{i}y^{i}$.

######  Lemma 22.

The coordinate of $\mu(g)$ with respect to $y^{i}\wedge y^{m+1}\wedge
y^{m+2}\wedge\cdots\wedge y^{n-1}$ is $b_{i}(b_{m})^{n-m-1}$.

###### Proof.

It is a straight-forward calculation. ∎

If $y\nmid f$ then $\mu(g)\neq 0$. For $a\in\mathbb{Q}^{*}$, we have
$\mu(ag)=a^{n-m}\mu(g)\neq 0$. Thus

$\mu([g]):=\mathbb{Q}\mu(g)\subset\bigwedge^{n-m}M$

is well-defined; $\mu([g])\in\mathbb{P}(\bigwedge^{n-m}M)$ does not depend on
the choice of representatives. 22 shows how to recover $[g]$ from $\mu([g])$.
This gives a map
$\eta:\mathbb{P}(\bigwedge^{n-m}M)\rightarrow\mathbb{P}(\mathbb{Q}[y])$ such
that $\eta\circ\mu$ is the identity on $\mathrm{Fact}_{m}(f)$.

######  Proposition 23.

Let

$E_{m}(f):=\\{\text{1-dimensional subspaces of }\bigwedge^{n-m}M\text{ that
are }\mathbb{Q}[y]\text{-modules}\\}.$

Then

$\mu(\mathrm{Fact}_{m}(f))\subseteq E_{m}(f).$

Now $\mathbb{Q}v$ is in $E_{m}(f)$ if and only if $v$ is an eigenvector of the
action of $y$ on $\bigwedge^{n-m}M$. 22 converts such $v$ to a polynomial $g$.
Then $g\in\mathrm{Fact}_{m}(f)$ (i.e. $g|f$) if and only if the so-called
Plücker relations hold, which in this case are equivalent to $(y^{i}g)\wedge
v=0$ for $i<n-m$. (We may omit $i=0$ since $g\wedge v=0$ will hold by
construction.) Indeed, if $g|f$, then the corresponding $v=\mu(g)$ is the
wedge product of $y^{i}g$ for $i<n-m$, and hence $(y^{i}g)\wedge v=0$.

######  Example 24.

Suppose $f=y^{4}+2y^{3}+3y^{2}+2y+2$. The goal is to illustrate Beke-Bronstein
by computing factors of $f$ of degree 2. In $M=\mathbb{Q}[y]/(f)$ the relation

$y^{4}=-2y^{3}-3y^{2}-2y-2$

holds. From that we compute the action of $y$ on a basis of $\bigwedge^{2}M$:

$\begin{split}{1\wedge y^{3}}&\mapsto y\wedge
y^{4}=y\wedge(-2y^{3}-3y^{2}-2y-2)\\\ {y\wedge y^{3}}&\mapsto y^{2}\wedge
y^{4}=y^{2}\wedge(-2y^{3}-3y^{2}-2y-2)\\\ {y^{2}\wedge y^{3}}&\mapsto
y^{3}\wedge y^{4}=y^{3}\wedge(-2y^{3}-3y^{2}-2y-2)\\\ 1\wedge y&\mapsto
y\wedge y^{2}\\\ 1\wedge y^{2}&\mapsto y\wedge y^{3}\\\ y\wedge y^{2}&\mapsto
y^{2}\wedge y^{3}\\\ \end{split}$

which can be represented by the matrix

$\begin{pmatrix}0&-2&0&2&0&-3\\\ 0&0&-2&0&2&2\\\ 2&2&3&0&0&0\\\ 0&0&0&0&1&0\\\
0&1&0&0&0&0\\\ 0&0&1&0&0&0\end{pmatrix}.$

The eigenvalues and eigenvectors are

$(2,\begin{pmatrix}-1\\\ 0\\\ -1\\\ -1\\\ 0\\\
1\end{pmatrix}),\quad(1,\begin{pmatrix}1\\\ 1\\\ 1/2\\\ 2\\\ 2\\\
1\end{pmatrix}).$

The Plücker relation(s) for $n=4$, $m=2$ are well known and are given by a
single polynomial

$X_{0,1}X_{2,3}-X_{0,2}X_{1,3}+X_{0,3}X_{1,2},$

where $X_{i,j}$ is the coordinate corresponding to the basis vector
$y^{i}\wedge y^{j}$. Both eigenvectors satisfy this Plücker relation. To
obtain the corresponding factors of $f$, we take the $y^{i}\wedge y^{3}$
coefficients of the above eigenvectors. We sorted the basis of
$\bigwedge^{2}M$ in such a way that these coefficients are the first three
entries of the eigenvectors. This leads to the factors

$-1-y^{2},\quad 1+y+\frac{1}{2}y^{2}.$

###  7.2 Difference case

Let $D=\mathbb{Q}(x)[\tau]$, let $L\in D$, and let $[R]=\mathbb{Q}(x)R$. We
call $[R]$ a _projective factor_ of $L$ if $R$ is a right-hand factor of $L$.
As before,

$\mathrm{Fact}_{m}(L):=\\{[R]:R\text{ is an order-}m\text{ right-hand factor
of }L\\}.$

Denote $M=D/DL$. Define

$\mu:\mathrm{Fact}_{m}(f)\rightarrow\bigwedge^{n-m}M,$
$[R]\mapsto\mathbb{Q}(x)\mu(R)$

where $\mu(R)=R\wedge\tau R\wedge\cdots\wedge\tau^{n-m-1}R$. Similar to 23 we
have

######  Proposition 25.

Let

$E_{m}(L):=\\{1\text{-dimensional }\mathbb{Q}(x)\text{-subspaces of
}\bigwedge^{n-m}M\text{ that are }D\text{-modules}\\}.$

Then $\mu(\mathrm{Fact}_{m}(L))\subseteq E_{m}(L)$.

As in the polynomial case, finding $\mathrm{Fact}_{m}(L)$ reduces to finding
1-dimensional submodules of $\bigwedge^{n-m}M$. Here this means finding
hypergeometric solutions.

######  Example 26.

Let $L=$

$(x^{3}+5x-6)\tau^{4}+(x^{4}+2x^{3}+6x^{2}-3x-18)\tau^{3}-(x^{4}+5x^{3}+14x^{2}+28x-12)\tau^{2}\\\
-(x^{5}+3x^{4}+8x^{3}+3x^{2}-21x-18)\tau+{3x^{2}(x^{2}+3x+8)}.$

To find order-$2$ factors of $L$, consider the $D$-module $\bigwedge^{2}D/DL$.
Like before, computing the action of $\tau$ on $\bigwedge^{2}D/DL$ with
respect to the basis

$e_{0}=1\wedge\tau^{3},e_{1}=\tau\wedge\tau^{3},e_{2}=\tau^{2}\wedge\tau^{3},e_{3}=1\wedge\tau,e_{4}=1\wedge\tau^{2},e_{5}=\tau\wedge\tau^{2}$

gives a matrix very similar to the one in the previous example. It has two
hypergeometric solutions, just like the previous example had two eigenvectors,
and again both satisfy the Plücker relation(s). From the first three entries
one finds:

$\tau^{2}-x,\quad\tau^{2}-x\tau-3.$

We have two implementations for this. The first implementation [1] computes
the hypergeometric solutions with Algorithm Version III.

The second factoring implementation computes generalized exponents, not of the
matrix, but of the operator $L$, and then constructs the generalized exponents
of the matrix from them. It also uses data computed from $L$ to determine, for
each individual entry: a bound for the denominator, a bound for the content of
the numerator, and a degree bound, similar to [18, Chapter 6]. This latter
implementation, written by van Hoeij, is available in the LREtools package in
Maple.

Both implementations are quite efficient, for example, we have successfully
tested the second implementation for factors of order $m=6$ of operators of
order $n=12$. This involves computing hypergeometric solutions of a system of
order $N$ by $N$, where $N={\rm binomial}(n,m)=924$.

## References

  * [1] Implementation.  http://www.math.fsu.edu/~hoeij/files/HypSols.
  * [2] Sergei Abramov.  Eg-eliminations.  Journal of Difference Equations and Applications, 5:393–433, 01 1999.
  * [3] Sergei A. Abramov and Moulay A. Barkatou.  Rational solutions of first order linear difference systems.  In Proceedings of ISSAC’98, Rostock, 1998.
  * [4] Sergei A. Abramov, Marko Petkovšek, and Anna A. Ryabenko.  Hypergeometric Solutions of First-Order Linear Difference Systems with Rational-Function Coefficients, pages 1–14.  Springer International Publishing, 2015.
  * [5] Moulay A. Barkatou.  Contribution à l’étude des équations différentielles et aux différences dans le champ complexe.  PhD thesis, l’Institut National Polytechnique de Grenoble, France, June 1989.
  * [6] Moulay A. Barkatou, Gary Broughton, and Eckhard Pflügel.  Regular systems of linear functional equations and applications.  In ISSAC ’08: Proceedings of the twenty-first international symposium on Symbolic and algebraic computation, pages 15–22, New York, NY, USA, 2008. ACM.
  * [7] Moulay A. Barkatou and G. Chen.  Computing the exponential part of a formal fundamental matrix solution of a linear difference system.  J. Difference Eq. Appl., 1996.
  * [8] Moulay A. Barkatou and G. Chen.  Some formal invariants of linear difference systems.  Pub. IRMA Lille, 44(2), 1997.
  * [9] Moulay A. Barkatou, Thomas Cluzeau, and Ali El-Hajj.  Simple forms and rational solutions of pseudo-linear systems.  In James H. Davenport, Dongming Wang, Manuel Kauers, and Russell J. Bradford, editors, Proceedings of the 2019 on International Symposium on Symbolic and Algebraic Computation, ISSAC 2019, Beijing, China, July 15-18, 2019, pages 26–33. ACM, 2019.
  * [10] Manuel Bronstein.  Personal communication.
  * [11] Manuel Bronstein.  An improved algorithm for factoring linear ordinary differential operators.  In Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC ’94, pages 336–340, New York, NY, USA, 1994. ACM.
  * [12] Yongjae Cha.  Closed Form Solutions of Linear Difference Equations.  PhD thesis, Florida State University, 2010.
  * [13] Thomas Cluzeau and Mark van Hoeij.  Computing hypergeometric solutions of linear recurrence equations.  Appl. Algebra Eng. Commun. Comput., pages 83–115, 2006.
  * [14] Beke E.  Die Irreduzibilität der homogenen linearen Differentialgleichungen.  Math. Ann., 45:278–294, 1894.
  * [15] M. van der Put and M.F. Singer.  Galois Theory of Linear Differential Equations.  Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 2003.
  * [16] Mark van Hoeij.  Finite singularities and hypergeometric solutions of linear recurrence equations.  Journal of Pure and Applied Algebra, 139(1):109 – 131, 1999.
  * [17] Robert J. Walker.  Algebraic Curves.  Springer-Verlag New York, 1978.
  * [18] Yi Zhou.  Algorithms for Factoring Linear Recurrence Relations.  PhD thesis, Florida State University, USA, April 2022.