Fast solution of toeplitz systems of equations and computation of Padé approximants

We present two new algorithms, ADT and MDT, for solving order-n Toeplitz systems of linear equations Tz = b in time O(n log²n) and space O(n). The fastest algorithms previously known, such as Trench's algorithm, require time $\text{Ω(n}{}^2\text{)}$ and require that all principal submatrices of T be nonsingular. Our algorithm ADT requires only that T be nonsingular. Both our algorithms for Toeplitz systems are derived from algorithms for computing entries in the Padé table for a given power series. We prove that entries in the Padé table can be computed by the Extended Euclidean Algorithm. We describe an algorithm EMGCD (Extended Middle Greatest Common Divisor) which is faster than the algorithm HGCD of Aho, Hopcroft and Ullman, although both require time O(n log²n), and we generalize EMGCD to produce PRSDC (Polynomial Remainder Sequence Divide and Conquer) which produces any iterate in the PRS, not just the middle term, in time O(n log²n). Applying PRSDC to the polynomials U₀(x) = x²ⁿ⁺¹ and U₁(x) = a₀ + a₁x + … + a_2nx²ⁿ gives algorithm AD (Anti-Diagonal), which computes any (m, p) entry along the antidiagonal m + p = 2n of the Padé table for U₁ in time O(n log²n). Our other algorithm, MD (Main-Diagonal), computes any diagonal entry (n, n) in the Padé table for a normal power series, also in time O(n log²n). MD is related to Schönhage's fast continued fraction algorithm. A Toeplitz matrix T is naturally associated with U₁, and the (n, n) Padé approximation to U₁ gives the first column of T^?1. We show how a formula due to Trench can be used to compute the solution z of Tz = b in time O(n log n) from the first row and column of T^?1. Thus, the Padé table algorithms AD and MD give O(n log²n) Toeplitz algorithms ADT and MDT. Trench's formula breaks down in certain degenerate cases, but in such cases a companion formula, the discrete analog of the Christoffel-Darboux formula, is valid and may be used to compute z in time O(n log²n) via the fast computation (by algorithm AD) of at most four Padé approximants. We also apply our results to obtain new complexity bounds for the solution of banded Toeplitz systems and for BCH decoding via Berlekamp's algorithm.     

Efficient partition trees

We prove a theorem on partitioning point sets inE ^d (d fixed) and give an efficient construction of partition trees based on it. This yields a simplex range searching structure with linear space,O(n logn) deterministic preprocessing time, andO(n ^1?1/d (logn) ^O(1)) query time. WithO(nlogn) preprocessing time, where δ is an arbitrary positive constant, a more complicated data structure yields query timeO(n ^1?1/d (log logn) ^O(1)). This attains the lower bounds due to Chazelle [C1] up to polylogarithmic factors, improving and simplifying previous results of Chazelleet al. [CSW]. The partition result implies that, forr ^d≤n ^1?δ, a (1/r)-approximation of sizeO(r ^d) with respect to simplices for ann-point set inE ^d can be computed inO(n logr) deterministic time. A (1/r)-cutting of sizeO(r ^d) for a collection ofn hyperplanes inE ^d can be computed inO(n logr) deterministic time, provided thatr≤n ^1/(2d?1).     

Essentially optimal computation of the inverse of generic polynomial matrices

We present an inversion algorithm for nonsingular n×n matrices whose entries are degree d polynomials over a field. The algorithm is deterministic and, when n is a power of two, requires O^∼(n³d) field operations for a generic input; the soft-O notation O^∼ indicates some missing log(nd) factors. Up to such logarithmic factors, this asymptotic complexity is of the same order as the number of distinct field elements necessary to represent the inverse matrix.     

Optimal multiway split trees

A class of multiway split trees is defined. Given a set of n weighted keys and a node capacity m, an algorithm is described for constructing a multiway split tree with minimum access cost. The algorithm runs in time O(mn⁴) and requires O(mn³) storage locations. A further refinement of the algorithm enables the factor m in the above costs to be reduced to log m.     

Localized and compact data-structure for comparability graphs

