Tellegen’s Principle into Practice

A. Bostan

Laboratoire STIX FRE CNRS 2341 � Ecole polytechnique

91128 Palaiseau, France

bostan@stix.polytechnique.fr

G. Lecerf

Laboratoire de Math� ematiques

UMR CNRS 8100 Universit� e de Versailles St-Quentin-en-Yvelines 45, avenue des � Etats-Unis 78035 Versailles, France

lecerf@math.uvsq.fr

� E. Schost

Laboratoire STIX FRE CNRS 2341 � Ecole polytechnique 91128 Palaiseau, France

schost@stix.polytechnique.fr

ABSTRACT

The transposition principle, also called Tellegen’s principle, is a set of transformation rules for linear programs. Yet, though well known, it is not used systematically, and few practical implementations rely on it. In this article, we propose explicit transposed versions of polynomial multiplication and division but also new faster algorithms for multipoint evaluation, interpolation and their transposes. We report on their implementation in Shoup’s NTL C++ library.

Categories and Subject Descriptors

F.2 [Theory of Computation]: Analysis of Algorithms and Problem Complexity; D.4.8 [Mathematics of Computing]: Mathematical Software

General Terms

Algorithm, Theory

Keywords

Computer algebra, interpolation, polynomial evaluation, Tellegen’s principle, transposition principle

1. INTRODUCTION

The transposition principle, sometimes referred to as Tellegen’s principle, asserts that a linear algorithm that performs a matrix-vector product can be transposed, producing an algorithm that computes the transposed matrix-vector product. Further, the transposed algorithm has almost the same complexity as the original one (see Section 3 for precise statements).

The following example illustrates this principle, using the computation graph representation. Taking x_{1}, x_{2} as input,

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro�t or commercial advantage and that copies bear this notice and the full citation on the �rst page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior speci�c permission and/or a fee. ISSAC’03, August 3�6, 2003, Philadelphia, Pennsylvania, USA. Copyright 2003 ACM 1-58113-641-2/03/0008 ...$5.00.

it computes y_{1} = ax_{1} + bx_{2}, y_{2} = cx_{1} + dx_{2}; edges perform multiplications by the constant values a, b, c, d.

Reversing all arrows and exchanging vertices + and • yield the following graph:

Taking y_{1}, y_{2} as input, it computes the transposed map x_{1} = ay_{1} + cy_{2}, x_{2} = by_{1} + dy_{2} (see [16] for details).

Such transformation techniques originate from linear circuit design and analysis [1, 3, 19, 24] and were introduced in computer algebra in [7, 8, 12, 16]. Since then, there has been a recurrent need for transposed algorithms [2, 4, 6, 10, 13, 14, 21, 22, 23, 25]. Yet, the transposition principle in itself is seldom applied, and specific algorithms were often developed to circumvent its use, with the notable exceptions of [10, 22].

Contributions. In this paper, we detail several linear algorithms for univariate polynomials and their transposes: multiplication, quotient, remainder, evaluation, interpolation and exemplify a systematic use of Tellegen’s principle.

Our first contribution concerns univariate polynomial remaindering: we show that this problem is dual to extending linear recurrence sequences with constant coefficients. This clarifies the status of algorithms by Shoup [21, 23], which were designed as alternatives to the transposition principle: they are actually the transposes of well-known algorithms.

Our second contribution is an improvement (by a constant factor) of the complexities of multipoint evaluation and interpolation. This is done by designing a fast algorithm for transposed evaluation, and transposing it backwards. We also improve the complexity of performing several multipoint evaluations at the same set of points: discarding the costs of the precomputations, multipoint evaluation and interpolation now have very similar complexities.

Figure 1: Transposed/direct ratios

Finally, we demonstrate that the transposition principle is quite practical. We propose (still experimental) NTL [20] implementations of all algorithms mentioned here and their transposes. They show that the expected time ratio of 1 between an algorithm and its transpose is well respected in practice. The source code can be downloaded from

http://www.math.uvsq.fr/~lecerf.

Contents. Section 2 introduces the notation used in the rest of this article. In Section 3, we state the transposition principle in two distinct computation models; we discuss the memory requirement question that was raised in Kaltofen’s Open Problem 6 in [13].

In Section 4, we transpose polynomial multiplication. This was already done in [10], where transposed versions of Karatsuba’s and Fast Fourier Transform (FFT) multiplications were given. We transpose NTL Karatsuba’s multiplication, which was optimized to treat unbalanced situations. We also transpose plain, as well as FFT multiplication.

In Section 5, we transpose the operations of polynomial quotient and remainder and show the duality between polynomial remaindering and extending linear recurrence sequences.

In Section 6, we finally consider polynomial multipoint evaluation and interpolation. We give algorithms for these operations and their transposes that improve previously known algorithms by constant factors.

Implementation. Our algorithms are implemented in the C++ library NTL [20]. Figures 1 and 2 describe the behavior of our implementation; the computations were done over ℤ/pℤ, with p a prime of 64 bit length, and the times were measured on a 32 bit architecture.

Figure 1 shows the time ratios between direct and transposed algorithms for polynomial multiplication and polynomial remainder, for degrees up to 200.

• For multiplication, the horizontal axis gives the degree m of the input polynomials. Multiplication in NTL uses three different algorithms: plain multiplication for m ≤ 19, Karatsuba’s multiplication for 20 ≤ m < 79 and FFT for larger m. The same thresholds are used for the transposed versions. Note that Karatsuba’s

Figure 2: Transposed/direct ratios

multiplication is slightly slower in its transposed version.

