Computing Gröbner and Involutive Bases for Linear Systems of Difference Equations

. The computation of involutive bases and Gröbner bases for linear systems of di ﬀ erence equations is solved and its importance for physical and mathematical problems is discussed. The algorithm and issues concerning its implementation in C are presented and calculation times are compared with the competing programs. The paper ends with consideration on the parallel version of this implementation and its scalability.


Introduction
Recently there was considerable interest in applying difference Gröbner bases to the reduction of multiloop Feynman integrals to minimal sets of basis integrals [1]. The related systems of linear difference polynomial equations are difficult to transform into the Gröbner basis form with available software. In this note we present a competitive software implementing the involutive approach to compute the basis form.

Elements of Difference Algebra [2]
Let {y 1 , . . . , y m } be the set of indeterminatesm functions (dependent variables) of n variables x 1 , . . . , x n (independent variables), and {θ 1 , . . . , θ n } be the set of mutually commuting difference operators: A difference ring R with differences θ 1 , . . . , θ n is a commutative ring: Let R be the difference ring of polynomials over field K in variables Let R L be the set of linear polynomials from R and let e-mail: yan@jinr.ru A difference ideal is an ideal J ⊆ R closed under the action of any operator from Θ. For F ∈ R, the smallest difference ideal containing F will be denoted by Id(F). If for an ideal J there is F ∈ R L such that J =Id(F), then J is a linear difference ideal.
A total ordering on the set of θ µ • y j is a ranking if for all i, j, k, µ, ν the following holds: If |µ| > |ν| implies θ µ • y j θ ν • y k for all j, k, then the ranking is orderly. If j > k implies θ µ • y j θ ν • y k for all µ, ν, then the ranking is elimination.
Given a ranking , a linear polynomial f ∈ R L \ {0} has the leading term aϑ A ranking acts in R L as a monomial order. If F ⊆ R L \ {0}, lm(F) will denote the set of the leading monomials and lm j (F) will denote its subset for the indeterminate y j . Thus,

Gröbner and Involutive Bases of a Difference Ideal [2-4]
Given a nonzero linear difference ideal J = Id(G) and a ranking , the ideal generating the set It follows that f ∈ J \ {0} is reducible modulo G: If f 0, then it is again reducible modulo G, and, by repeating the reduction, in finitely many steps we obtain Similarly, a nonzero polynomial h ∈ R L , the terms of which are reducible (if any) modulo a set F ⊂ R L , can be reduced to an irreducible polynomialh, which is said to be in normal form modulo F (denotation:h = NF(h, F)).
Replacing the conventional division with the involutive one [3] are obtained definitions of the involutive basis and involutive normal form which have several major advantages as compared to the Gröbner basis computation: automatic avoidance of some unnecessary reductions, smoother arithmetics growth during calculation and natural parallelism.
The author used the Janet division as involutive, but there exist theoretically better divisions, for example Janet-like [2].

Algorithm
In the Algorithm 1 below the computations at the lines 3-8, which consume about 90% of programm runtime are mutually independent, so they may be started simultaneously. This is a source of natural parallelism of the involutive algorithm [5][6][7].

Implementation Issues
This algorithm was implemented as a program written in C language. The initial system of linear difference equations rewrites the input data into the form of a polynomial system with parametric coefficients (independent variables) and power products of the forward shift operators θ µ , |µ| ≥ 0 applied to the dependent variable (sec. 1.1). For computations of the greatest common divisor for reduced polynomials with parametric coefficients our implementation uses the well-known multimodular Brown algorithm [8], arbitrary precision arithmetics uses libgmp.

Benchmarking
Due to their high practical value for the computation of Feynman integrals [1] we shall take some systems of difference equations as benchmarks and examples. As a competitive software, the Maple package LDA [2] was chosen. The comparative timings of sequental programs are shown in Table 1. The most difficult computed example, Smirnov3, was used as a benchmark for the parallel version of our program (Table 2). Smirnov1 [1]: Integral relations for the one-loop massless propagator − a 2 f (a 1 − 1, a 2 + 1) + a 2 q 2 f (a 1 , a 2 + 1).

Conclusion
Using a similar algorithmic background, our implementation demonstrates significantly better timings in comparison with the Maple implementation. It has thus the potential for computing relatively large problems. The parallel version shows rather good CPU-utilization, but it is not the final result. Hope to achieve better results in the future.