We show that every comparability graph of any two-dimensional poset over n elements (a.k.a. permutation graph) can be preprocessed in O(n) time, if two linear extensions of the poset are given, to produce an O(n) space data-structure supporting distance queries in constant time. The data-structure is localized and given as a distance labeling, that is each vertex receives a label of O(logn) bits so that distance queries between any two vertices are answered by inspecting their labels only. This result improves the previous scheme due to Katz, Katz and Peleg [M. Katz, N.A. Katz, D. Peleg, Distance labeling schemes for well-separated graph classes, Discrete Applied Mathematics 145 (2005) 384-402] by a log n factor.As a byproduct, our data-structure supports all-pair shortest-path queries in O(d) time for distance-d pairs, and so identifies in constant time the first edge along a shortest path between any source and destination.More fundamentally, we show that this optimal space and time data-structure cannot be extended for higher dimension posets. More precisely, we prove that for comparability graphs of three-dimensional posets, every distance labeling scheme requires Ω(n^1/3) bit labels.     

Tensor-Structured Preconditioners and Approximate Inverse of Elliptic Operators in ℝ d

In the present paper we analyze a class of tensor-structured preconditioners for the multidimensional second-order elliptic operators in ? ^d , d≥2. For equations in a bounded domain, the construction is based on the rank-R tensor-product approximation of the elliptic resolvent ? _R ≈(??λ I)^?1, where ? is the sum of univariate elliptic operators. We prove the explicit estimate on the tensor rank R that ensures the spectral equivalence. For equations in an unbounded domain, one can utilize the tensor-structured approximation of Green’s kernel for the shifted Laplacian in ? ^d , which is well developed in the case of nonoscillatory potentials. For the oscillating kernels e ^{?i κ‖x‖}/‖x‖, x∈? ^d , κ∈?₊, we give constructive proof of the rank-O(κ) separable approximation. This leads to the tensor representation for the discretized 3D Helmholtz kernel on an n×n×n grid that requires only O(κ?|log?ε|²? n) reals for storage. Such representations can be applied to both the 3D volume and boundary calculations with sublinear cost O(n ²), even in the case κ=O(n). Numerical illustrations demonstrate the efficiency of low tensor-rank approximation for Green’s kernels e ^?λ‖x‖/‖x‖, x∈?³, in the case of Newton (λ=0), Yukawa (λ∈?₊), and Helmholtz (λ=i κ,?κ∈?₊) potentials, as well as for the kernel functions 1/‖x‖ and 1/‖x‖ ^d?2, x∈? ^d , in higher dimensions d>3. We present numerical results on the iterative calculation of the minimal eigenvalue for the d-dimensional finite difference Laplacian by the power method with the rank truncation and based on the approximate inverse ? _R ≈(?Δ)^?1, with 3≤d≤50.     

Factoring matrices with a tree-structured sparsity pattern

Let A be a matrix whose sparsity pattern is a tree with maximal degree d_max. We show that if the columns of A are ordered using minimum degree on |A|+|A|^∗, then factoring A using a sparse LU with partial pivoting algorithm generates only O(d_maxn) fill, requires only O(d_maxn) operations, and is much more stable than LU with partial pivoting on a general matrix. We also propose an even more efficient and just-as-stable algorithm called sibling-dominant pivoting. This algorithm is a strict partial pivoting algorithm that modifies the column preordering locally to minimize fill and work. It leads to only O(n) work and fill. More conventional column pre-ordering methods that are based (usually implicitly) on the sparsity pattern of |A|^∗|A| are not as efficient as the approaches that we propose in this paper.     

A Baby Steps/Giant Steps Probabilistic Algorithm for Computing Roadmaps in Smooth Bounded Real Hypersurface

We consider the problem of constructing roadmaps of real algebraic sets. This problem was introduced by Canny to answer connectivity questions and solve motion planning problems. Given s polynomial equations with rational coefficients, of degree D in n variables, Canny’s algorithm has a Monte Carlo cost of sⁿlog(s)D^O(n2)s^{n}\log(s)D^{O(n^{2})} operations in ℚ; a deterministic version runs in time sⁿlog(s)D^O(n4)s^{n}\log(s)D^{O(n^{4})} . A subsequent improvement was due to Basu, Pollack, and Roy, with an algorithm of deterministic cost s^d+1D^O(n2)s^{d+1}D^{O(n^{2})} for the more general problem of computing roadmaps of a semi-algebraic set (d≤n is the dimension of an associated object).     

Dynamic point location in arrangements of hyperplanes