• The division deals with polynomials of degrees 2m and m. NTL provides the plain remainder algorithm and a faster one based on Sieveking-Kung’s algorithm, with further optimizations. The threshold is 90, we used the same value for the transposed algorithms.

Figure 2 presents the ratios between the direct and transposed versions for our new fast multipoint evaluation and interpolation algorithms, with degree on the horizontal axis (Section 6 gives their complexity estimates).

Remark. Transposed algorithms were used for modular multiplication and power series inversion in NTL as soon as [22], yet only for FFT multiplication. More generally the transposed product is used in [10] to speed up some algorithms for power series.

2. DEFINITIONS AND NOTATION

Let R be a commutative ring with unity and n ≥ 0. By R[x]_{n} we denote the free R-module of polynomials of degree at most n. We use the monomial basis 1, x, …, x^{n} on R[x]_{n} to represent polynomials by vectors and linear maps by matrices. It will also be convenient to consider elements of R[x] as infinite sequences with finite support: we will write a polynomial a as ∑_{i≥0} a_{i}x^{i}, where almost all a_{i} vanish. The degree of a is denoted by deg(a), with deg(0) = −∞. For l ≥ 0 and h ≥ 0, we introduce the following maps

⌈.⌉^{h}, ⌊.⌋_{l}, [.]_{l}^{h} on R[x] and the power series ring R[[x]]:

⌈a⌉^{h} = ∑_{i=0}^{h−1} a_{i}x^{i}, ⌊a⌋_{l} = ∑_{i≥0} a_{i+l}x^{i}, [a]_{l}^{h} = ∑_{i=0}^{h−l−1} a_{i+l}x^{i}.

Observe that these maps satisfy the following relations:

a = ⌈a⌉^{h} + x^{h} ⌊a⌋_{h} , [a]_{l}^{h} = ⌈⌊a⌋_{l}⌉^{h−l} = ⌊⌈a⌉^{h}⌋_{l} .

We define the reversal endomorphism rev(n, .) of R[x]_{n} by rev(n, a) = ∑_{k=0}^{n} a_{n−k}x^{k}, for all a ∈ R[x]_{n}.

By ⌊q⌋ we denote the integer part of a rational q. Finally, we use a block matrix notation: 0_{h,l} denotes a h × l block filled with zeros and 1_{h} the h × h identity matrix.

3. TELLEGEN’S PRINCIPLE

Tellegen’s principle is usually stated using linear computation graphs [16, 24]. In this section, we first briefly recall it in terms of linear straight-line programs. We then define another computational model, and prove Tellegen’s principle in that setting. The first model can be thought as taking only time complexity into account, while the second one considers both time and space.

Linear Straight-Line Programs. Linear straight-line pro-

grams can be thought as “ordinary” straight-line programs, but with only linear operations, see [5, Chapter 13] for precise definitions. Their complexity is measured by their number of operations, the size, which reflects a time complexity. In this model, Tellegen’s principle can be formulated as:

Proposition 1. [5, Th. 13.20] Let ϕ : R^{n} → R^{m} be a linear map that can be computed by a linear straight-line program of size L and whose matrix in the canonical bases has z_{r} zero rows and z_{c} zero columns. Then the transposed map ϕ^{t} can be computed by a linear straight-line program of size L − n + m − z_{r} + z_{c}.

Let us remark that Theorem 13.20 in [5] actually expresses this result in terms of linear complexity of ϕ, that is, the minimum size of a linear straight-line program computing ϕ. A quick look at the proof shows that the statement in the proposition above is implicitely proved therein. We will use Proposition 1 in Section 6, as the complexity estimates in that section are given in the linear straight-line program model.

Another Computational Model. We now introduce a mea-

sure of space complexity. Informally speaking, our model derives from the random access memory model by restricting to linear operations. In this context, a program P is given by:

• A finite set D of registers, with two distinguished subsets (non necessarily disjoint), I for the input and O for the output: before execution I contains the input values, at the end O contains the output.

• A finite sequence of instructions (ϕ_{i})_{i∈{1,... ,L}} of length L. Each instruction is a linear endomorphism of R^{D} of the following type, with p, q in D, a in R and f in R^{D}:

– p + = q denotes the map that sends f to the function g defined by g(r) = f(r) if r ≠ p and g(p) = f(p) + f(q).

– p ∗ = a denotes the map that sends f to the function g defined by g(r) = f(r) if r ≠ p and g(p) = af(p).

– p = 0 denotes the map that sends f to the function g defined by g(r) = f(r) if r ≠ p and g(p) = 0.

Let I (resp. P) be the injection R^{I} → R^{D} (resp. the projection R^{D} → R^{O}). Then we say that the program P computes the linear map Φ : R^{I} → R^{O} defined by P ∘ ϕ_{L} ∘ ϕ_{L−1} ∘ ⋯ ∘ ϕ_{1} ∘ I.

We can now define the transpose P^{t}(D^{t}, I^{t}, O^{t}) of P:

• The set of registers of P^{t} is still D^{t} = D but we let I^{t} = O and O^{t} = I: input and output are swapped.

• The instructions of P^{t} are ϕ_{L}^{t}, ϕ_{L−1}^{t}, …, ϕ_{1}^{t}.

Let us verify that P^{t} is well-defined. We examine the transpose of each type of instructions:

• The transpose ϕ^{t} of an instruction ϕ of type p += q is the instruction q += p.

• The last two instructions are symmetric maps.

The equality (ϕ_{L} ∘ ϕ_{L−1} ∘ ⋯ ∘ ϕ_{1})^{t} = ϕ_{1}^{t} ∘ ϕ_{2}^{t} ∘ ⋯ ∘ ϕ_{L}^{t} then shows that P^{t} computes the linear map Φ^{t}. This yields the following form of Tellegen’s principle:

Proposition 2. According to the above notation, let ϕ : R^{n} → R^{m} be a linear map that can be computed by a program with D registers and L instructions. Then ϕ^{t} can be computed by a program with D registers and L instructions.

As an example, consider 5 registers, x_{1}, x_{2} for input, y_{1}, y_{2} for output and r for temporaries. Given a, b, c, d in R, consider the instructions: y_{1} += x_{1}, y_{1} ∗= a, r += x_{2}, r ∗= b, y_{1} += r, r = 0, y_{2} += x_{1}, y_{2} ∗= c, r += x_{2}, r ∗= d, y_{2} += r. We let the reader check that the linear map computed by this program is the 2 × 2 matrix-vector product presented in the introduction. The transposed program is r += y_{2}, r ∗= d, x_{2} += r, y_{2} ∗= c, x_{1} += y_{2}, r = 0, r += y_{1}, r ∗= b, x_{2} += r, y_{1} ∗= a, x_{1} += y_{1}.

We could have used more classical instructions such as p = q + r and p = aq, a ∈ R. We let the reader check that such programs can be rewritten in our model within the same space complexity but a constant increase of time complexity. For such programs, Tellegen’s principle would be stated with no increase in space complexity, but a constant increase in time complexity. This is why we use straight-line programs for our complexity estimates in Section 6.

Open Problem 6 in [13] asks for a transposition theorem without space complexity swell. The above Proposition sheds new light on this problem: In the present computational model, it is immediate to observe that memory consumption is left unchanged under transposition.

Comments. The above models compute functions of fixed input and output size: in the sequel, we will actually write families of programs, one for each size. We also use control instructions as for and if, and calls to previously defined subroutines. Last, we will consider algorithms mixing linear and non-linear precomputations; the transposition principle leaves the latter unchanged.

4. POLYNOMIAL MULTIPLICATION

In this section a is a fixed polynomial of degree m. For an integer n ≥ 0, consider the multiplication map by a:

\begin{aligned}mul(a, .) : & R[x]_{n} & → & R[x]_{m+n} \\ & b & ↦ & ab. \end{aligned}

The transposed map is denoted by mul^{t}(n, a, .); to write our pseudo-code, it is necessary to consider n as an argument of the transposed function. The next subsection shows that this map is:

\begin{aligned}mul^{t}(n, a, .) : & R[x]_{m+n} & → & R[x]_{n} \\ & c & ↦ & [rev(m, a)c]_{m}^{n+m+1} . \end{aligned}

We adopt the following convention: if deg(c) > m + n then mul^{t}(n, a, c) returns an error.

The above formula explains why the transposed multiplication is also called the middle product [10]. Performing this operation fast is the first task to accomplish before transposing higher level algorithms. Observe that computing mul^{t}(n, a, c) by multiplying rev(m, a) by c before extracting the middle part requires to multiply two polynomials of degrees m and m + n.

Tellegen’s principle implies that this transposed computation can be performed for the cost of the multiplication of two polynomials of degrees m and n only. In what follows, we make this explicit for the plain, Karatsuba and Fast Fourier Transform multiplications.

4.1 Plain multiplication

In the canonical monomial bases, the matrix of the linear map mul(a, .) : R_{n}[x] → R_{m+n}[x] is the following Toeplitz matrix T with m + n + 1 rows and n + 1 columns:

T = [ \begin{matrix}a_{0} & 0 & … & 0 \\ a_{1} & a_{0} & ... & ⋮ \\ ⋮ & a_{1} & ... & 0 \\ ⋮ & ⋮ & ... & a_{0} \\ a_{m} & ⋮ & ... & a_{1} \\ 0 & a_{m} & ... & ⋮ \\ ⋮ & ... & ... & ⋮ \\ 0 & … & 0 & a_{m} \end{matrix} ] .

Basically, the plain multiplication algorithm corresponds to performing the matrix vector product c = T b naively using the following sequence of instructions:

c ← 0; for i from 0 to m + n do

for j from max(0, i − n) to min(m, i) do c_{i} ← c_{i} + a_{j}b_{i−j};

Reversing the computation flow, b = mul^{t}(n, a, c) is computed by the following program:

b ← 0; for i from m + n downto 0 do

for j from min(m, i) downto max(0, i − n) do b_{i−j} ← b_{i−j} + a_{j}c_{i};

Observing that the entries of b are given by:

b_{i} = ∑_{j=i}^{i+m} a_{j−i}c_{j}, i ∈ {0, …, n},

we deduce that the map mul^{t}(n, a, .) reformulates in terms of the middle product [10] by the relation

mul^{t}(n, a, c) = [rev(m, a)c]_{m}^{n+m+1} ,

as announced above.

4.2 Karatsuba’s algorithm

The next paragraphs are devoted to the transposition of Karatsuba’s multiplication algorithm. Concretely, we consider the NTL implementation and present the transpose

Figure 3: KarMul(a, b)

m ← deg(a); n ← deg(b); if n = 0 then return b_{0}a; if m ≤ 0 then return a_{0}b; μ ← ⌊m/2⌋ + 1; ν ← ⌊n/2⌋ + 1; if μ > n then

u ← KarMul(⌈a⌉^{μ} , b); v ← KarMul(⌊a⌋_{μ} , b); return u + x^{μ}v;

if ν > m then u ← KarMul(a, ⌈b⌉^{ν}); v ← KarMul(a, ⌊b⌋_{ν}); return u + x^{ν}v;

