Integration of multivariate rational functions given by straight-line programs GUILLERMO MATERA∗†‡ Departamento de Matem´ aticas, Fac. de Ciencias Exactas Univ. de Buenos Aires (1428) Buenos Aires, ARGENTINA. [email protected]

Abstract We present a well-parallelizable algorithm which, taking a straight-line program for the evaluation of a vectorial field of rational functions of Q(X1 , . . . , Xn ) as input, decides whether they allow a rational potential function and, in case of affirmative answer, computes it as output. We introduce a mixed model of representation of polynomials to allow the application of integration techniques and show how to perform some basic operations with it. The algorithm is presented as a family of arithmetic networks of polynomial size and polylogarithmic depth in the degree of the occurring polynomials.

1

Introduction

This paper is concerned with the complexity problem of performing integration on polynomials (or rational functions) represented by straight-line programs. The choice of the representation of multivariate polynomials and rational functions is one of the most important tasks in computer algebra systems. The “classical” model is the dense representation: the data structure which represents a polynomial is supplied with one unit of memory space for each possible coefficient in it; in other words we a ¡ represent ¢ polynomial F ∈ k[X1 , . . . , Xn ] of degree at most d by the vector of all its d+n possible n ¡ ¢ coefficients. Since d+n ≤ dn , we observe that to the polynomial F corresponds in dense n representation a data structure of size O(dn ). The amount of space requiered by this exponential behaviour may become intractable some problems from a computational point of view. Hence, several different models have been suggested. A valid alternative is given by the straight-line programs. They constitute an algebraic “modelization” of the concept of algorithm without branchings and have proved to be a very useful tool to improve known ∗ Partially

supported by Research Project PB93-0472-C02-02 supported by GDR de Calcul Formel Medicis (CNRS) France ‡ EC Project ESPRIT BRA Contract 6846 PoSSo † Partially

1