We present algorithms for maintaining data structures supporting fast (polylogarithmic) point-location and ray-shooting queries in arrangements of hyperplanes. This data structure allows for deletion and insertion of hyperplanes. Our algorithms use random bits in the construction of the data structure but do not make any assumptions about the update sequence or the hyperplanes in the input. The query bound for our data structure isÕ(polylog(n)), wheren is the number of hyperplanes at any given time, and theÕ notation indicates that the bound holds with high probability, where the probability is solely with respect to randomization in the data structure. By high probability we mean that the probability of error is inversely proportional to a large degree polynomial inn. The space requirement isÕ(n ^d). The cost of update isÕ(n ^d?1 logn. The expected cost of update isO(n ^d?1); the expectation is again solely with respect to randomization in the data structure. Our algorithm is extremely simple. We also give a related algorithm with optimalÕ(logn) query time, expectedO(n ^d) space requirement, and amortizedO(n ^d?1) expected cost of update. Moreover, our approach has a versatile quality which is likely to have further applications to other dynamic algorithms. Ford=2, 3 we also show how to obtain polylogarithmic update time in the CRCW PRAM model so that the processor-time product matches (within a polylogarithmic factor) the sequential update time.     

Parallel solution of certain Toeplitz least-squares problems

We describe a systolic algorithm for solving a Toeplitz least-squares problem of special form. Such problems arise, for example, when Volterra convolution equations of the first kind are solved by regularization. The systolic algorithm is based on a sequential algorithm of Eldén, but we show how the storage requirements of Eldén's algorithm can be reduced from O(n²) to O(n). The sequential algorithm takes time O(n²); the systolic algorithm takes time O(n) using a linear systolic array of O(n) cells. We also show how large problems may be decomposed and solved on a small systolic array.     

Range Searching and Point Location among Fat Objects

We present a data structure that can store a set of disjoint fat objects ind-space such that point location and bounded-size range searching with arbitrarily shaped ranges can be performed efficiently. The structure can deal with either arbitrary (fat) convex objects or nonconvex (fat) polytopes. The multipurpose data structure supports point location and range searching queries in timeO(log^d−1 n) and requiresO(n log^d−1 n) storage, afterO(n log^d−1 n log log n) preprocessing. The data structure and query algorithm are rather simple.     

On the complexity of calculating factorials

It is shown that n! can be evaluated with time complexity O(log log nM (n log n)), where M(n) is the complexity of multiplying two n-digit numbers together. This is effected, in part, by writing n! in terms of its prime factors. In conjunction with a fast multiplication this yields an O(n(log n log log n)²) complexity algorithm for n!. This might be compared to computing n! by multiplying 1 times 2 times 3, etc., which is ω(n² log n) and also to computing n! by binary splitting which is O(log nM(n log n)).     

Fast construction of irreducible polynomials over finite fields

We present a randomized algorithm that on inputting a finite field K with q elements and a positive integer d outputs a degree d irreducible polynomial in K[x]. The running time is d ^1+?(d)×(log q)^5+?(q) elementary operations. The function ? in this expression is a real positive function belonging to the class o(1), especially, the complexity is quasi-linear in the degree d. Once given such an irreducible polynomial of degree d, we can compute random irreducible polynomials of degree d at the expense of d ^1+?(d) × (log q)^1+?(q) elementary operations only.     

A Parallel Reduction of Hamiltonian Cycle to Hamiltonian Path in Tournaments

We propose a parallel algorithm which reduces the problem of computing Hamiltonian cycles in tournaments to the problem of computing Hamiltonian paths. The running time of our algorithm is O(log n) using O(n²/log n) processors on a CRCW PRAM, and O(log n log log n) on an EREW PRAM using O(n²/log n log log n) processors. As a corollary, we obtain a new parallel algorithm for computing Hamiltonian cycles in tournaments. This algorithm can be implemented in time O(log n) using O(n²/log n) processors in the CRCW model and in time O(log²n) with O(n²/log n log log n) processors in the EREW model.     

Sorting under partial information (without the ellipsoid algorithm)

We revisit the well-known problem of sorting under partial information: sort a finite set given the outcomes of comparisons between some pairs of elements. The input is a partially ordered set P, and solving the problem amounts to discovering an unknown linear extension of P, using pairwise comparisons. The information-theoretic lower bound on the number of comparisons needed in the worst case is log e(P), the binary logarithm of the number of linear extensions of P. In a breakthrough paper, Jeff Kahn and Jeong Han Kim (J. Comput. System Sci. 51 (3), 390–399, 1995) showed that there exists a polynomial-time algorithm for the problem achieving this bound up to a constant factor. Their algorithm invokes the ellipsoid algorithm at each iteration for determining the next comparison, making it impractical. We develop efficient algorithms for sorting under partial information. Like Kahn and Kim, our approach relies on graph entropy. However, our algorithms differ in essential ways from theirs. Rather than resorting to convex programming for computing the entropy, we approximate the entropy, or make sure it is computed only once, in a restricted class of graphs, permitting the use of a simpler algorithm. Specifically, we present:

1. an O(n ²) algorithm performing O(logn·log e(P)) comparisons
2. an O(n ^2:5) algorithm performing at most (1+?) log e(P)+O _? (n) comparisons
3. an O(n ^2:5) algorithm performing O(log e(P)) comparisons.

All our algorithms can be implemented in such a way that their computational bottleneck is confined in a preprocessing phase, while the sorting phase is completed in O(q)+O(n) time, where q denotes the number of comparisons performed.     

Binary Segmentation for Multivariate Polynomials

This paper develops a fast method of binary segmentation for multivariate integer polynomials and applies it to polynomial multiplication, pseudodivision, resultant and gcd calculations. In the univariate case, binary segmentation yields the fastest known method for multiplication of integer polynomials, even slightly faster than fast Fourier transform methods (however, the underlying fast integer multiplication method is based on fast Fourier transforms). Fischer and Paterson ("SIAM-AMS Proceedings, 1974," pp. 113-125) seem to be the first authors who suggest binary segmentation for polynomials with coefficients in {0, 1}. Schönhage (J. Complexity1 (1985), 118-137), as well as Bini and Pan (J. Complexity2 (1986), 179-203) have applied the binary segmentation to univariate polynomial division and gcd calculation. In the multivariate case, well- known methods are the iterated application of univariate binary segmentation and Kronecker′s map. Our method of binary segmentation to the contrary is based on a one-step algebra homomorphism mapping multivariate polynomials directly into integers. More precisely, the variables x₁, . . . , x_n of a multivariate polynomial are substituted by integers 2^ν₁, . . . , 2^ν_n, where ν₁ < · · · < ν_n are chosen as small as possible. If n is the number of variables, d bounds the total degrees, and l bounds the bit lengths of the input, we achieve the bit complexities O(ψ((2d)ⁿ (n log d + l))) for multivariate multiplication, O(d ψ(d²ⁿ(n log d + l))) for multivariate pseudo-division (provided n = O(d log d)), and O(d ψ((2d²)ⁿ (n log d + l))) for multivariate subresultant chain calculation. Here ψ(l) = l log l log log l denotes the complexity of integer multiplication.     

Difference sets and shifted primes

We show that if A is a subset of {1, …, n} which has no pair of elements whose difference is equal to p ? 1 with p a prime number, then the size of A is O(n(log log n)^{?c(log log log log log n)}) for some absolute c > 0.     

Fast enumeration of words generated by Dyck grammars

The problem of enumerating and denumerating words generated by Dyck grammars arises in the work of compilers for high-level programming languages and a number of other applications. The present paper proposes an algorithm for the fast enumeration and denumeration of words of Dyck languages; the complexity of this algorithm per one symbol of enumerated words is O(log³ n log log n) bit operations, provided that the Schönhage-Strassen multiplication and division algorithmis used. The well-knownmethods applied earlier possess complexityO(n) per one symbol of enumerated words. The construction of the proposed algorithm is based on the Ryabko method.     

Systolic implementation of the lattice algorithm for least squares linear prediction problems

We present a linear systolic array of O(n) cells that solves the least squares linear prediction problem in time O(n) via an algorithm based on the so-called lattice algorithm. The total storage required is O(n) words, i.e., only a constant number of words are needed at each cell.     

An O(n log n) unidirectional distributed algorithm for extrema finding in a circle

In this paper we present algorithms, which given a circular arrangement of n uniquely numbered processes, determine the maximum number in a distributive manner. We begin with a simple unidirectional algorithm, in which the number of messages passed is bounded by 2 n log n + O(n). By making several improvements to the simple algorithm, we obtain a unidirectional algorithm in which the number of messages passed is bounded by 1.5nlogn + O(n). These algorithms disprove Hirschberg and Sinclair's conjecture that O(n²) is a lower bound on the number of messages passed in undirectional algorithms for this problem. At the end of the paper we indicate how our methods can be used to improve an algorithm due to Peterson, to obtain a unidirectional algorithm using at most 1.356nlogn + O(n) messages. This is the best bound so far on the number of messages passed in both the bidirectional and unidirectional cases.