h ← max(μ, ν); r ← KarMul(⌈a⌉^{h} , ⌈b⌉^{h}); s ← KarMul(⌊a⌋_{h} , ⌊b⌋_{h}); t ← KarMul(⌈a⌉^{h} + ⌊a⌋_{h} , ⌈b⌉^{h} + ⌊b⌋_{h}); return r + x^{h}(t − s − r) + x^{2h}s;

version we made of it. Figures 3 and 4 describe Karatsuba’s algorithm and its transpose. As for the mul^{t} function, the transposed Karatsuba multiplication takes n as an additional argument.

In order to prepare the proof of Algorithm TKarMul we first decompose Algorithm KarMul into linear maps. Following Algorithm KarMul we enter the procedure with polynomials a of degree m and b of degree n. We let μ = ⌊m/2⌋+1, ν = ⌊n/2⌋ + 1 and h = max(μ, ν) and distinguish three exclusive cases: Normal case: (m < n and ν ≤ m) or (m ≥ n and μ ≤ n). Note that deg(⌊a⌋_{h}) = m − h, let ρ = deg(⌈a⌉^{h}) and λ = deg(⌈a⌉^{h} + ⌊a⌋_{h}). Let M_{l} : R[x]_{h−1} → R[x]_{ρ+h−1}, M_{h} : R[x]_{n−h} → R[x]_{m+n−2h} and M_{f} : R[x]_{h−1} → R[x]_{λ+h−1} be the linear maps of multiplication by resp. ⌈a⌉^{h}, ⌊a⌋_{h} and ⌈a⌉^{h} + ⌊a⌋_{h}. Noticing that n − h ≤ h − 1, we construct the following block matrix M, where we take r = m + n and q = ρ + h:

M = [ \frac{\begin{matrix}M_{l} & 0_{q,n−h+1} \\ 0_{r−2h+1,h} & M_{h} \end{matrix}}{M_{f} · [ 1_{h} \frac{1_{n−h+1}}{0_{2h−n−1,n−h+1}} ]} ] .

Remark that (r, s, t) = Mb, when considering that r ∈ R[x]_{ρ+h−1}, s ∈ R[x]_{m+n−2h} and t ∈ R[x]_{λ+h−1}. In order to recover the product c = ab it remains to compute c = (N_{+} − N_{−})(r, s, t), with

N_{+} = [ \frac{1_{ρ+h}}{0_{r+1−q,q}}  \frac{0_{2h,r−2h+1}}{1_{r−2h+1}}  \begin{matrix}0_{h,λ+h} \\ 1_{λ+h} \\ 0_{r+1−λ−2h,λ+h} \end{matrix} ] ,

N_{−} = [ \begin{matrix}0_{h,ρ+h} \\ 1_{ρ+h} \\ 0_{r−ρ−2h+1,ρ+h} \end{matrix}  \begin{matrix}0_{h,r−2h+1} \\ 1_{r−2h+1} \\ 0_{h,r−2h+1} \end{matrix}  0_{r+1,λ+h} ] .

Figure 4: TKarMul(n, a, c)

m ← deg(a); if deg(c) > m + n then Error; if n = 0 then return ∑_{k=0}^{m} a_{k}c_{k}; if m ≤ 0 then return a_{0}c; μ ← ⌊m/2⌋ + 1; ν ← ⌊n/2⌋ + 1; ρ ← deg(⌈a⌉^{h}); if μ > n then

u ← TKarMul(n, ⌈a⌉^{μ} , ⌈c⌉^{n+ρ}); v ← TKarMul(n, ⌊a⌋_{μ} , ⌊c⌋_{μ}); return u + v;

if ν > m then u ← TKarMul(ν − 1, a, ⌈c⌉^{m+ν}); v ← TKarMul(n − ν, a, ⌊c⌋_{ν}); return u + v;

h ← max(μ, ν); λ ← deg(⌈a⌉^{h} + ⌊a⌋_{h}); r ← TKarMul(h − 1, ⌈a⌉^{h} , ⌈c⌉^{ρ+h} − [c]_{h}^{ρ+2h}); s ← TKarMul(n − h, ⌊a⌋_{h} , ⌊c⌋_{2h} − [c]_{h}^{m+n−h+1}); t ← TKarMul(h − 1, ⌈a⌉^{h} + ⌊a⌋_{h} , [c]_{h}^{λ+2h}); return r + t + x^{h}(s + ⌈t⌉^{n−h+1});

Degenerate case 1: μ > n. We let ρ = deg(⌈a⌉^{μ}) and consider u ∈ R[x]_{ρ+n}, v ∈ R[x]_{m+n−μ} so that we have ab = N_{1}(u, v), where N_{1} is:

N_{1} = [ \frac{1_{ρ+n+1}}{0_{m−ρ,ρ+n+1}} \frac{0_{μ,r−μ+1}}{1_{r−μ+1}} ] .

Degenerate case 2: ν > m. We consider u ∈ R[x]_{m+ν−1}, v ∈ R[x]_{m+n−ν} so that we have ab = N_{2}(u, v), with

N_{2} = [ \frac{1_{m+ν}}{0_{n−ν+1,m+ν}} \frac{0_{ν,r−ν+1}}{1_{r−ν+1}} ] .

Proposition 3. Algorithm TKarMul is correct.

Proof. We use the notation of Figure 4 and distinguish the same three cases. The normal case follows from the equalities

N_{+}^{t}c = (⌈c⌉^{ρ+h} , ⌊c⌋_{2h} , [c]_{h}^{λ+2h}),
N_{−}^{t}c = ([c]_{h}^{ρ+2h} , [c]_{h}^{m+n−h+1} , 0).