bounds of some algorithmic algebraic problem. A straight-line program can be defined as a (finite) sequence of rational functions of Q(X1 , . . . , Xn ) namely Q1 , . . . , Qα (intermediate results), such that every Qi of the sequence is either a variable X1 , . . . , Xn or an element of Q (parameter), or either the sum, the difference, the product or the division of two preceeding polynomials in the sequence (for precisions and elementary properties of straightline programs we refer to [7, 8, 13, 14]). An example of the convenience of this representation is the symbolic n × n determinant: its dense representation has size n! but it can be computed with O(n3 ) operations by techniques of gaussian elimination. It is often convenient to use the representation of a straight-line program by a directed acyclic graph. The inputs are the sources of the graph, the node for every Qi is labelled by its operation and it has 2 input edges numbered from the predecessors of Qi . Two natural measures for a straight-line program β are its size or sequential time (total number of gates) and its depth or parallel complexity (length of the longest path in β). We shall use the nonscalar measures which counts the steps where an essential operation is performed, that is, the product between two nonconstant polynomials and the division with nonconstant denominator. A justification of the choice of this measure would be the fact that, from a complexity point of view, nonscalar time “dominates” other known measures. Besides, nonscalar depth provides the best bounds for the degree of the represented polynomials. Naturally a question arises: whether polynomials given by straight-line program can be manipulated in order to perform operations like greatest common divisor, differentation, etc. In this sense, several important transformations are known (see [11, 15]). A fundamental one by Baur-Strassen [1] allows to compute all first derivatives with growth in length only by a constant factor (without even the need of a degree bound). It has enormous consequences in complexity theory and computer science. Thus, we ask for the “inverse” problem, that is, the integration: given a straight-line program to evaluate a rational function (or more generally, a gradient field of rational functions from Q(X1 , . . . , Xn )), construct a straight-line program to evaluate its integral (their potential function). Here we can hardly expect the Baur-Strassen’s kind of bounds. There is some evidence that integration is conected with some problems involving straight-line programs which are known to be NP-hard. Hence, an upper bound of order (nL)O(1) for the process of integration would probably imply that P=NP. On the other side, it is well known that the dense representation of an n-variate polynomial of degree d can be obtained from a straight-line program of nonscalar length L by means of dn L nonscalar operations. Since we are interested in the advantages of the straight-line program structure, this is a natural upper bound to the results to be obtained. In this work we solve this problem by means of a well-parallelizable algorithm with (ndL)O(1) nonscalar size. More precisely, we shall prove the following: Theorem 1 Let F1 , . . . , Fn , G1 , ..., Gn be polynomials in Q[X1 , . . . , Xn ] of highest degree d, given by an arithmetic circuit without divisions of nonscalar size L and nonscalar depth `. Suppose further that Gi 6= 0 for i = 1, ..., n. 2

Let P1 :=

Fn F1 , ..., Pn := . G1 Gn

Then there exists a randomized algorithm with nonscalar sequential time complexity (Lnd)O(1) and nonscalar parallel time complexity O(` + log3 d) which: 1. decides whether the vector field (P1 , ..., Pn ) allows a potential function Q 2. being (1) true, whether there exists a potential function Q belonging to Q(X1 , . . . , Xn ) 3. in case of affirmative answer in (2), constructs a straight-line program that evaluates Q on the domain of P1 , . . . , Pn . The network can be constructed by a probabilistic algorithm of nonscalar size O(Lnd)O(1) and nonscalar depth O(` + log3 d).

1.1

On the complexity model

Our algorithms will be represented by randomly generated (alternatively, non uniform) families of arithmetic networks over Q (see [7]). These arithmetic networks are directed acyclic graphs whose internal nodes are labelled by arithmetic operations of Q, comparisons in Q, selections of elements of Q (associated to comparisons) and Boolean operations corresponding to propositional logic. We consider the nonscalar measures as real valued functions depending on input parameters of the given straight-line program such that d, n, L and ` (highest degree, number of variables, nonscalar size and nonscalar depth) and we try to analyse their asymptotic behaviour.

2

The mixed representation model

All the algorithmic tools developed for rational integration seem strongly dependent on the knowledge of the coefficients of the polynomials to integrate. However, the integration of a multivariate rational function with respect to one variable only requires the dense structure with respect to that variable. This is a key observation for this work, because it suggests a model well adapted to our requirements.

2.1

A mixed model: dense representation and straight-line programs

As we mentioned in the introduction, the dense representation is too expensive to compute. On the other side, the straight-line program code does not seem to allow univariate integration. The mixed representation is an intermediate model which avoids these two difficulties. It consists of the representation of n-variate polynomials with respect to one fixed variable, giving the coefficients (which are themselves n − 1-variate polynomials) by means 3

of a straight-line program. It is less explicit than dense representation, which allows its computation from a straight-line program, but enough to allow univariate integration. Let Xi be a fixed indeterminate (the variable of “integration”), and let F be a polynomial d X Fk ·Xik , where F0 , . . . , Fd ∈ Q[X1 , .., Xi−1 , Xi+1 , .., Xn ]. in Q[X1 , . . . , Xn ]. Let us write F = k=0

Then the mixed representation of F with respect to Xi is defined as a straight-line program β in Q[X1 , .., Xi−1 , Xi+1 , .., Xn ] which computes F0 , . . . , Fd . The next lemma computes the nonscalar cost required to transform a straight-line program into a mixed representation. Lemma 1 Let X1 , . . . , Xn , T be distinct indeterminates over Q and let F ∈ Q[X1 , ..., Xn , T ] a polynomial of degree d given by a division-free straight-line program β in Q[X1 , ..., Xn , T ] of nonscalar size L and nonscalar depth `. Then there exists a division-free straight-line program in Q [X1 , ..., Xn ] of nonscalar size O(dL) and nonscalar depth ` computing the mixed representation of F with respect to T . Proof Since deg(F ) = d, in particular degT (F ) ≤ d, and therefore we can find the writing in the variable T by interpolating in d+1 points of Q. Doing all the interpolations in parallel, with (nonscalar) sequential cost (d + 1)L and parallel cost `, we obtain the a straight-line program in Q[X1 , ..., Xn ] which computes the coefficients of F with respect to T . 2

2.2

Correct test sequences

A crucial point in algorithms treating with polynomial expressions is the “equality checking” between two polynomials which appear as intermediate results. Suppose for the moment that one wishes to know if a n-variate polynomial F with degree d produced by a suitable straight-line program is the zero polynomial. A straightforward procedure for this purpose consists in the succesive evaluation of the straight-line program in the integer points of a box of cardinality dn (interpolation) and is therefore exponential in the degree of F . This main difficulty can be partially avoided by means of a fine result due to J.Heintz & C.P. Schnorr [9] (see also [10]). Let d, L, ` be non negative integers. We consider the set W (d, n, L, `) of all polynomials in Q[X1 , . . . , Xn ] of degree bounded by d which can be evaluated by an straight-line program in Q(X1 , . . . , Xn ) of size and depth bounded by L and `. Let m ∈ IN, we say that a sequence γ := (γ1 , . . . , γm ) (where γi ∈ Qn for 1 ≤ i ≤ m) is a correct test sequence for W (d, n, L, `) if the polynomial 0 is the unique polynomial of the set W (d, n, L, `) that annihilates γ1 , γ2 , . . . and γm . With these notations we have: Theorem 2 ([9, Theorem 4.4]) Let Γ be a subset of Q of cardinality (2`+1 − 2)(2` + 1)2 , m := 6(L + n)(L − (n + 1)) and denote by σ(d, n, L, `, Γ) the subset of Qnm of all correct test sequences of W (d, n, L) contained in Γnm . Then the following inequality holds: m (#Γ)nm (1 − (#Γ)− 6 ) ≤ #σ(d, n, L, `, Γ).

4

Suppose now that F ∈ Q[X1 , . . . , Xn ] is an arbitrary polynomial of degree bounded by d and which can be computed by a straight-line program in Q[X1 , . . . , Xn ] without divisions of size L and depth `. From the previous proposition one deduces: Corollary 1 [6, Corollaire 2.1]) There exists an arithmetical network over Q working in size O(Lm) = O(L(L+n)2 ) and depth O(`) which decides if F is the polynomial zero. Moreover, this network can be constructed by a probabilistic algorithm in sequential time O(L(L + n)2 ) 1 and parallel time O(`); its failure probability is smaller than ε = 262144 . Let us remark that the proposition guarantees the existence of at least one correct test sequence of size m which depends polynomially on the parameters (and not exponentially as in interpolation). The construction of correct test sequences depends only on the parameters; therefore if one knows a priori upper bounds for all polynomials in the algorithm, it is quite natural to consider the construction of the correct test sequences as a preprocessing computed once for ever. The corollary introduces a probabilistic algorithm (choosing randomly a correct test sequence) which unfortunately has a failure probability (ε < 12 ): if one is interested in the application of this procedure in order to check if a given polynomial is zero, the answer NO is always true but YES involves a failure probability smaller than ε.

2.3

The manipulation of polynomials given by mixed representations

This section is devoted to perform some manipulation with the mixed model, keeping polynomial-size polylogarithmic-depth complexity bounds. We shall pay attention to the computation of greatest common divisor (gcd) and square-free factorisation because the integration techniques to be applied are based on them. Let us fix some notations. Throughout the remaining of the section, let d, e, n be fixed d X natural numbers. We shall deal with polynomials F (X1 , ..., Xn , T ) = Fk (X1 , ..., Xn )T k and G(X1 , ..., Xn , T ) =

e X

k=0

Gk (X1 , ..., Xn )T k of Q[X1 , ..., Xn ][T ] with degree d,e in T re-

k=0

spectively, given by their mixed representations. Their coefficients {Fk , Gk } will be given by a straight-line program β of nonscalar size L and nonscalar depth `. We shall assume furthermore that we are given an n-tuple α ∈ Qn such that Fd (α) 6= 0 and Ge (α) 6= 0. The idea of the next algorithms is the following: we start running β on input α ∈ Qn to evaluate the coefficients of the input polynomials, and then we operate with those coefficients in order to compute the required output. The philosophy is to carry our problems to a linear algebra setting where they can be solved by means of well-parallelizable networks of polynomial size.

5

2.3.1

Computation of greatest common divisor

We shall follow here the ideas of [3]: given polynomials f ,g ¡in Q[T ], with deg(f ) := d ≥ ¢ deg(g) := e, if f is not a constant multiple of g, then deg gcd(f, g) := j = min{k ∈ IN; ∃r, s in Q[T ], deg(r) < e − k, deg(s) < d − k and deg(rf + sg) = k}. This can be translated into a linear system, whose matrix, namely Pj (f, g), is known as the j-th principal submatrix of the Sylvester matrix of f and g. It can be obtained from this latter by deleting its last k columns of coefficients of f , its last k columns of coefficients of g and its last¡2k rows. ¢ Since deg gcd(f, g) = ¡min{k; det(P k ) 6= 0} ([7, theorem 2.3]), we can multiply both ¢ sides of the system by Adj Pj (f, g) and conclude that ¡ ¢ ¡ ¢ det Pj (f, g) (re−j−1 , · · · , r0 , sd−j−1 , · · · , s0 )t = Adj Pj (f, g) (0, · · · , 0, 1)t ¡ ¢ (let us remark that the right hand¡ side of the last equation is the last column of Adj P (f, g) ). j ¢ This allows us to compute a det Pj (f, g) -multiple of the coefficients of r(T ) and s(T ) in a suitable construction. ¡ ¢ In our setting, since Fd (α)·Ge (α) 6= 0 it follows that Pk F (α, T ), G(α, T ) = Pk (F, G)(α) (where Pk (F, G) is the k-th principal submatrix of F and G as polynomials in Q(X1 , . . . , Xn )[T ]). We design the following algorithm: - Input: A straight-line program β and an n-tuple α ∈ Qn . ¡ ¢ ¡ ¢ - Run β on α to obtain the coefficients Fd (α), . . . , F0 (α) and Ge (α), . . . , G0 (α) . µ ³ e−i−1 d−i−1 ´¶ X (i) X (i) - Compute det Pi F (α, T ), G(α, T ) , A(i) (T ) = rk T k and B (i) (T ) = sk T k (i)

k=0 (i)

(i)

(i)

k=0

in parallel for every i, where the vector (re−i−1 , · · · , r0 , sd−i−1 , · · · , s0 )t (α) is the last µ ³ ´¶ column of Adj Pi F (α, T ), G(α, T ) . - Compute (in parallel) the coefficients of A(i) (T ) · F (α, T ) + B (i) (T ) · G(α, T ) n¯ ³ ´ o ¯ - Set j := min i¯Pi F (α, T ), G(α, T ) 6= 0 . - Compute µ ³ ´¶ ³ ´ det Pj F (α, T ), G(α, T ) ·gcd F (α, T ), G(α, T ) = A(j) (T )·F (α, T )+B (j) (T )·G(α, T ) 2 During its execution, the algorithm computes in parallel – the prin¡ O(d ) determinants ¢ cipal subresultants and the last column of Adj Pi (F, G) (α) – of matrices whose maximal size is d × d and performs O(d) comparisons (only O(log d) sequentially) and O(d) products (in parallel). The determinants are computed by means of effective methods of linear algebra due to Berkowitz [2] with sequential time d3 and parallel time log2 d. In conclusion, we have:

6

Proposition 1 There exists a division-free ¢ arithmetic network which, on input α, computes ¡ a scalar multiple of gcd F (α, T ), G(α, T ) by means of O(L + dO(1) ) sequential and O(` + log2 d) parallel nonscalar steps. 2.3.2

Square-free Factorisation

Let f be a polynomial in Q[T ], we can write f (T ) =

d Y ¡

fk (T )

¢k

where, for every k, the

k=1

fk ’s belong to Q[T ], are square-free and satisfy: gcd(fi , fj ) = 1 for every 1 ≤ i < j ≤ d. This is known as the square-free factorisation of f . In order to get low parallel complexity, we take the ideas of [7] and compute the factors fk in the following way: for 0 ≤ k ≤ d, let us define gk := gcd(f, f 0 , . . . , f (k) ) (hence we have 2 gk = fk+1 fk+2 . . . fdd−k ); therefore the factors fk of the square-free factorisation satisfy the equality: gk−1 gk µ ¶ fk = gk−1 0 gcd , gk−1 gk Following the above scheme, we design an algorithm to compute square-free factorisation on polynomials in mixed representation: - Input: A straight-line program β and an n-tuple α ∈ Qn ¡ ¢ - Run β on α to obtain Fd (α), . . . , F0 (α) . d

X ∂kF T ) = i · · · (i − k + 1)Fi (α, T )T i−k and (α, ∂T k i=k ³ ´ ∂kF all Gk (α, T ) := gcd F (α, T ), . . . , (α, T ) by a tree of gcd’s: start computing (in k ³ ∂ 2j−1 F ∂T ´ ∂ 2j F parallel) all possible gcd (α, T ), (α, T ) and, every next step, compute ∂T 2j−1 ∂T 2j gcd of two (successive) previous results. (Let us observe that only O(log d) parallel steps of computation of greatest common divisors are needed). ³ G (α, T ) ∂G ´ k k−1 - Compute (in parallel) all gcd , (α, T ) . Gk−1 (α, T ) ∂T - For every k, compute (in parallel)

Gk (α, T ) Gk−1 (α, T ) - Compute (in parallel) all Fk (α, T ) = ³ G (α, T ) ∂G ´. k−1 k , (α, T ) gcd Gk−1 (α, T ) ∂T Let us remark that the last step of the algorithm consists of a division of two polynomials whose result is known to be a polynomial a priori. This task can be rewritten as a linear system of suitable size, whose parameters depend on the coefficients of the polynomials 7

involved. Hence, the size and depth of the whole procedure are increased in a way according to our requirements. Since derivation in T is not counted in the nonscalar model, applying the result of proposition 1 we have: Proposition 2 There exists a division-free arithmetic network of nonscalar size O(L + dO(1) ) and nonscalar depth O(` + log3 d) which, on input α, performs the square-free factorisation F (α, T ).

3 3.1

On the integration of rational functions Algorithmic tools of integration

We start this section with a short account of the techniques we shall apply to integrate rational functions, which are taken from [5]. They are a strongly simplified fragment of the Liouville’s theory of formal integration. A differential field is a field F , together with a derivation on F , that is, a map of F into itself, usually denoted by a → a0 , such that (a + b)0 = a0 + b0 and (ab)0 = a0 b + ab0 . If a, b are elements of a differential field F , we say that b is aRprimitive of a, or b is an integral of a, if b0 = a. Sometimes we denote a primitive of a by a. Likewise, if a is nonzero, b is a a0 logarithm of a if b0 = . By a differential extension field of a differential field F we mean a a differential field which is an extension field of F whose derivation extends the derivation on F . Let us remark that if F is a differential field of characteristic zero and K is an algebraic extension field of F , then the derivation on F can be extended to a derivation on K in a unique way. We shall assume all algebraic extension fields provided with this uniquely determined differential structure. Let E be an algebraic extension of Q and X an indeterminate over Q. In the sequel we shall deal with Q(T ) provided with its natural differential structure, given by the unique derivation defined by the rule T 0 = 1. We say that a primitive of a rational function f is expressible in E (by means of adjunctions of logarithms) if there exists a primitive of f on a differential extension field of E(T ), which is obtained by adjoining some logarithms of rational functions of E(T ). A first observation is that an algebraic extension is sufficient to express the integral of a rational function. In fact, the decomposition field of its denominator, namely E, is enough. Furthermore, we can write the primitive as a sum of a rational function g in E(T ) plus a sum of logarithms of rational functions from E(T ). This construction is effective and provides a fundamental tool to this work: the rational function g can be choosen in Q(T ). We summarize all this in the following theorem. p(T ) where p, q ∈ Q[T ] are q(T ) coprimes. Let E be the Rdecomposition field Rof q(T ). Then, there exists a primitive of f which can be written as f (T )dT = g(T ) + h(T )dT where g, h ∈ Q(T ) and the primitive of h is a sum of logarithms of elements of E(T ) (up to constants).

Theorem 3 Let f be a rational function de Q(T ), let f (T ) =

8

Proof ([5, Chapter 1,Theorems 1 and 2]) d Y ¡ ¢j qj (T ) be the square-free factorisation of q. Let q(T ) = j=1

The main tool in the proof is the partial fraction decomposition of f with respect to the square-free factorisation of q: this consists of the writing f (T ) = a0 (T ) +

d X aj (T ) ¡ ¢j j=1 qj (T )

where the aj ’s are polynomials of Q[T ] verifying deg(aj ) < deg(qj ). By integration in the expression above we have Z

Z f (T )dT =

a0 (T )dT +

d Z X j=1

aj (T ) ¡ ¢j dT qj (T )

Since a0 (T ) is a polynomial from Q[T ], we know that it has a primitive in Q[T ]. aj (T ) r(T ) On the other hand, all ¡ where r, s are ¢j are rational functions of type j s (T ) qj (T ) polynomials of Q[T ] coprime and s is square-free. We write r(T ) = b1 (T )·s(T )+b2 (T )·s0 (T ) where deg(b1 ) < deg(s0 ) and deg(b2 ) < deg(s). Then we have µ ¶0 r b 2 · s0 b1 1 b2 b = + = − + j−1 sj sj sj−1 j − 1 sj−1 s where b = b1 +

1 0 b . Thus j−1 2 Z

r 1 b2 =− · j−1 + j s j−1 s

Z

b sj−1

Iterating this procedure we obtain an expression of the from: Z

r(T ) dX = g˜(T ) + sj (T )

Z ˜ b(T ) dT s(T )

where g˜ ∈ Q(T ) and ˜b ∈ Q[T ]. Since s is square-free, its integral is a sum of logarithms (of rational functions of E(T)) and hence we have proven the theorem. 2

3.2

Decomposition in partial fractions

The proof of theorem 3 contains all essential methodical points we need in our algorithmic setting. As it can be seen, we need the decomposition of rational functions given by mixed representations into partial fractions with respect to square-free factorisation. We start with the following: 9

Lemma 2 Let p, q be polynomials in Q[T ] such that deg(p) < deg(q) and let q = q1 · · · qs be the square-free factorisation of q. Then, there exist unique polynomials r1 , . . . , rs with deg(ri ) < deg(qi ) for i = 1, . . . , s verifying: s

p X ri = q q i=1 i Proof [12, Theorem 8, Chapter 5, Paragraph 5].

2

We can replace the condition above by the following: p=

s X i=1

ri

q qi

This latter can be rewritten as a linear system with the coefficients of r1 , . . . , rn as unknowns, whose matrix has size d×d. The uniqueness forces this matrix to be nonsingular, and applying the same ideas as in section 2.3.1 we compute the coefficients of r1 , . . . , rs by means of dO(1) arithmetic sequential operations in O(log2 d) parallel steps. Now, let F, G ∈ Q[X1 , . . . , Xn , T ] polynomials of degree d, e in T , with d < e, given by a mixed representation of nonscalar sequential time L and nonscalar parallel time `. Let α ∈ Qn such that Fd (α) 6= 0 and Ge (α) 6= 0. Set p(T ) := F (α, T ) and q(T ) := G(α, T ). First, we can compute the square-free factorisation of q(T ) and then, we can apply the scheme explained above in order to compute decomposition in partial fractions of pq with respect to the square-free factorisation of q. The results obtained are summarized in the following: Theorem 4 There exists an arithmetic network working in nonscalar sequential time O(L+ dO (1)) and nonscalar parallel time O(`+log3 d) which on input α computes the decomposition F (α, T ) in partial fractions of . G(α, T )

3.3

The “rational” decomposition

After performing decomposition in partial fractions we come into a situation of integration r of a rational function of type j , where r, s ∈ Q[T ], s is square-free and deg(r) < deg(sj ). s We reduce this problem to the computation of a polynomial b ∈ Q[T ]Zand a rational Z function r b g ∈ Q(T ), such that deg(b) < deg(s) and an equality of the type =g+ holds. j s s r This expression, which we call the rational decomposition of sj , completes theRprocess of rational integration. g and h are called the rational and the logarithmic part of srj . We shall perform this rational decomposition in two stages: first, we compute a kind of n−1 X “Taylor” representation of r as r = ri ·si where deg(ri ) < deg(s) for every i = 1, . . . , n−1. i=1

The existence and uniqueness of this representation is straightforward, and an efficient parallel computation can be performed if one takes into account the following observations: 10

- the remainder of the division of r by si is ri−1 si−1 + · · · + r0 . - the quotient of ri−1 si−1 + · · · + r0 by si−1 is ri−1 . - these two operations can be done in parallel for every i = 1, . . . , j − 1. Further, we obtain the rational decompositions under the additional assumption that deg(r) < deg(s) which is guaranteed by the previous stage. We first prove the existence and uniqueness of this representation: Lemma 3 Let r, s polynomials of Q[T ] such that s Then, there exists unique polynomials r0 , r1 , . . . , rj−1 i = 0, . . . , j − 1 and Z Z j−1 X r ri = + sj si i=1

is square-free and deg(r) < deg(s). such that deg(ri ) < deg(s) for every r0 s

Proof Let us prove the uniqueness. Given two rational decompositions, we can write their difference in the form: Z j−1 X ri r0 0= + i s s i=1 where deg(ri ) < deg(s) for every i = 0, . . . , j − 1. Taking derivative and multiplying the last expression by sj we have:

0=

j−2 X £ 0 j−i ¤ £ 0 ¤ ri · s − i · ri · s0 · sj−i−1 − rj−1 · s − (j − 1) · rj−1 · s0 + r0 · sj−1 i=1

This identity implies that s divides rj−1 and since deg(rj−1 ) < deg(s), rj−1 = 0. By inductive arguments, it follows that rj−2 = · · · = r1 = 0, and finally that r0 = 0. The existence is contained in the proof of theorem 3. 2 Now we consider the problem of parallelizing this procedure. The proof of the Theorem 3 shows that Z

r 1 b =− · + sj j − 1 sj−1

Z a+

b0 j−1 sj−1

where a, b satisfies the relation r = a · s + b · s0 and deg(a) < deg(s0 ) and deg(b) < deg(s). We obtain the rational decomposition repeating this process. Let r(T ) := rd−1 · T d−1 + · · · + r0 and let us denote by P0 (s, s0 ) the Sylvester matrix of s and s0 . In section 2.3.1 is shown that the coefficients ad−2 , . . . , a0 , bd−1 , . . . , b0 of a, b satisfy

11

the following equation:

a





 0  ..   ...   .       0   a0      = P0 (s, s0 )−1 ·  r   bd−1   d−1   .   .   .   ..  . r0 b0 d−2

The next step consists of the computation of polynomials a(1) , b(1) of suitable degree such that b0 a(1) · s + b(1) · s0 = a + j−1 This can be done in the following way: let A be a fixed (2d − 1) × (2d − 1) matrix such that  a   0 d−2 ..   ..   .   .       0    a   d−1 A· 0 = · b a +   d−1 d−1 j−1  bd−1     .    .  .   ..  . 1 a0 + j−1 · b1 b0 (1)

(1)

(1)

(1)

Then the coefficients ad−2 , . . . , a0 , bd−1 , . . . , b0

of a(1) , b(1) can be computed as:



   (1)   ad−2 0 0  .  .. .    .   ..  .    .       (1)   0  0    a0    0 −1 0 −1 0 −1    (1)  = P0 (s, s ) ·  bd−1 + d−1 · bd−1  = P0 (s, s ) · A · P0 (s, s ) ·  r  j−1  bd−1   d−1          ..  .   ...    .  ..  1 (1) r0 a0 + j−1 · b1 b 0

The same ideas lead us to compute every a(i) , b(i) as a product of matrices. Since the iterated product of matrices works well in parallel, we compute the rational decomposition by means of dO(1) products in O(` + log3 d) parallel steps. These ideas, in a mixed representation setting, yield the following result: Theorem 5 With notations as in theorem 4, on input α the rational decomposition of F (α, T ) can be performed by a division-free arithmetic network of nonscalar size LdO(1) G(α, T ) and nonscalar depth O(` + log3 d).

12

3.4

The existence of a rational primitive

An important point in this work is to decide whether a rational function allows a rational primitive. In order to answer this question, we shall deal with the residues of the function to be integrated. From Theorem 3, we know that all the logarithms which appear in a primitive of a given rational function f ∈ Q(T ) can be expressed by means of the integral of a rational p1 (T ) function h(T ) = where gcd(p1 , q1 ) = 1, deg(p1 ) < deg(q1 ) and q1 is square-free. Let q1 (T ) d d Y X Resci (f ) q1 (T ) = (T − ci ) the factorisation of q1 over C, since h(T ) = , we have: (T − ci ) i=1 i=1 Z d Z d X X Resci (f ) h(T )dT = dT = Resci (f ) log(T − ci ) (T − ci ) i=1 i=1 1 1 ,..., are linearly inde(T − c1 ) (T − cd ) pendent over Q(c1 , . . . , cd ). Hence, the same is true for log(T − c1 ), . . . , log(T − cd ). Thus, R h(T )dT = 0 iff Resc1 (f ) = · · · = Rescd (f ) = 0. That is, we have proved: Since ci 6= cj if i 6= j, one easily deduces that

Lemma 4 f (T ) has a rational primitive in Q(T ) iff Resci (f ) = 0 for i = 1, . . . , d We can check this condition by means of the following result, due to Bronstein [4]: there exists a polynomial B(f ) whose roots are exactly the residues of f at all its poles. Furthermore, it can be computed in a completely rational way, as follows: p write f = where p, q ∈ Q[T ], gcd(p, q) = 1 and q is monic. Let q = q1 · q22 · · · qdd be q the square-free factorisation of q and let m be a natural number smaller than d such that deg(qm ) ≥ 1. Let U be a differential indeterminate over Q(T ). Compute h=

(f qm U −m )(m−1) ∈ Q(T )hU i (m − 1)!

and write h as h=

P (T, U, U 0 , . . . , U (m−1) ) Q(T, U, U 0 , . . . , U (m−1) )

where P, Q are (d + 1)-variate polynomials with coefficients in Q. Then compute µ (m) ¶ q 00 qm 0 P˜ = P T, qm , m,..., ∈ Q[T ] 2 m µ (m) ¶ 00 qm qm 0 ˜ Q = Q T, qm , ,..., ∈ Q[T ] 2 m

13

˜ qm ) ∈ Q[T ]. Hence the roots of Bm (P ) are the residues of Let Bm (P ) = ResT (P˜ − T Q, f at all its poles of order m in Q. Thus, the residues of f at all its poles are exactly the d Y roots of B(f ) := Bm (f ) over Q. m=1 R Hence we replace the condition of existence of a rational primitive f (T )dT by nonexistence of nonzero roots of B(f ). Altogether, we have: Proposition 3 With notations as in theorem 4, on input α it can be decided whether F (α, T ) allows a rational primitive by a division-free arithmetic network of nonscalar size G(α, T ) O(L + dO(1) ) in O(` + log3 d) nonscalar parallel steps.

4

The main results

Throughout this section, F1 , . . . , Fn , G1 , . . . , Gn will denote polynomials in Q[X1 , . . . , Xn ], represented by a division-free straight-line program with nonscalar size L and nonscalar F1 Fn depth `. Assume that the Gi 6= 0 for i = 1, . . . , n and set P1 := , . . . , Pn := . G1 Gn We start studying the first question of this work: whether the vectorial field (P1 , . . . , Pn ) is a gradient field. Theorem ¡ ¢ 6 There exists a randomized division-free arithmetic network of nonscalar size O (Ln)3 and nonscalar depth is O(`) which decides whether (P1 , . . . , Pn ) is a gradient field. This network can be produced by a probabilistic algorithm of sequential time O(L(L + n)2 ) and parallel time O(`). Proof: The statement of the Theorem is equivalent to prove that ∂Pj ∂Pi = ∂Xj ∂Xi for every 1 ≤ i, j ≤ n.

∂Pi ∂Pj = (Gi Gj )2 · . This can be written ∂Xj ∂Xi as a polynomial equality involving only Fi , Fj , Gi , Gj and their first order derivatives. From [1, Theorem 1], the first derivatives of all Fi , Gi ’s can be computed (in parallel) with nonscalar sequential time O(nL) and nonscalar parallel time O(`). Thus the problem is reduced to checking equalities between polynomials given by straight-line programs. Hence corollary 1 asserts there exists an arithmetical network in the conditions of the theorem which performs this task. 2 We can rephrase this condition by (Gi Gj )2 ·

Next problem is to decide whether (P1 , . . . , Pn ) – which is now supposed to be a gradient field – allows a rational potential function Q, i.e., whether there exist a potential function Q ∈ Q(X1 , . . . , Xn ) for (P1 , . . . , Pn ). 14

By means of an affine change of coordinates we can assume that the n-dimensional origin (0, . . . , 0) is in the domain of (P1 , . . . , Pn ) and that Gi (X1 , . . . , Xn ) is monic in Xi (the variable of integration) for i = 1, . . . , n. The coefficients of the matrix of this change of coordinates can be chosen from a correct test sequence of suitable size. Since Gi (0, . . . , 0) 6= 0 for every i = 1, . . . , n, there exists a positive real number r and a neighborhood D := {α ∈ Qn / k α k< r} ⊆ Qn of (0, . . . , 0) such that Gi (α) 6= 0 holds over D. This means that Q is defined over D. As it is well known, Q can be computed over D by the formula: Q(X1 , . . . , Xn ) =

n Z X i=1

Xi

Pi (X1 , . . . , Xi−1 , T, 0, . . . , 0)dT

(1)

0

Remark 1 The property of having (P1 , . . . , Pn ) a rational potential function Q(X1 , . . . , Xn ) does not depend on the point we integrate around of, but on the formal integral of the Pi ’s. In other words, if a rational function Q is a potential of (P1 , . . . , Pn ) in a neighborhood of (0, . . . , 0), it follows that Q will be a potential over the whole domain of (P1 , . . . , Pn ). The following lemma simplifies the condition of rationality of the potential function: Lemma 5 Let P1 , · · · , Pn ∈ Q(X1 , . . . , Xn ) such that (P1 , · · · , Pn ) is a gradient field. Then, there exists a rational potential function Q ∈ Z Q(X1 , . . . , Xn ) for (P1 , · · · , Pn ) iff for every i = 1, . . . , n there exists a rational primitive

Pi (X1 , . . . , Xi , 0, . . . , 0)dXi .

∂Q Proof Suppose first that there exists a rational potential Q. Since = Pi holds for ∂Xi ¡ ¢ ∂ i = 1, . . . , n, then Q(X1 , . . . , Xi , 0, . . . , 0) = Pi (X1 , . . . , Xi , 0, . . . , 0). But this means ∂Xi that Q(X1 , . . . , Xi , 0, . . . , 0) is the rational primitive required. The converse is clearly true by equation (1) and remark 1. 2 The next theorem solves the second question of this work: the estimation of upper bounds of the complexity of deciding whether our vectorial field allows a rational potential. Theorem 7 Suppose (P1 , . . . , Pn ) to be a gradient field. Then, there exists a randomized division-free arithmetic network of nonscalar size (Lnd)O(1) and nonscalar depth O(` + log3 d) which decides whether (P1 , . . . , Pn ) allows a potential function Q ∈ Q(X1 , . . . , Xn ). The network can be constructed by a probabilistic algorithm of nonscalar size (Lnd)O(1) and nonscalar depth O(` + log3 d). Proof Let r be a positive real number such that (P1 , . . . , Pn ) is defined over the open set D := {α ∈ Qn / k α k< r}. Then Q can be computed over D as: Q(X1 , . . . , Xn ) =

n Z X i=1

Xi

Pi (X1 , . . . , Xi−1 , T, 0, . . . , 0)dT

0

15

Z From remark 1 and lemma 5, it is enough to check whether there exist rational primitives Pi (X1 , . . . , Xi−1 , T, 0, . . . , 0)dT for every i = 1, . . . , n. By lemma 4 we know that the rationality of an integral of a rational function depends on whether the residues at all its poles are zero. This latter condition can be analysed applying the method explained in section 3.4. Let us fix an index i, 1 ≤ i ≤ n and an n-tuple α := (α1 , . . . , αn ) ∈ D. We are going to study the nullity of the residues of Pi (α1 , . . . , αi−1 , T, 0, . . . , 0). In order to simplify notations we set: α(i) := (α1 , . . . , αi−1 ) (X, T ) := (X1 , . . . , Xi−1 , T ) F˜i (X, T ) := Fi (X1 , · · · , Xi−1 , T, 0, . . . , 0) ˜ i (X, T ) := Gi (X1 , · · · , Xi−1 , T, 0, . . . , 0) G P˜i (X, T ) := Pi (X1 , · · · , Xi−1 , T, 0, . . . , 0) ˜ i are represented by a straight-line program in Q[X, T ] of nonscalar size L Since F˜i , G and nonscalar depth `, by lemma 1 we can compute their mixed representation with respect to T with nonscalar sequential time O(dL) and nonscalar parallel time `. ˜ i (α(i) , T ) Hence the conditions of proposition 2 hold, and the square-free factorisation of G O(1) can be computed for every α ∈ D in nonscalar sequential time O(dL + d ) and nonscalar parallel time O(` + log3 d). ¡ ¢ Finally, the construction of the polynomial B P˜i (α(i) , T ) defined in section 7 (whose roots are the residues of P˜i (α(i) , T ) at all its poles) is performed by means of O(dL + dO(1) ) 3 (nonscalar) arithmetic operations ¡ ¢ in O(` + log d) nonscalar parallel steps. (i) We compute B P˜i (α , T ) and check whether all their nonprincipal coefficients vanish in parallel for i = 1, . . . , n by means of a correct test sequence of polynomials produced by straight-line programs of nonscalar size O(dL + dO(1) ) and nonscalar depth O(` + log3 d). Combining the results of proposition 3 and corollary 1 the theorem follows. 2 It remains to answer the third question: the problem of computing a rational potential whenever it is possible. We solve it in the following theorem: Theorem 8 If (P1 , . . . , Pn ) allows a rational potential, there exists an arithmetic network computing and evaluating Q of nonscalar sequential complexity LndO(1) and nonscalar parallel complexity O(` + log3 d). Z αi Proof With notations as in theorem 3.4, fix α ∈ Qn . We have to compute P˜i (α(i) , T )dT . 0

As it has shown, there exists a arithmetic network of nonscalar size O(dL + dO(1) ) and ˜ i and the nonscalar depth (` + log3 d) which computes the mixed representation of F˜i , G (i) ˜ square-free factorisation of Gi (α , T ).

16

˜ i (α(i) , T ) = Let G

d Y ¡ ¢j ˜ i (α(i) , T ). gj (α(i) , T ) the square-free factorisation of G j=1

By lemma 4, in (nonscalar) sequential time LdO(1) and parallel time O(` + log3 d) we can obtain the decomposition in partial fraction of P˜i (α(i) , T ) with respect to the square-free facd X Aj (α(i) , T ) (i) (i) (i) ˜ ˜ torisation of Gi (α , T ), that is, the expression Pi (α , T ) = A0 (α , T ) + ¡ ¢j (i) j=1 gj (α , T ) ¡ ¢ ¡ ¢ where Aj ∈ Q(X1 , . . . , Xi−1 )[T ] and degT Aj (α(i) , T ) < degT gj (α(i) , T ) . Aj (α(i) , T ) Further, the rational decomposition of every ¡ ¢j , i.e. gj (α(i) , T ) Z

Bj (α(i) , T ) Aj (α(i) , T ) ¡ ¢j dT = ¡ ¢j−1 + gj (α(i) , T ) gj (α(i) , T )

Z

Cj (α(i) , T ) dT gj (α(i) , T )

can be computed with (nonscalar) sequential complexity LdO(1) and parallel complexity O(` + log3 d). Since our integral is rational by lemma 5 and gj (α(i) , T ) is square-free, its logarithmic Z Cj (α(i) , T ) part, dT , is null. Therefore, gj (α(i) , T ) Z

αi

Z

αi

P˜i (α(i) , T )dT =

0

A0 (α(i) , T )dT +

0

d X j=1

Bj (α(i) , T ) ¡ ¢j−1 gj (α(i) , T )

αi 0

Z αi First we observe that the computation of A0 (α(i) , T )dT does not have cost in the 0 Z αi nonscalar model. Besides, the integrals P˜i (α(i) , T )dT can be computed in parallel 0

for every i = 1, . . . , n. Thus, this yields an arithmetic network with nonscalar sequential complexity LndO(1) and nonscalar parallel complexity O(` + log3 d) which computes Q(X1 , . . . , Xn ). 2

5

Conclusions

We have shown that integration on straight-line programs can be performed by means of well parallelizable algorithms if one allows polynomial-size degree growth of bounds. It can be seen that the mixed representation is a key model. Let us observe that this is not an improvement of the known complexity bounds in the one-variate case, because there our methods consist on interpolation and classical techniques of dense integration. Nevertheless, the most interesting examples arise in this case.

17

Acknowledgments The author is deeply indebted to Joos Heintz, for having strongly inspirated this work.

References [1] Baur W., Strassen V.: The complexity of partial derivatives. Theoret. Comput. Sci. 22 (1982) 317–330. [2] Berkowitz S.J.: On computing the determinant in small parallel time using a small number of processors. Information Processing Letter 18 (1984) 147–150. [3] Borodin A., von zur Gathen J., Hopcroft J.: Fast parallel matrix and GCD computations. Information and Control 52 (1982) 241–256. [4] Bronstein M.: Formulas for series computations. Applied Algebra in Engineering, Communication and Computing, AAECC 2, Springer-Verlag (1992) 195–206. [5] Davenport J.H.: Int´egration Formelle. IMAG Reserch Report No.375 (1983) 1–23. [6] Fitchas N., Giusti M. and Smietanski F.: Sur la complexit´e du th´eor`eme des z´eros. Preprint Ecole Polytechnique Palaiseau (1992). [7] von zur Gathen J.: Parallel algorithms for algebraic problems. Proc. 13-th. Conf. MFCS, Springer LN Comput. Sci. 356 (1989) 269–300. [8] Heintz J.: On the computational complexity of polynomials and bilinear mappings. A survey. Applied Algebra, Algebraic Algorithms and Error–Correcting Codes, 5th Intern. Conf. AAECC-5, Menorca 1987, L. Huguet and A. Poli, eds., Springer LN Comput. Sci. 356 (1989) 269–300. [9] Heintz J. and Schnorr C.P.: Testing polynomials which are easy to compute, in : 12-th Ann. ACM Symp. Theory of Computing (1980) 262–280. [10] Ibarra O.H. and Moran S.: Probabilistic algorithms for deciding equivalence of straight-line programs. J.ACM 30, 1 (1983) 217–228. [11] Kaltofen E.: Greatest common divisors of polynomials given by Straight-line Programs. J.ACM 35 No. 1 (1988) 234–264. [12] Lang S.: Algebra. Adisson-Wesley Publ. Comp., Reading, Massachusetts (1969). [13] Stoss H.J.: On the representation of rational functions of bounded complexity. Theoret. Comput. Sci. 64 (1989) 1–13. [14] Strassen V.: Berechnung und Programm I. Acta Inform.1 (1972) 320–334. [15] Strassen V.: Vermeidung von Divisionen. J. reine u. angew. Math. vol. 264 (1973) 182–202.

18