Automated Generation of Non-Linear Loop Invariants Utilizing Hypergeometric Sequences

TL;DR

This paper introduces a novel method for automatically generating all polynomial invariants for loops with complex arithmetic, including multiplication, by transforming programs into recurrence relations and applying algebraic techniques.

Contribution

The paper presents a new approach that extends the class of P-solvable loops handled by invariant generation, incorporating multiplication with loop counters and using Gr"obner basis computation.

Findings

Successfully computes polynomial invariants for complex loops

Extends the class of loops for which invariants can be automatically generated

Implemented in the Aligator Mathematica package

Abstract

Analyzing and reasoning about safety properties of software systems becomes an especially challenging task for programs with complex flow and, in particular, with loops or recursion. For such programs one needs additional information, for example in the form of loop invariants, expressing properties to hold at intermediate program points. In this paper we study program loops with non-trivial arithmetic, implementing addition and multiplication among numeric program variables. We present a new approach for automatically generating all polynomial invariants of a class of such programs. Our approach turns programs into linear ordinary recurrence equations and computes closed form solutions of these equations. These closed forms express the most precise inductive property, and hence invariant. We apply Gr\"obner basis computation to obtain a basis of the polynomial invariant ideal, yielding thus a finite representation of all polynomial invariants. Our work significantly extends the class of so-called P-solvable loops by handling multiplication with the loop counter variable. We implemented our method in the Mathematica package Aligator and showcase the practical use of our approach.

Equations

$$Sa=\sigma (a)S\text{for all }a\in \text{mathbbm}K(x),$$

$$s\sim t:\iff s-t\text{ has finitely many non-zero elements}.$$

$$\tau :\text{mathbbm}K(x)[S;\sigma ,0]\times \mathcal{S}\to \mathcal{S}$$

$$\displaystyle\tau(L(S,x),t)(n)=\tau\biggl{(}\sum_{i=0}^{d}l_{i}(x)S^{i},t\biggr{)}(n):=\sum_{i=0}^{d}l_{i}(n)t(n+i),$$

$$L_1=S-\frac{p(x+1)}{p(x)}.$$

$$f(n)=\frac{1}{\sqrt{5}}\left({\left(\frac{1+\sqrt{5}}{2}\right)}^n-{\left(\frac{1-\sqrt{5}}{2}\right)}^n\right),$$

$$h(n)=\sum _{i=1}^n\frac{1}{i},$$

$$c(n)=p_1(n){\theta }_1^n+p_2(n){\theta }_2^n+\cdots +p_s(n){\theta }_s^n,$$

$$h(n)={\theta }^{n}r(n)((n+{\zeta }_1{)}^{\underline{n}}{)}^{k_1}((n+{\zeta }_2{)}^{\underline{n}}{)}^{k_2}\cdots ((n+{\zeta }_{\ell }{)}^{\underline{n}}{)}^{k_{\ell }},$$

$$pred(v_1,\dots ,v_m)$$

$$I({\tilde{v}}_1,\dots ,{\tilde{v}}_k)=I(v_1,\dots ,v_m)\cap \text{mathbbm}K[{\tilde{v}}_1,\dots ,{\tilde{v}}_k].$$

$$\text{mathbbm}K=\text{mathbbm}F(v_{1,0},\dots ,v_{1,k},v_{2,0},\dots ,v_{m,\ell }),$$

$$I(quo,rem,x,y)=⟨rem+quo\cdot y-y\cdot quo(0)-rem(0)⟩.$$