In Degenerate cases 1 and 2, we respectively compute

N_{1}^{t}c = (⌈c⌉^{n+ρ} , ⌊c⌋_{μ}) and N_{2}^{t}c = (⌈c⌉^{m+ν} , ⌊c⌋_{ν}).

□

Some refinements may be done when implementing Algorithm TKarMul. First observe that it saves memory to compute ⌈c⌉^{m+n−h+1} − ⌊c⌋_{h} in order to deduce r and −s and then one propagates this change of sign by returning r +t+x^{h}(⌈t⌉^{n−h+1} −s). Another observation is that Karatsuba’s multiplication can be done in a different manner using the identity ⌈a⌉^{h} ⌊b⌋_{h} + ⌊a⌋_{h} ⌈b⌉^{h} = ⌈a⌉^{h} ⌈b⌉^{h} + ⌊a⌋_{h} ⌊b⌋_{h} − (⌈a⌉^{h} − ⌊a⌋_{h})(⌈b⌉^{h} − ⌊b⌋_{h}). When transposing this slightly different version, we obtain the middle product algorithm presented in [10].

4.3 The Fast Fourier Transform

Multiplication algorithms using the Fast Fourier Transform are quite easy to transpose since the matrices involved are symmetric. We only give a brief presentation and refer to [9] for more details about the discrete Fourier transform. Let l ∈ ℕ such that m + n + 1 ≤ 2^{l} and that R contains a primitive 2^{l}-th root ω of unity. The discrete Fourier transform DFT(ω, a) of the polynomial a is the vector

(a(1), a(ω), …, a(ω^{2^{l}−1})) ∈ R^{2^{l}}.

Let DFT^{−1}(ω, .) : R^{2^{l}} → R[x]_{2^{l}−1} be the inverse map of DFT(ω, .) and H the diagonal matrix of diagonal DFT(ω, a). Then we have the equality

ab = DFT^{−1}(ω, H DFT(ω, b)).

Since DFT(ω, .) and H are symmetric, we deduce that:

mul^{t}(n, a, c) = ⌈DFT(ω, H DFT^{−1}(ω, c))⌉^{n+1} ,

for any polynomial c of degree at most m + n. Letting H̃ be the diagonal DFT(ω, rev(m, a)) matrix and using

DFT^{−1}(ω, .) = \frac{1}{2^{l}} DFT(ω^{−1}, .),

we deduce the equalities

mul^{t}(n, a, c) = ⌈DFT^{−1}(ω, ω^{−m}H̃ DFT(ω, c))⌉^{n+1}
= [DFT^{−1}(ω, H̃ DFT(ω, c))]_{m}^{m+n+1} ,

which can also be obtained from the middle product formulation.

5. POLYNOMIAL DIVISION

We come now to the transposition of the Euclidean division. In this section we are given a polynomial a ≠ 0 of degree m whose leading coefficient a_{m} is invertible. For a polynomial b, we write the division of b by a as b = aq + r with deg(r) < m and define the maps

quo(a, .) : R[x]_{n} →  R[x]_{n−m}
b  ↦  q,
rem(a, .) : R[x]_{n} →  R[x]_{m−1}
b  ↦  r,
quorem(a, .) : R[x]_{n} → R[x]_{n−m} × R[x]_{m−1}
b  ↦  (q, r).

The transposed operations are written quorem^{t}(n, a, q, r), quo^{t}(n, a, q) and rem^{t}(n, a, r), with the convention that these functions return an error if deg(r) ≥ deg(a) or deg(q) > n − deg(a).

The next paragraphs are devoted to transpose the quorem map through the plain and Sieveking-Kung’s division algorithms. We will prove that the transposed remainder is given by:

\begin{aligned}rem^{t}(n, a, .) : & R^{m} & → & R^{n+1} \\ & (r_{0}, …, r_{m−1}) & ↦ & (b_{0}, …, b_{n}), \end{aligned}

where the b_{j} are defined by the linear recurrence

(∗) b_{j} = − \frac{1}{a_{m}} (a_{m−1}b_{j−1} + ⋯ + a_{0}b_{j−m}), m ≤ j

with initial conditions b_{j} = r_{j}, j ∈ {0, …, m − 1}, so that linear recurrence sequence extension is dual to remainder computation.

5.1 Plain division

We enter the plain division procedure with the two polynomials a and b and compute q and r by

q ← 0; r ← b; for i from 0 to n − m do

q ← xq + r_{n−i}/a_{m}; r ← r − r_{n−i}/a_{m}x^{n−m−i}a;

To transpose this algorithm we introduce the sequences

q^{{i}} ∈ R[x]_{i−1}, r^{{i}} ∈ R[x]_{n−i} i = 0, …, n − m + 1.

They are defined by q^{{0}} = 0, r^{{0}} = b and for i ≥ 1

q^{{i+1}} = xq^{{i}} + r_{n−i}^{{i}}/a_{m},
r^{{i+1}} = r^{{i}} − r_{n−i}^{{i}}/a_{m}x^{n−m−i}a,

so that the relation b = ax^{n−m+1−i}q^{{i}} + r^{{i}} holds. For i = n − m + 1, we have q = q^{{n−m+1}} and r = r^{{n−m+1}}.

In order to formulate the algorithm in terms of linear maps we introduce v^{{i}} = (q^{{i}}, r^{{i}}) ∈ R^{n+1}. Then we have v^{{i+1}} = Mv^{{i}}, where

M = [ \begin{matrix}0_{1,n} & 1/a_{m} \\ 1_{n} & 0_{n−m,1} ; −a_{0}/a_{m} ; ⋮ ; −a_{m−1}/a_{m} \end{matrix} ] .

Now reverse the flow of calculation: we start with a vector v^{{n−m+1}} = (q^{{n−m+1}}, r^{{n−m+1}}) ∈ R^{n+1} with q^{{n−m+1}} ∈ R[x]_{n−m} and r^{{n−m+1}} ∈ R[x]_{m−1}. Then we compute v^{{i}} = M^{t}v^{{i+1}} by the formulae

q^{{i}} = ⌊q^{{i+1}}⌋_{1} ,

r^{\{i\}} = r^{\{i+1\}} + \frac{1}{a_{m}} x^{n−i}(q_{0}^{\{i+1\}} − ∑_{j=0}^{m−1} a_{m−1−j}r_{n−i−j−1}^{\{i+1\}}).

We deduce the following transposed algorithm for computing b = quorem^{t}(n, a, q, r):

b ← r; for i from m to n do

b_{i} ← (q_{i−m} − ∑_{j=0}^{m−1} a_{m−1−j}b_{i−j−1}) /a_{m};

The maps quo(a, .) and rem(a, .) can be obtained by composing a projection after quorem(a, .) but in practice it is better to implement specific optimized procedures for each. It is easy to implement these optimizations; we refer to our NTL implementation for details.

Since b = rem^{t}(n, a, r) = quorem^{t}(n, a, 0, r), the coefficients of b satisfy the linear recurrence relation with constant coefficients (∗), as announced above.

5.2 Sieveking-Kung’s division

Sieveking-Kung’s algorithm is based on the formula

rev(n, b) = rev(n − m, q)rev(m, a) + x^{n−m+1}rev(m − 1, r), see [9] for details. This yields

q = rev(n − m, ⌈mul(α, ⌈rev(n, b)⌉^{n−m+1})⌉^{n−m+1}),

where α = ⌈rev(m, a)^{−1}⌉^{n−m+1}. Then we deduce r from q using r = b − aq.

Transposing these equalities, we see that for q ∈ R[x]_{n−m}

quo^{t}(n, a, q) = rev(n, mul^{t}(n − m, α, rev(n − m, q))).

For r ∈ R[x]_{m−1}, it follows:

rem^{t}(n, a, r) = r − quo^{t}(n, a, mul^{t}(n − m, a, r)).

Let s = rem^{t}(n, a, r) and p = mul^{t}(n − m, a, r). Using the middle product formula to express this last quo^{t} expression yields

s = r − rev(n, [rev(n − m, α)rev(n − m, p)]_{n−m}^{2(n−m)+1}),

which simplifies to

s = r − rev(n, rev(n − m, ⌈αp⌉^{n−m+1})).

Last we obtain

rem^{t}(n, a, r) = r − x^{m} ⌈αmul^{t}(n − m, a, r)⌉^{n−m+1} ,

which can also be rewritten this way, using the middle product formula again:

rem^{t}(n, a, r) = r − x^{m} ⌈α [rev(m, a)r]_{m}^{n+1}⌉^{n−m+1} .

This actually coincides with the algorithm given in [21] for extending linear recurrence sequences.

Remark. More generally, the following formula holds:

quorem^{t}(n, a, q, r) = r−x^{m} ⌈α([rev(m, a)r]_{m}^{n+1} − q)⌉^{n−m+1} .

5.3 Modular multiplication

As a byproduct, our algorithms enable to transpose the modular multiplication. Given a monic polynomial a of degree m with invertible leading coefficient and a polynomial b of degree at most m − 1, consider the composed map

\begin{matrix}R[x]_{m−1} & → & R[x]_{2m−2} & → & R[x]_{m−1} \\ c & ↦ & bc & ↦ & bc mod a. \end{matrix}

An ad hoc algorithm for the transpose map is detailed in [23]. Using our remark on the transpose of polynomial remaindering, it is seen to essentially coincide with the one obtained by composing the algorithms for transposed multiplication and transposed remainder. A constant factor is lost in [23], as no middle product algorithm is used to transpose the first map; this was already pointed out in [10].

6. TRANSPOSED VANDERMONDE

We now focus on algorithms for Vandermonde matrices; we assume that R = k is a field and consider m + 1 pairwise distinct elements a_{0}, …, a_{m} in k. Even if not required by the algorithms, we take m = 2^{l} − 1, with l ∈ ℕ, in order to simplify the complexity estimates.

The Vandermonde matrix V_{a} is the square matrix:

V_{a} = [ \begin{matrix}1 & a_{0} & a_{0}^{2} & … & a_{0}^{m} \\ 1 & a_{1} & a_{1}^{2} & … & a_{1}^{m} \\ ⋮ & ⋮ & ⋮ & & ⋮ \\ 1 & a_{m} & a_{m}^{2} & … & a_{m}^{m} \end{matrix} ] .

If b is in k[x]_{m} then V_{a}b is the vector (b(a_{0}), …, b(a_{m})). This computation is commonly referred to as multipoint evaluation. The inverse problem is interpolation: for c ∈ k^{m+1}, V_{a}^{−1}c corresponds to the polynomial b ∈ k[x]_{m} satisfying b(a_{i}) = c_{i} for i ∈ {0, …, m}.

We first recall the notion of subproduct tree and refer to [9, §10.1] for details and historical notes. Then we design a fast algorithm for transposed evaluation. We transpose it backwards to obtain an improved evaluation algorithm, and deduce similar improvements for interpolation and its transpose.