$$\{\begin{array}{l}a(n+2)-5(n+2)\cdot a(n+1)-6(n^2+3n+2)\cdot a(n)=0\\ b(n+1)-2\cdot b(n)=0\\ c(n+1)-3(n+1)\cdot c(n)=0\\ d(n+1)-(n+1)\cdot d(n)=0.\end{array}$$

$$v_i(n)=\sum _{k\in \text{mathbbm}{Z}^{\ell }}{p}_{i,k}(n,{\theta }_1^n,\dots ,{\theta }_s^n)((n+{\zeta }_1{)}^{\underline{n}}{)}^{k_1}\cdots ((n+{\zeta }_{\ell }{)}^{\underline{n}}{)}^{k_{\ell }}$$

$$v(n)=h_1(n)+\cdots +h_w(n),\text{ where}$$

$$h_i(n)=q_i(n){\tilde{\theta }}_i^n((n+{\zeta }_1{)}^{\underline{n}}{)}^{k_{i,1}}\cdots ((n+{\zeta }_{\ell }{)}^{\underline{n}}{)}^{k_{i,\ell }},$$

$$G=$$

$$=$$

$$⋮$$

$$=$$

$$0=G(v)=G(h_1+\cdots +h_w)=G(h_1)+\cdots +G(h_w)=r_1{h}_1+\dots r_w{h}_w.$$

$$h_1=(-1{)}^{n}n!,h_2=6^{n}n!$$

$$a(n)=(k_1(-1{)}^n+k_2{6}^n)n!$$

$$b(n)=2^{n}b(0),c(n)=3^{n}n!c(0),d(n)=n!d(0).$$

$$\sum _{i=1}^w{h}_i(n+1)=r_h(n)\sum _{i=1}^w{h}_i(n).$$

$$\sum _{i=1}^w(r_i(n)-r_h(n))h_i(n)=0.$$

$$\sum _{i=1}^{w-1}(r_i(n)-r_h(n))h_i(n)=(r_w(n)-r_h(n))h_w(n).$$

$${\theta }_1^n\cdot ((n-{\zeta }_1{)}^{\underline{n}}{)}^{k_1}\cdots ((n-{\zeta }_{\ell }{)}^{\underline{n}}{)}^{k_{\ell }}=$$

$$r(n)\cdot {\theta }_2^n\cdot ((n-{\zeta }_1{)}^{\underline{n}}{)}^{c_1}\cdots ((n-{\zeta }_{\ell }{)}^{\underline{n}}{)}^{c_{\ell }},$$

$$\underset{\text{hypergeometric}}{\underbrace{{\left(\frac{{\theta }_1}{{\theta }_2}\right)}^n((n-{\zeta }_1{)}^{\underline{n}}{)}^{k_1-c_1}\cdots ((n-{\zeta }_{\ell }{)}^{\underline{n}}{)}^{k_{\ell }-c_{\ell }}}}=r(n).$$

$$q(x)=\frac{g(x)f(x+1)}{g(x+1)f(x)},$$

$$(n+{\zeta }_i{)}^{\underline{n}}\cdot \prod _{w=1}^k({\zeta }_j-w)=(n+{\zeta }_j{)}^{\underline{n}}\cdot \prod _{w=1}^k(n+w+{\zeta }_j).$$

Full text

program analysis, loop invariants, recurrence relations, hypergeometric sequences

††copyright: rightsretained††conference: The 42nd International Symposium on Symbolic and Algebraic Computation; July 2017; Kaiserslautern, Rheinland-Pfalz, Germany††journalyear: 2017††ccs: Theory of computation Invariants††ccs: Theory of computation Automated reasoning††ccs: Theory of computation Program verification††ccs: Mathematics of computing Discrete mathematics

1. Introduction

1.1. Overview

Analysis and verification of software systems requires non-trivial automation. Automatic generation of program properties describing safety and/or liveness is a key step to such automation, in particular in the presence of program loops (or recursion). For programs with loops one needs additional information, in the form of loop invariants or conditions on ranking functions.

In this paper we focus on loop invariant generation for programs with assignments implementing numeric computations over scalar variables. Our programming model extends the class of so-called P-solvable loops. Our work is based on and extends results of (Rodríguez-Carbonell and Kapur, 2007; Kovács, 2007), in particular it relies on the fact that the set of polynomial invariants of P-solvable loops form a polynomial ideal and we employ reasoning about C-finite and hypergeometric sequences to determine algebraic dependencies. We show how to compute the ideal of polynomial invariants of extended P-solvable loops as follows: we model programs as a system of recurrence equations and compute closed form sequence solutions of these recurrences. If these sequences are of a certain type, which includes, among others, polynomials, rational functions, exponential and factorial sequences, then we compute a set of generators of the polynomial invariant ideal via Gröbner bases. We implemented our approach in the Mathematica package Aligator (Kovács, 2008) that is able to compute polynomial loop invariants for programs that, to the best of our knowledge, no other approach is able to handle.

This paper is organized as follows. In Section 2, we state basic definitions and facts about the algebra of linear ordinary recurrence operators as well as C-finite and hypergeometric sequences. We also give a precise definition of the programming model we take into consideration, particularly the notion of imperative loops with assignment statements only. This is followed by a description of the class of P-solvable loops and its reach and limitations in Section 3. In Section 4 we present our main contribution, an extension of P-solvable loops by reasoning about hypergeometric sequences and we derive the necessary theoretical and algorithmical results to offer fully automated polynomial invariant generation therein. We conclude the paper with a presentation of our implementation in the Mathematica package Aligator in Section 5 and a summary of possible future research directions in Section 6.

1.2. Related Work

Many classical data flow analysis problems, such as constant propagation and finding definite equalities among program variables, can be seen as problems about polynomial identities expressing loop invariants. In (Müller-Olm and Seidl, 2004; Sankaranarayanan et al., 2004) a method built upon linear and polynomial algebra is developed for computing polynomial equalities of a bounded degree. A related approach was also proposed by (Rodriguez-Carbonell and Kapur, 2007) using abstract interpretation. Abstract interpretation is also used in (Farzan and Kincaid, 2015; de Oliveira et al., 2016) for computing polynomial invariants of programs whose assignments can be described by C-finite recurrences. In our work we do not rely on abstract interpretation but use algebraic reasoning about holonomic sequences. For program loops with assignments only, our technique can handle programs with more complex arithmetic than the previously mentioned methods. Our work is currently restricted though to single-path loops.

Without an a priori fixed polynomial degree, in (Rodríguez-Carbonell and Kapur, 2007) the polynomial invariant ideal is approximated by a fixed point procedure based on polynomial algebra and abstract interpretation. In (Kovács, 2007), the author defines the notion of P-solvable loops which strictly generalizes the programming model of (Rodríguez-Carbonell and Kapur, 2007). Given a P-solvable loop with assignments and nested conditionals, the results in (Kovács, 2007) yield an automatic approach for computing all polynomial loop invariants. Our work extends (Kovács, 2007; Rodríguez-Carbonell and Kapur, 2007) in new ways: it handles a richer class of P-solvable loops where multiplication with the loop counter is allowed. Our technique relies on manipulating hypergeometric sequences and relaxes the algebraic restrictions of (Kovács, 2007; Rodríguez-Carbonell and Kapur, 2007) on program operations. To the best of our knowledge, no other method is able to derive polynomial invariants for extended P-solvable loops. Unlike (Kovács, 2007; Rodríguez-Carbonell and Kapur, 2007), we however only treat loops with assignments; that is, invariants for extended P-solvable loops with conditionals are not yet treated by our approach.

2. Preliminaries

In this section we give a brief overview of the algebra of linear ordinary recurrence operators as well as C-finite and hypergeometric sequences that we use further on. We also describe our programming model in detail.

2.1. Recurrence Operators and Holonomic Sequences

Let $\mathbbm{K}$ be a computable field of characteristic zero.

The algebra of linear ordinary recurrence operators in one variable will serve as the algebraic foundation to deal with recurrence equations. For details on general Ore algebras, see (Bronstein and Petkovšek, 1996; Ore, 1933).

Definition 2.1.

Let $\mathbbm{K}(x)[S]$ be the set of univariate polynomials in the variable $S$ over the set of rational functions $\mathbbm{K}(x)$ in $x$ and let $\sigma\colon\mathbbm{K}(x)\rightarrow\mathbbm{K}(x)$ be the forward shift operator in $x$, i.e. $\sigma(r(x))=r(x+1)$ for $r(x)\in\mathbbm{K}(x)$. We define the Ore polynomial ring of ordinary recurrence operators $(\mathbbm{K}(x)[S],+,\cdot)$ with component-wise addition and the unique distributive and associative extension of the multiplication rule

[TABLE]

to arbitrary polynomials in $\mathbbm{K}(x)[S]$. To clearly distinguish this ring from the commutative polynomial ring over $\mathbbm{K}(x)$, we denote it by $\mathbbm{K}(x)[S;\sigma,0]$. The order of an operator $L\in\mathbbm{K}(x)[S;\sigma,0]$ is its degree in $S$.

Without loss of generality, we assume that the leading coefficient of any operator $L\in\mathbbm{K}(x)[S;\sigma,0]$ is equal to 1. Otherwise, we can divide by the leading coefficient of $L$ from the left. $\mathbbm{K}(x)[S;\sigma,0]$ is a right Euclidean domain, i.e. we have the notion of the greatest common right divisor and the least common left multiple of operators and we are able to determine both algorithmically. Consequently, $\mathbbm{K}(x)[S;\sigma,0]$ is a principal left ideal domain and every left ideal is generated by the greatest common right divisor of a given set of generators.

Consider the ring $\mathbbm{K}^{\mathbbm{N}}$ of all sequences in $\mathbbm{K}$ with component-wise addition and the Hadamard product (i.e. component-wise product) as multiplication. We follow (Petkovšek et al., 1996) in identifying sequences as equal if they only differ in finitely many terms. This will prove beneficial in two ways. Firstly, it allows us to define the action of operators on sequences in a natural way. Secondly, disregarding finitely many starting values makes it possible to identify unnecessary loop variables, whose values are eventually equal to the values of another variable, and therefore can be computed outside of any while loop. Let $\sim$ be the equivalence relation on $\mathbbm{K}^{\mathbbm{N}}$ defined by

[TABLE]

We then set $\mathcal{S}$ to be the quotient ring $\mathbbm{K}^{\mathbbm{N}}/\sim$. Subsequently, it will not be necessary to distinguish between $t\in\mathbbm{K}^{\mathbbm{N}}$ and $\pi(t)\in\mathcal{S}$, where $\pi:\mathbbm{K}^{\mathbbm{N}}\rightarrow\mathcal{S}$ is the canonical homomorphism. The field $\mathbbm{K}$ can be embedded in $\mathcal{S}$ via the map $c\mapsto(c)_{n\in\mathbbm{N}}$. The action of an operator in $\mathbbm{K}(x)[S;\sigma,0]$ on an element in $\mathcal{S}$ is defined by the map

[TABLE]

where the evaluation is well defined for all $n\geq n_{0}$ for some $n_{0}\in\mathbbm{N}$, and we set $L(t):=\tau(L,t)\in\mathcal{S}$. If $L(t)\equiv 0$, then we say that $L$ is an annihilator of $t$ ($L$ annihilates $t$) and $t$ is a solution of $L(t)=0$. A sequence that is annihilated by a non-zero operator in $\mathbbm{K}(x)[S;\sigma,0]$ is called holonomic sequence. For a given sequence $t$, the set of all its annihilators forms a left ideal in $\mathbbm{K}(x)[S;\sigma,0]$. We call it the annihilator ideal of $t$ and denote it by $\operatorname{ann}(t)$.

Example 2.2.

Let $p(x)$ be a polynomial in $\mathbbm{K}[x]$. The polynomial sequence $(p(n))_{n\in\mathbbm{N}}$ is annihilated by the operator

[TABLE]

$L_{1}$ is a generator of the annihilator ideal of $p$. Set $\Delta:=S-1$. Then $\tilde{p}=\Delta(p)$ is again a polynomial sequence with $\deg(\tilde{p})<\deg(p)$. It follows that $L_{2}=\Delta^{\deg(p)+1}$ is another annihilator of $p$ in $\mathbbm{K}(x)[S;\sigma,0]$ and its coefficients are independent of $x$. Since $L_{1}$ generates $\operatorname{ann}(p)$, there exists an operator $Q$ with $L_{2}=QL_{1}$.

In our work, we focus on two different special kinds of holonomic sequences:

Definition 2.3.

Let $t\in\mathcal{S}$. Then

• •

$t$ is called C-finite if it is annihilated by an operator in $\mathbbm{K}(x)[S;\sigma,0]$ with only constant coefficients. $(l_{i}\in\mathbbm{K}$)

• •

$t$ is called hypergeometric if it is annihilated by an order 1 operator in $\mathbbm{K}(x)[S;\sigma,0]$.

Example 2.4.

We give some examples of commonly encountered sequences.

• •

As was shown in Example 2.2, polynomial sequences are both, C-finite and hypergeometric.

• •

Rational function sequences $(r(n))_{n\in\mathbbm{N}}$, $r\in\mathbbm{K}(x)\setminus\mathbbm{K}[x]$, are hypergeometric but not C-finite.

• •

The factorial sequence $(n!)_{n\in\mathbbm{N}}$ is hypergeometric but not C-finite.

• •

The Fibonacci sequence $(f(n))_{n\in\mathbbm{N}}$ with

[TABLE]

is C-finite but not hypergeometric.

• •

The sequence of harmonic numbers $(h(n))_{n\in\mathbbm{N}}$ with

[TABLE]

is neither hypergeometric nor C-finite.

In a sufficiently large algebraic field extension $\overline{\mathbbm{K}}/\mathbbm{K}$, every C-finite sequence $(c(n))_{n\in\mathbbm{N}}$ can be uniquely written (up to reordering) in the form

[TABLE]

for some $s\in\mathbbm{N}$ and $p_{i}\in\mathbbm{K}[x],\theta_{i}\in\overline{\mathbbm{K}}$ for $i=1,\dots,s$ with $\theta_{i}\neq\theta_{j}$ for $i\neq j$. For any $r\in\mathbbm{K}(x)$ and $n\in\mathbbm{N}$, $r(x)^{\underline{n}}$ is defined as $\prod_{i=0}^{n-1}r(x-i)$. Then every hypergeometric sequence $(h(n))_{n\in\mathbbm{N}}$ can be uniquely written (up to reordering) in the form

[TABLE]

for some $\ell\in\mathbbm{N}$, $r(x)\in\mathbbm{K}(x)$, $\theta\in\overline{\mathbbm{K}}$, $\zeta_{i}\in\overline{\mathbbm{K}}$ and $k_{i}\in\mathbbm{Z}$ for $i=1,\dots,\ell$, and the difference $\zeta_{i}-\zeta_{j}$ is not an integer for $i\neq j$. From these closed forms it is immediate that finite sums and products of C-finite sequences are again C-finite and finite products of hypergeometric sequences are again hypergeometric. Sums of hypergeometric sequences are not necessarily hypergeometric, see Lemma 4.3. Subsequently, we will assume that $\mathbbm{K}$ is large enough so that all occurring C-finite and hypergeometric sequences have a closed form representation in $\mathbbm{K}$.

For more details on C-finite and hypergeometric sequences, as well as proofs for the facts given in this section, see (Kauers and Paule, 2011).

For functions $f_{1},\dots,f_{m}:\mathbbm{U}\rightarrow\mathbbm{K}$ with $\mathbbm{N}\subset\mathbbm{U}\subset\mathbbm{K}$ that are algebraically independent over $\mathbbm{K}$, we distinguish between the polynomial ring $\mathbbm{K}[f_{1},\dots,f_{m}]$, where $f_{1},\dots,f_{m}$ are used as variables, and the ring $\mathbbm{K}[f_{1}(n),\dots,f_{m}(n)]\subset\mathcal{S}$ of all sequences $(t(n))_{n\in\mathbbm{N}}$ of the form $t(n)=p(f_{1}(n),\dots,f_{m}(n))$ with $p\in\mathbbm{K}[f_{1},\dots,f_{m}]$. This distinction is important, as e.g. the function $\sin(x\cdot\pi)$ is algebraically independent over $\mathbbm{K}$, but the sequence $(\sin(n\cdot\pi))_{n\in\mathbbm{N}}=(0,0,0,\dots)$ is not, and thus $\mathbbm{K}[\sin(n\cdot\pi)]$ is isomorphic to $\mathbbm{K}$, but $\mathbbm{K}[\sin(x\cdot\pi)]$ is not.

Remark 0.

In the context of this paper, since the operators in question emerge from program loops, we can safely assume that the rational function coefficients of any operator do not have poles in $\mathbbm{N}$. Otherwise, a division by zero error would occur for some program input.

2.2. Programming Model

We consider a simple programming model of single-path loops with rational function assignments. That is, nested loops and/or loops with conditionals are not yet handled in our work. Our programming model is thus given by the following loop pattern, written in a C-like syntax:

[TABLE]

where $v_{1},\dots,v_{m}$ are (scalar) variables with values from $\mathbbm{K}$, the $f_{i}$ are rational functions over $\mathbbm{K}$ in $m$ variables and $pred$ is a a Boolean formula (loop condition) over $v_{1},\dots,v_{m}$. In our approach however we ignore loop conditions and treat program loops as non-deterministic programs. In (Müller-Olm and Seidl, 2004), it is shown that the set of all affine equality invariants is not computable if the programming model includes affine equality tests/conditions. With this consideration, our programming model from (1) becomes:

[TABLE]

Due to particular importance in our reasoning, we suppose that there is always a variable $n$ denoting the loop iteration counter. The initial value of $n$ will always be $n=0$ and $n$ will be incremented by $1$ at the end of each iteration.

Each program variable gives rise to a sequence $(v_{i}(n))_{n\in\mathbbm{N}}$. For a program variable $v$, we allow ourselves to abuse the notation and also use the identifier $v$ as a variable in polynomial rings as well as an identifier for the sequence $(v(n))_{n\in\mathbbm{N}}$.

A polynomial loop invariant is a non-zero polynomial $p$ over $\mathbbm{K}$ in $m$ variables such that $p(v_{1}(n),\dots,v_{m}(n))=0$ for all $n$. As observed in (Rodríguez-Carbonell and Kapur, 2007; Kovács, 2007), the set of all polynomial invariants forms a polynomial ideal in $\mathbbm{K}[v_{1},\dots,v_{m}]$, called the polynomial invariant ideal and is denoted by $I(v_{1},\dots,v_{m})$. For a subset $\{\tilde{v}_{1},\dots,\tilde{v}_{k}\}\subset\{v_{1},\dots,v_{m}\}$, we define

[TABLE]

In general, polynomial loop invariants depend on the initial values of program variables. To simplify the presentation, we fix $\mathbbm{K}$ to be

[TABLE]

for a computable field $\mathbbm{F}$ of characteristic zero that allows us to represent all occurring C-finite and hypergeometric sequences in closed form, and sufficiently many variables $v_{1,0},\dots,v_{m,\ell}$ that represent the initial values of the program variables $v_{1},\dots,v_{m}$.

3. Polynomial Invariants for P-Solvable Loops

We now turn our attention to the class of P-solvable loops introduced in (Kovács, 2007) that allows for computing all polynomial loop invariants..

Definition 3.1.

An imperative loop with assignment statements only is called P-solvable if the sequence of each recursively changed program variable $v$ is C-finite and the ideal of all polynomial invariants over $\mathbbm{K}$ is not the zero ideal.

Example 3.2.

In (Kovács, 2007), it is shown that the Euclidean algorithm is P-solvable. Given the program:

[TABLE]

The ideal of polynomial loop invariants is shown to be

[TABLE]

With $quo(0)=0$ and $rem(0)=x$, this gives $\langle rem+quo\cdot y-x\rangle$.

While P-solvable loops cover a wide class of program loops, there are several significant cases which do not fall into this class. Notably, multiplication with the loop counter $n$ will generally result in loops that are not P-solvable.

Example 3.3.

Consider the following loop with relevant loop variables $a,b,c,d$. The variables $t_{1},t_{2}$ are temporary variables used to access previous values of $a$. Along with the loop counter $n$, we will not take them into consideration for the loop invariants in this example.

[TABLE]

The program then satisfies the following system of recurrences:

[TABLE]

This loop is not P-solvable as, for example, the variable $c$ is updated by a sequence that is not C-finite (due to the multiplication between the program variables $n$ and $c$). To the best of our knowledge, none of the existing invariant generation techniques is able to to compute polynomial invariants for this loop. In the next section, we extend the class of P-solvable loops, covering also programs as the one above, and introduce an automated approach to derive all polynomial invariants of such loops.

4. Extension of P-Solvable Loops

4.1. Definition of Extended P-Solvable Loops

Consider the sequences $(v_{1}(n))_{n\in\mathbbm{N}},\dots,(v_{m}(n))_{n\in\mathbbm{N}}$ with values in $\mathbbm{K}$ given by

[TABLE]

where $s,\ell\in\mathbbm{N}$, the $p_{i,k}$ are polynomials in $\mathbbm{K}(x)[y_{1},\dots,y_{s}]$, not identically zero for finitely many $k\in\mathbbm{Z}^{\ell}$, and the $\theta_{i}$ and $\zeta_{j}$ are elements of $\mathbbm{K}$ for $i=1,\dots,s$, $j=1,\dots,\ell$ with $\theta_{i}\neq\theta_{j}$ and $\zeta_{i}-\zeta_{j}\notin\mathbbm{Z}$ for $i\neq j$.

In particular, this class of sequences comprises C-finite sequences as well as hypergeometric sequences and Hadamard products of C-finite and hypergeometric sequences, which could not be handled in automated invariant generation before. We give an extension of Definition 3.1 based on this class of sequences

Definition 4.1.

An imperative loop with assignment statements only is called extended P-solvable if the sequence of each recursively changed program variable $v$ is of the form (3).

Note that in Definition 4.1, we drop the requirement of Definition 3.1 that the ideal of algebraic relations is not the zero ideal. This change is just for convenience.

While it is obvious that the inclusion of hypergeometric terms in extended P-solvable loops allows assignments of the form $v\mathrel{:=}r(n)v$, where $r$ is a rational function in $\mathbbm{K}[x]$, it also allows assignments that turn into higher order recurrences, as illustrated in Example 4.2. It also allows for assignments of the form $v_{2}\mathrel{:=}r(v_{1})v_{2}$, with $r\in\mathbbm{K}(x)$, as long as the closed form of $v_{1}$ is a rational function in $n$.

4.2. Detecting Extended P-Solvable Loops

In order to employ the ideas we develop in Section 4.3 for finding algebraic relations in extended P-solvable loops, we have to be able to detect sequences of the form (3). This means, given a recurrence operator $R$ of order $d$ and starting values $s_{0},\dots,s_{d-1}$, compute, if possible, $p_{k}$,$\theta_{i}$ and $\zeta_{j}$ as in (3) such that $v$ is a solution of $R(v)=0$ with $v(n)=s_{n}$ for $n\in\{0,\dots,d-1\}$. We can write $v$ as a sum of hypergeometric sequences:

[TABLE]

with $q_{i}\in\mathbbm{K}(x)$, $\tilde{\theta}_{i}\in\mathbbm{K}$, and $k_{i}\in\mathbbm{Z}^{\ell}$. Note that we use $\tilde{\theta}_{i}$ instead of $\theta_{i}$ since the exponential sequence for each summand can be a product of several $\theta_{i}^{n}$. We can assume without loss of generality that the $h_{i}$ are linearly independent over $\mathbbm{K}(n)$. In fact, if $h_{1}(n)=r_{2}(n)h_{2}(n)+\dots+r_{w}(n)h_{w}(n)$, we can set $\tilde{h}_{1}=(1+r_{2})h_{2},\dots,\tilde{h}_{w-1}=(1+r_{w})h_{w}$ and get $v(n)=\tilde{h}_{1}(n)+\dots+\tilde{h}_{w-1}(n)$. Let $L$ be the least common left multiple of the first order operators $L_{1},\dots,L_{w}$ that annihilate $h_{1},\dots,h_{w}$ respectively in the Ore algebra $\mathbbm{K}(x)[S;\sigma,0]$ and let $G$ be a generator of $\operatorname{ann}(v)$. We show that $G$ and $L$ are equal. (Note that we required all operators to have leading coefficient $1$.)

By right division with remainder, we can write $G$ as

[TABLE]

with $Q_{1},\dots,Q_{w}\in\mathbbm{K}(x)[S;\sigma,0]$ and some $r_{1},\dots,r_{w}\in\mathbbm{K}(x)$. We then get

[TABLE]

Since the $h_{i}$ are linearly independent, we have $r_{1}=\dots=r_{w}=0$, and so, $L_{1},\dots,L_{w}$ are right factors of $G$. This proves the claim.

Since every annihilator of $v$ is a multiple of $G$ and therefore also an annihilator of $h_{i}$, we can use Petkovšek’s algorithm (Petkovšek, 1992) to determine $p_{k}$,$\theta_{i}$ and $\zeta_{j}$ as in (3). More precisely, given an operator $R\in\mathbbm{K}(x)[S;\sigma,0]$ of order $d$ and starting values $s_{0},\dots,s_{d-1}$, we compute $v$ as in (3) such that $R(v)=0$ (if possible), by computing all hypergeometric solutions of $R$. This gives $\theta_{i},\zeta_{i}$ and $p_{i}$, linearly dependent on parameters $c_{1},\dots,c_{w}$. Next, we solve the linear system $v(i)=s_{i}$ in terms of $c_{i}$. Any solution then gives rise to a sequence $(v(n))_{n\in\mathbbm{N}}$ with the desired properties.

Example 4.2.

For the recurrence for $a$ in Example 3.3, we compute two hypergeometric solutions using Petkovšek’s algorithm:

[TABLE]

Thus, we get

[TABLE]

with the relations $a(0)=k_{1}+k_{2}$ and $a(1)=6k_{2}-k_{1}$ stemming from the starting values of $a$. Since $b,c,d$ are given by first order recurrences, their closed forms can be easily computed:

[TABLE]

It follows that the program loop given in Example 3.3 is extended P-solvable.

4.3. The Ideal of Algebraic Relations

We now turn to the problem of, given sequences $v_{1},\dots,v_{m}$ as in (3), how to compute a basis for the ideal $I(v_{1},\dots,v_{m})$ of all algebraic relations among the $v_{i}$. We proceed by identifying the terms $(n+\zeta_{i})^{\underline{n}}$ that are algebraically independent over $\mathbbm{K}(n,\theta_{1}^{n},\dots,\theta_{s}^{n})$. For this, we use basic properties of sums and products of hypergeometric terms. First, we state a necessary condition for a finite sum of hypergeometric terms to be again hypergeometric.

Lemma 4.3.

Let $h_{1},\dots,h_{w}$ be hypergeometric sequences. If the sum $h_{1}+\dots+h_{w}$ is hypergeometric, then there exist integers $i,j\in\{1,\dots,w\}$, $i\neq j$, and a rational function $r(x)\in\mathbbm{K}(x)$ such that $h_{i}(n)=r(n)h_{j}(n)$.

Proof.

We prove the claim by induction on $w$. For the case $w=1$, there is nothing to show. Now suppose the claim holds for some $(w-1)\in\mathbbm{N}^{*}$. There is a rational function $r_{h}(x)\in\mathbbm{K}(x)$ such that

[TABLE]

Let $r_{i}\in\mathbbm{K}(x)$ be such that $h_{i}(n+1)=r_{i}(n)h_{i}(n)$. We then get

[TABLE]

We first treat the case in which for all $i$, $(r_{i}(x)-r_{h}(x))$ is not zero. Then, bringing $(r_{w}(n)-r_{h}(n))h_{w}(n)$ in (4) to the other side yields

[TABLE]

The sequence $(r_{w}(n)-r_{h}(n))h_{w}(n)$ is hypergeometric, and by the induction hypothesis it follows that there are $i,j$ and a rational function $\tilde{r}$ with $(r_{i}(n)-r_{h}(n))h_{i}(n)=\tilde{r}(n)(r_{j}(n)-r_{h}(n))h_{j}(n)$. Dividing by $r_{i}(n)-r_{h}(n)$ proves the claim. For the case that there is an $i$ with $(r_{i}(x)-r_{h}(x))=0$, the left hand side of (4) is a sum of fewer than $w$ hypergeometric terms and the right hand side is hypergeometric. The induction hypothesis then again yields suitable $i,j$ and $r(x)$. ∎

Example 4.4.

The sums $2n!+(n+3)!$ and $n!+(n+\frac{1}{2})^{\underline{n}}-n!$ are hypergeometric, whereas $1+n!$ is not.

The next lemma gives a characterization of when the quotient of two hypergeometric sequences is a rational function sequence. Together with Lemma 4.3, this then will yield the algebraic independence of certain hypergeometric sequences in Lemma 4.6.

Lemma 4.5.

Let $\zeta_{1},\dots,\zeta_{\ell}\in\mathbbm{K}$ be such that for all $i,j=1,\dots,\ell$ with $i\neq j$, we have $\zeta_{i}-\zeta_{j}\notin\mathbbm{Z}$. Then for $k_{1},\dots,k_{\ell}\in\mathbbm{N}$, $c_{1},\dots,c_{\ell}\in\mathbbm{N}$, and $\theta_{1},\theta_{2}\in\mathbbm{K}$, there is a rational function $r(x)\in\mathbbm{K}(x)$ such that

[TABLE]

if and only if $\theta_{1}=\theta_{2}$ and $(k_{1},\dots,k_{\ell})=(c_{1},\dots,c_{\ell})$.

Proof.

If $\theta_{1}=\theta_{2}$ and $(k_{1},\dots,k_{\ell})=(c_{1},\dots,c_{\ell})$, then we can set $r(x)=1$. For the other direction, we have

[TABLE]

A hypergeometric term $h$ is a rational function if and only if its shift quotient $h(x+1)/h(x)$ can be written in the form

[TABLE]

with $f,g\in\mathbbm{K}[x]$. Therefore, for any root in the numerator of $q(x)$ there is a root in integer distance in the denominator of $q(x)$, which, by the condition on the $\zeta_{i}$, is not possible if $\theta_{1}\neq\theta_{2}$ or $(k_{1},\dots,k_{\ell})\neq(c_{1},\dots,c_{\ell})$ ∎

Lemma 4.6.

Let $\theta_{1},\dots,\theta_{s}\in\mathbbm{K}$ and $\zeta_{1},\dots,\zeta_{\ell}\in\mathbbm{K}$. The sequences $(n+\zeta_{1})^{\underline{n}},(n+\zeta_{2})^{\underline{n}},\dots,(n+\zeta_{\ell})^{\underline{n}}$ are algebraically independent over $\mathbbm{K}(n,\theta_{1}^{n},\dots,\theta_{s}^{n})$ if and only if there are no $i,j\in\{1,\dots,\ell\}$, $i\neq j$ such that $\zeta_{i}-\zeta_{j}\in\mathbbm{Z}$.

Proof.

If there are $i,j\in\{1,\dots,\ell\}$, $i\neq j$ with $\zeta_{i}-\zeta_{j}=k\in\mathbbm{Z}$, then we get the algebraic relation

[TABLE]

Conversely, let $p$ be a nonzero polynomial over $\mathbbm{K}(n,\theta_{1}^{n},\dots,\theta_{s}^{n})$ in $\ell$ variables. We can write $\operatorname{denominator}(p)\cdot p((n+\zeta_{1})^{\underline{n}},\dots,(n+\zeta_{\ell})^{\underline{n}})$ as a sum of the form

[TABLE]

Assume that $p((n+\zeta_{1})^{\underline{n}},\dots,(n+\zeta_{\ell})^{\underline{n}})=0.$ Then, by Lemma 4.3, there have to be terms $(i,k),(j,c)\in\mathbbm{N}\times\mathbbm{Z}^{\ell}$, $(i,k)\neq(j,c)$ and a rational function $r(x)\in\mathbbm{K}(x)$ with

[TABLE]

By Lemma 4.5, this can only be the case if there are $\zeta_{i}\neq\zeta_{j}$ in integer distance, which contradicts the condition on the $\zeta_{i}$. ∎

Example 4.7.

Let $h_{1},h_{2},h_{3}$ be hypergeometric sequences given by $h_{1}(0)=h_{2}(0)=h_{3}(0)=1$ and

[TABLE]

The closed forms then are

[TABLE]

From Lemma 4.6 it follows that $h_{1},h_{2}$ are algebraically independent over $\mathbbm{K}$, but $h_{1},h_{3}$ are not.

Lemma 4.6 allows us to represent the sequences arising in extended P-solvable loops as rational function sequences over the field $\mathbbm{K}(n,\theta_{1}^{n},\dots,\theta_{s}^{n})$ as follows: Let $v_{1},\dots,v_{m}$ be of the form (3) and let $\tilde{Z}=\{\tilde{\zeta}_{1},\dots\tilde{\zeta}_{k}\}$ be a subset of $Z=\{\zeta_{1},\dots,\zeta_{\ell}\}$ such that there are no $i,j=1,\dots,k$, $i\neq j$, with $\tilde{\zeta}_{i}-\tilde{\zeta}_{j}\in\mathbbm{Z}$ and for each $\zeta\in Z\setminus\tilde{Z}$ there exists an $i$ such that $\tilde{\zeta}_{i}-\zeta\in\mathbbm{Z}$. Let $z_{1},\dots,z_{\ell}\in\mathbbm{K}[x,y_{1},\dots,y_{k}]$ be such that

[TABLE]

for all $n\in\mathbbm{N}$ and $i=1,\dots,\ell$. Then there exist $k_{1},\dots,k_{m}\in\mathbbm{Z}^{\ell}$ with

[TABLE]

Substituting variables $v_{i}$ for $v_{i}(n)$, $h_{i}$ for $(n-\tilde{\zeta}_{i})^{\underline{n}}$, $e_{i}$ for $\theta_{i}^{n}$ and $x$ for $n$ then gives

[TABLE]

where $r_{i}$ is a rational function over $\mathbbm{K}$ in $1+s+k$ variables. We now can compute the ideal of all algebraic dependencies among the program variables of a P-solvable loop as the ideal of algebraic relations among rational functions.

Proposition 4.8.

Let $(v_{1}(n))_{n\in\mathbbm{N}},\dots,(v_{m}(n))_{n\in\mathbbm{N}}$ be sequences of the form (3) and consider the corresponding rational functions $r_{1},\dots,r_{m}$ in $\mathbbm{K}(x,e_{1},\dots,e_{s},h_{1},\dots,h_{k})$ as above. For each $i=1,\dots,m$, write $r_{i}=f_{i}/g_{i}$ with coprime polynomials $f_{i},g_{i}$ over $\mathbbm{K}$. Denote by $I(\theta_{1}^{n},\dots,\theta_{s}^{n})$ the ideal of algebraic relations among $\theta_{1}^{n},\dots,\theta_{s}^{n}$ in $\mathbbm{K}[e_{1},\dots,e_{s}]$. Then the ideal of algebraic relations among the sequences $(v_{1}(n))_{n\in\mathbbm{N}},\dots,(v_{m}(n))_{n\in\mathbbm{N}}$ in $\mathbbm{K}[v_{1}\dots,v_{m}]$ is given by

[TABLE]

Proof.

The proposition follows immediately from the fact that the ideal of algebraic dependencies among a set of rational functions

[TABLE]

in the polynomial ring $\mathbbm{K}[y_{1},\dots,y_{m}]$ is given by

[TABLE]

and that by Lemma 4.6 there are no algebraic relations over the field $\mathbbm{K}(n,\theta_{1}^{n},\dots,\theta_{s}^{n})$ among the terms $(n-\tilde{\zeta}_{i})^{\underline{n}}$ with $\tilde{\zeta}_{i}$ as above for $i=1,\dots,k$. ∎

Example 4.9.

We compute the ideal of algebraic relations among $a,b,c,d$ given in Example 3.3. First, we compute the ideal of algebraic relations among $(-1)^{n},2^{n},3^{n}$ and $6^{n}$ with corresponding variables $e_{-1},e_{2},e_{3},e_{6}$. We get

[TABLE]

Now we can compute the ideal of algebraic relations among $a,b,c,d$ by adding the relations $a-(k_{1}e_{-1}-k_{2}e_{6})f,k_{1}+k_{2}-a(0),-k_{1}+6k_{2}-a(1),b-b(0)e_{3},c-c(0)e_{2}f,d-d(0)f$, where $f$ is used to model $n!$, and eliminate the variables $k_{1},k_{2},e_{-1},e_{2},e_{3},e_{6}$ and $f$.

[TABLE]

For instance, with the starting values $a(0)=2,a(1)=5$ and $b(0)=c(0)=d(0)=1$ we get the relation

[TABLE]

with

[TABLE]

Remark 0.

Proposition 4.8 can easily be turned into an algorithm with the help of Gröbner bases, which allow computing a set of generators for the sum of ideals and also the elimination of variables. While computationally demanding, the use of Gröbner bases is viable in part because of the highly optimized tools that are available in modern computer algebra systems and in part because, as observed empirically in our experiments, the polynomial systems arising in practice in this context are typically small and easy to compute.

5. Implementation

The techniques presented in this paper are implemented in the open source Mathematica software package Aligator111Aligator requires the Mathematica packages Hyper (Petkovšek, 1998), Dependencies (Kauers and Zimmermann, 2008) and FastZeil (Paule and Schorn, 1995), where the latter two are part of the compilation package ErgoSum (for Symbolic Computation., 2016). (Kovács, 2008), available for download at

https://ahumenberger.github.io/aligator/

We give an illustrative example of the provided facilities.

Example 5.1.

We compute the ideal of algebraic relations among the program variables $a,b,c,d,e,f$ as given in the following loop. The loop exhibits two first-order and two second-order recurrence relations ($a,e$ and $b,d$ resp.), which Aligator could not handle before. Furthermore we have two first-order C-finite recurrence relations ($c,f$).

• • In[1]:=

Aligator[**  WHILE[True,****    a := 3(n + $\displaystyle\frac{\text{3}}{\text{2}}$)a;****    s1 := s2; s2 := b;****    b := 5($\displaystyle\frac{\text{3}}{\text{2}}$ + n)s2 - $\displaystyle\frac{\text{3}}{\text{2}}$(1 + 2n)(3 + 2n)s1;****    c := -3c + 2;****    t1 := t2; t2 := d;****    d := 4(4 + n)t2 - 3(3 + n)(4 + n)t1;****    e := (n + 4)e;****    f := 2f],****  LoopCounter -> n,****  IniVal -> {****    t1 := 1; t2 := 1;****    s1 := 1; s2 := 2;****    a := 3; b := 1;****    c := 1; d := 3;****    e := 2; f := 5}]**

The input is given to Aligato\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert**’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:** in form of a while loop, and two optional arguments: LoopCounte\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert**’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:** (default: \lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert**’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:) and IniVa\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:** (default: {\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert**’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:). The former is for specifying which variable within the loop corresponds to the loop counter, whereas the latter is for specifying the initial values of the program variables. If no initial values are given, then the invariants contain the starting values in the form of a[0\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:**, representing the initial value of $a$.

The following output of Aligato\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert**’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:** is a conjunction of the elements of the Gröbner basis of the ideal of all algebraic relations among $a,b,c,d,e$ and $f$. The loop counter is eliminated via Gröbner basis computation.

• • Out[1]=

2 d == 3 e &&(a + b) (2 d - 3 e) == 0 &&2 a d f == 3 a e f &&450 a b $\displaystyle\text{(1\ -\ 2\ c)}^{\text{2}}$ + 225 $\displaystyle\text{b}^{\text{2}}$ $\displaystyle\text{(1\ -\ 2\ c)}^{\text{2}}$ +  $\displaystyle\text{a}^{\text{2}}$ (225 $\displaystyle\text{(1\ -\ 2\ c)}^{\text{2}}$ - 16 $\displaystyle\text{f}^{\text{2}}$) == 0

Note that the second and third invariant are consequences of the first one. By setting the option GroebnerReduce -¿ Tru\lsthk@PreSet\lsthk@TextStyle__mmacells_lst_init:n\lst@FVConvert**’missing@nil\lst@ReenterModes\lst@PrintToken\lst@InterruptModes__mmacells_lst_deinit:** a reduced Gröbner basis is computed which does not contain redundant elements.

6. Conclusion and Future Work

We extended the class of P-solvable loops to include sums and products of hypergeometric and C-finite sequences. This was made possible by identifying algebraically independent factors in hypergeometric terms and then viewing the sequences in question as rational function sequences over a transcendental field extension. The implementation in Mathematica underlines the practicality of the approach.

There are several promising directions in which we plan to expand this line of research. Obviously, it is very desirable to include more types of recurrences in P-solvable loops. These include further subclasses of the class of holonomic sequences as well as partial and non-linear recurrence equations. It is advisable to conduct a careful study on which kind of recurrences are relevant in practice and also good-natured from a mathematical perspective. Uncoupling techniques for systems of recurrence equations can also prove to be helpful in this context.

Another possible extension is to consider nested loops. With the help of $\Pi\Sigma^{*}$-theory (Schneider, 2017), it might be possible to derive invariants for the outermost loop, although the inner loops are not P-solvable by themselves.