For complexity analysis we use the straight-line program model [5, Chapter 4]. The function M(n) denotes the complexity of multiplying a polynomial of degree less than n by a fixed polynomial of degree less than n, restricting to linear straight-line programs only. This way Proposition 1 implies that the complexity of the transpose of this multiplication is at most M(n) + n operations in k. Restricting to linear straight-line programs is actually not embarrassing since all known multiplication algorithms [9] for the straight-line program model fit into this setting.

As in [9, §8.3] we assume that M(n_{1} + n_{2}) ≥ M(n_{1}) + M(n_{2}) for any positive integers n_{1} and n_{2}. For the sake of simplicity we also assume that n log(n) ∈ O(M(n)).

6.1 Going up the subproduct tree

A common piece of the following algorithms is the computation of the subproduct tree T with leaves x − a_{0}, x − a_{1}, …, x − a_{m}. It is defined recursively, together with the sequence of integers m_{i}, by

T_{0,j} = x − a_{j}, for j ∈ {0, …, m}, m_{0} = m + 1,

and for i ≥ 1 by

T_{i,j} = T_{i−1,⌊j/2⌋}T_{i−1,⌊j/2⌋+1}, for j < h_{i} = ⌊m_{i}/2⌋ .

If m_{i} = 2h_{i} + 1 we let T_{i,h_{i}} = T_{i−1,m_{i}−1} and m_{i+1} = h_{i} + 1, otherwise we just let m_{i+1} = h_{i}. Let d be the smallest integer such that m_{d} = 1; we get T_{d,0} = ∏_{j=0}^{m}(x−a_{j}). By [9, Exercise 10.3], T can be computed within 1/2M(m) log(m) + O(m log(m)) operations.

The following algorithm describes our basic use of T. On input c = (c_{0}, …, c_{m}), it computes the polynomial b = ∑_{j=0}^{m} c_{j} \frac{T_{d,0}}{x−a_{j}}.

UpTree(c)

b ← c; for i ← 0 to d − 1 do

for j ← 0 to h_{i} − 1 do

b_{j} ← T_{i,2j+1}b_{2j} + T_{i,2j}b_{2j+1}; if m_{i} = 2h_{i} + 1 then b_{h_{i}} ← b_{m_{i}−1};

return b_{0};

We deduce its transpose as follows:

TUpTree(b)

c_{0} ← b; for i ← d − 1 downto 0 do

if m_{i} = 2h_{i} + 1 then c_{m_{i}−1} ← c_{h_{i}}; for j ← h_{i} − 1 downto 0 do

n ← deg(T_{i+1,j}) − 1; c_{2j+1} ← mul^{t}(n, T_{i,2j}, c_{j}); c_{2j} ← mul^{t}(n, T_{i,2j+1}, c_{j});

return c; On input b ∈ k[x]_{m}, it outputs the coefficients of x^{m} in the polynomial products rev(m, b) \frac{T_{d,0}}{x−a_{j}}, for j = 0, …, m.

Once T is computed, using Proposition 1 and [9, Th. 10.10], both algorithms require M(m) log(m)+O(m log(m)) operations.

6.2 Multipoint evaluation

We first treat the transposed problem. Let c_{0}, …, c_{m} be in k. A direct computation shows that the entries of b = V_{a}^{t}c are the first m + 1 coefficients of the Taylor expansion of S(x) =

∑_{j=0}^{m} \frac{c_{j}}{1 − a_{j}x} = \frac{1}{rev(m + 1, T_{d,0})} ∑_{j=0}^{m} \frac{c_{j}rev(m + 1, T_{d,0})}{1 − a_{j}x} .

The last sum is obtained by computing UpTree(c) and reversing the result. Computing the Taylor expansion of S requires one additional power series inversion and one multiplication.

α ← 1/rev(m + 1, T_{d,0}) mod x^{m+1}; s ← UpTree(c); t ← rev(m, s); b ← mul(α, t) mod x^{m+1};

We deduce the following algorithm for evaluating a polynomial b:

α ← 1/rev(m + 1, T_{d,0}) mod x^{m+1}; t ← mul^{t}(m, α, b); s ← rev(m, t); c ← TUpTree(s);

By the above complexity results, these algorithms require 3/2M(m) log(m)+O(M(m)) operations, since the additional operations have negligible cost. We gain a constant factor on the usual algorithm of repeated remaindering, of complexity 7/2M(m) log(m) + O(M(m)), see [18] or [9, Exercise 10.9]. We also gain on 13/6M(m) log(m) + O(M(m)) given in [18, §3.7] for base fields k allowing Fast Fourier Transform in k[x].

Moreover, if many evaluations at the same set of points a_{i} have to be performed, then all data depending only on the evaluation points (the tree T and α) may be precomputed and stored, and the cost of evaluation drops to essentially M(m) log(m)+O(M(m)). This improves [9, Exercise 10.11] by a factor of 2.

As a remark, our interpretation of the transposed remainder shows that the algorithm in [21] for transposed evaluation is the exact transposition of the classical evaluation algorithm as given in [9, §10.1], so ours is faster. The algorithm of [6] is still slower by a constant factor, since it uses an interpolation routine.

6.3 Interpolation

The fast interpolation algorithm c = V_{a}^{−1}b proceeds this way, see [9, §10.2]:

p ← dT_{d,0}/dx; z ← (p(a_{0}), …, p(a_{m})); c ← (b_{0}/z_{0}, …, b_{m}/z_{m}); c ← UpTree(c);

Here, z is computed using fast multipoint evaluation. At the end, c contains the interpolating polynomial. Reversing the computation flow, we get the transposed algorithm for computing b = (V_{a}^{−1})^{t}c: p ← dT_{d,0}/dx; z ← (p(a_{0}), …, p(a_{m})); b ← TUpTree(c); b ← (b_{0}/z_{0}, …, b_{m}/z_{m});

Using the results given above, these algorithms require 5/2M(m) log(m) + O(M(m)) operations. This improves again the complexity results of [9].

The algorithm of [15] for transposed interpolation does the same precomputation as ours, but the call to TUpTree is replaced by (mainly) a multipoint evaluation. Using precomputations and our result on evaluation, that algorithm also has complexity 5/2M(m) log(m) + O(M(m)).

7. CONCLUSION, FUTURE WORK

A first implementation of our improved evaluation algorithm gains a factor of about 1.2 to 1.5 over the classical one [9, §10.1], for degrees about 10000. But let us mention that the crossover point between repeating Horner’s rule and the classical fast evaluation we have implemented is about 32. In consequence, it seems interesting to explore the constant factors hidden behind the above quantities O(M(m)).

Our results of Section 6 can be generalized to solve the simultaneous modular reduction and the Chinese remainder problems faster than in [9, Chapter 10]. Theoretical and practical developments related to these problems are work in progress.

Transposed computations with Vandermonde matrices are the basis of fast multiplication algorithms for sparse and dense multivariate polynomials [2, 6, 25] and power series [11, 17]. Up to logarithmic factors, over a base field of characteristic zero, these multiplications have linear complexities in the size of the output. Their implementation is the subject of future work.

Acknowledgments. We thank G. Hanrot, M. Quercia and P. Zimmermann for useful discussions and the anonymous referees for their constructive remarks.

8. REFERENCES

[1] A. Antoniou. Digital Filters: Analysis and Design. McGraw-Hill Book Co., 1979.

[2] M. Ben-Or and P. Tiwari. A deterministic algorithm for sparse multivariate polynomial interpolation. In 20th Annual ACM Symposium on the Theory of Computing, pages 301–309, 1988.

[3] J. L. Bordewijk. Inter-reciprocity applied to electrical networks. Appl. Sci. Res. B., 6:1–74, 1956.

[4] A. Bostan, B. Salvy, and É. Schost. Fast Algorithms for Zero-Dimensional Polynomial Systems Using Duality. Applicable Algebra in Engineering, Communication and Computing. To appear.

[5] P. Bürgisser, M. Clausen, and A. Shokrollahi. Algebraic Complexity Theory. Springer, 1997.

[6] J. Canny, E. Kaltofen, and L. Yagati. Solving systems of non-linear polynomial equations faster. In Proceedings of ISSAC’89, pages 121–128. ACM, 1989.

[7] C. M. Fiduccia. On obtaining upper bounds on the complexity of matrix multiplication. In Complexity of computer computations (Proc. Sympos., IBM Thomas J. Watson Res. Center, Yorktown Heights, N.Y., 1972), pages 31–40. Plenum, New York, 1972.

[8] C. M. Fiduccia. On the algebraic complexity of matrix multiplication. PhD thesis, Brown Univ., Providence, RI, Center Comput. Inform. Sci., Div. Engin., 1973.

[9] J. von zur Gathen and J. Gerhard. Modern computer algebra. Cambridge University Press, 1999.

[10] G. Hanrot, M. Quercia, and P. Zimmermann. The middle product algorithm, I. Speeding up the division and square root of power series. Technical report, RR INRIA 4664, 2002.

[11] J. van der Hoeven. Relax, but don’t be too lazy. Journal of Symbolic Computation, 34(9):479 – 542, 2002.

[12] J. Hopcroft and J. Musinski. Duality applied to the complexity of matrix multiplication and other bilinear forms. SIAM Journal on Computing, 2:159–173, 1973.

[13] E. Kaltofen. Challenges of symbolic computation: my favorite open problems. Journal of Symbolic Computation, 29(6):891–919, 2000.

[14] E. Kaltofen and V. Shoup. Subquadratic-time factoring of polynomials over finite fields. Mathematics of Computation, 67(223):1179–1197, 1998.

[15] E. Kaltofen and L. Yagati. Improved sparse multivariate polynomial interpolation algorithms. In P. Gianni, editor, Proceedings of ISSAC’88, volume 358 of LNCS, pages 467–474. Springer Verlag, 1989.

[16] M. Kaminski, D. Kirkpatrick, and N. Bshouty. Addition requirements for matrix and transposed matrix products. Journal of Algorithms, 9(3):354–364, 1988.

[17] G. Lecerf and É. Schost. Fast multivariate power series multiplication in characteristic zero. SADIO Electronic Journal on Informatics and Operations Research. To appear, manuscript of December 2002.

[18] P. L. Montgomery. An FFT extension of the elliptic curve method of factorization. PhD thesis, University of California, Los Angeles CA, 1992.

[19] P. Penfield, Jr., R. Spencer, and S. Duinker. Tellegen’s theorem and electrical networks. The M.I.T. Press, Cambridge, Mass.-London, 1970.

[20] V. Shoup. NTL: A library for doing number theory.

http://www.shoup.net.

[21] V. Shoup. A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic. In Proceedings of ISSAC’91, pages 14–21. ACM, 1991.

[22] V. Shoup. A new polynomial factorization algorithm and its implementation. Journal of Symbolic Computation, 20(4):363–397, 1995.

[23] V. Shoup. Efficient computation of minimal polynomials in algebraic extensions of finite fields. In Proceedings of ISSAC’99, pages 53–58. ACM, 1999.

[24] B. Tellegen. A general network theorem, with applications. Philips Research Reports, 7:259–269, 1952.

[25] R. Zippel. Interpolating polynomials from their values. Journal of Symbolic Computation, 9(3):375–403, 1990.