In this paper, we consider the problem of deciding the existence of real solutions to a system of polynomial equations having real coefficients, and which are invariant under the action of the symmetric group. We construct and analyze a Monte Carlo probabilistic algorithm which solves this problem, under some regularity assumptions on the input, by taking advantage of the symmetry invariance property.

A Las Vegas randomized algorithm is given to compute the Hermite normal form of a nonsingular integer matrix A of dimension n . The algorithm uses quadratic integer multiplication and cubic matrix multiplication and has running time bounded by O(n 3 (log n + log ||A||) 2 (log n ) 2 ) bit operations, where || A ||= max ij | A ij | denotes the largest entry of A in absolute value. A variant of the algorithm that uses pseudo-linear integer multiplication is given that has running time (n 3 log || A ||) 1+ o (1) bit operations, where the exponent “ + o (1)” captures additional factors c 1 (log n ) c2 (loglog|| A ||) c3 for positive real constants c 1 ,c 2 ,c 3 .

A Bohemian matrix family is a set of matrices all of whose entries are drawn from a fixed, usually discrete and hence bounded, subset of a field of characteristic zero. Originally these were integers---hence the name, from the acronym BOunded HEight Matrix of Integers (BOHEMI)---but other kinds of entries are also interesting. Some kinds of questions about Bohemian matrices can be answered by numerical computation, but sometimes exact computation is better. In this paper we explore some Bohemian families (symmetric, upper Hessenberg, or Toeplitz) computationally, and answer some open questions posed about the distributions of eigenvalue densities.

Consider a matrix F ε K [x]^mxn of univariate polynomials over a field K. We study the problem of computing the column rank profile of F. To this end we first give an algorithm which improves the minimal kernel basis algorithm of Zhou, Labahn, and Storjohann (Proceedings ISSAC 2012). We then provide a second algorithm which computes the column rank profile of F with a rank-sensitive complexity of O~ (rw-2n(m+d)) operations in K. Here, D is the sum of row degrees of F, w is the exponent of matrix multiplication, and O~ (.) hides logarithmic factors.

We present a Las Vegas randomized algorithm to compute the Smith normal form of a nonsingular integer matrix. For an A ∈ Zn×n, the algorithm requires O(n3(log n + log ||A||)2 (log n)2) bit operations using standard integer and matrix arithmetic, where ||A|| = maxij |Aij | denotes the largest entry in absolute value. Fast integer and matrix multiplication can also be used, establishing that the Smith form can be computed in about the same number of bit operations as required to multiply two matrices of the same dimension and size of entries as the input matrix.

We present a deterministic reduction to matrix multiplication for the problem of linear system solving: given as input a nonsingular A \in \Z^n \times n and b \in \Z^n \times 1 , compute A^-1 b. We give an algorithm that computes the minimal integer e such that all denominators of the entries in 2^eA^-1 are relatively prime to 2. Then, for a b that has entries with bitlength O(n) times as large as the bitlength of entries in A, we give an algorithm to produce the 2-adic expansion of 2^eA^-1 b up to a precision high enough such that A^-1 b over \Q can be recovered using rational number reconstruction. Both e and the 2-adic expansion can be computed in O(\MM(n,łog n + łog ||A||) \times (łog n) (łog n + łoglog ||A||)) bit operations. Here, ||A||= \max_ij |A_ij | and \MM(n,d) is the cost to multiply together, modulo 2^d, two n \times n integer matrices. Our approach is based on the previously known reductions of linear system solving to matrix multiplication which use randomization to find an integer lifting modulus X that is relatively prime to \det A. Here, we derandomize by first computing a permutation P, a unit upper triangular M, and a diagonal S with \det S a power of two, such that U := APMS^-1 is an integer matrix with 2 \perp \det U. This allows our modulus X to be chosen a power of 2.

We present a new algorithm for the computation of the integer-linear decomposition of a multivariate polynomial. Such a decomposition is used in Ore-Sato theory and discrete creative telescoping, for example to detect applicability of Zeilberger's algorithm to a hypergeometric term. Our algorithm is quite straightforward, requiring only basic polynomial arithmetic along with efficient rational root finding. We present complete complexity analyses for both our and previous algorithms in the case of bivariate integer polynomials, and show that our method has a better theoretical performance. We also provide a Maple implementation which shows that our method is faster in practice than previous algorithms.

We consider the problem of computing the nearest matrix polynomial with a non-trivial Smith Normal Form. We show that computing the Smith form of a matrix polynomial is amenable to numeric computation as an optimization problem. Furthermore, we describe an effective optimization technique to find a nearby matrix polynomial with a non-trivial Smith form. The results are later generalized to include the computation of a matrix polynomial having a maximum specified number of ones in the Smith Form (i.e., with a maximum specified McCoy rank). We discuss the geometry and existence of solutions and how our results can used for a backwards error analysis. We develop an optimization-based approach and demonstrate an iterative numerical method for computing a nearby matrix polynomial with the desired spectral properties. We also describe the implementation of our algorithms and demonstrate the robustness with examples in Maple.

Matrix polynomials appear in many areas of computational algebra, control systems theory, differential equations, and mechanics, typically with real or complex coefficients. Because of numerical error and instability, a matrix polynomial may appear of considerably higher rank (generically full rank), while being very close to a rank-deficient matrix. "Close" is defined naturally under the Frobenius norm on the underlying coefficient matrices of the matrix polynomial. In this paper we consider the problem of finding the nearest rank-deficient matrix polynomial to an input matrix polynomial, that is, the nearest square matrix polynomial which is algebraically singular. We prove that such singular matrices at minimal distance always exist (and we are never in the awkward situation having an infimum but no actual matrix polynomial at minimal distance). We also show that singular matrices at minimal distance are all isolated, and are surrounded by a basin of attraction of non-minimal solutions. Finally, we present an iterative algorithm which, on given input sufficiently close to a rank-deficient matrix, produces that matrix. The algorithm is efficient and is proven to converge quadratically given a sufficiently good starting point. An implementation demonstrates the effectiveness and numerical robustness in practice.

In this paper, we solve the existence problem of telescopers for rational functions in three discrete variables. We reduce the problem to that of deciding the summability of bivariate rational functions, a problem which has recently been solved. This existence criteria is used, for example, for detecting the termination of Zeilberger's algorithm to the function classes studied in this paper.

Given a rectangular matrix F ∈ K[x]mxn with m < n of univariate polynomials over a field K. we give an efficient algorithm for computing a unimodular completion of F. Our algorithm is deterministic and computes such a completion, when it exists, with cost O~ (nωs) field operations from K. Here s is the average of the m largest column degrees of F and ω is the exponent on the cost of matrix multiplication. Here O~ is big-O but with log factors removed. If a unimodular completion does not exist for F, our algorithm computes a unimodular completion for a right cofactor of a column basis of F, or equivalently, computes a completion that preserves the generalized determinant.

Given a matrix of univariate polynomials over a field K, its columns generate a K[x]-module. We call any basis of this module a column basis of the given matrix. Matrix gcds and matrix normal forms are examples of such module bases. In this paper we present a deterministic algorithm for the computation of a column basis of an m x n input matrix with m ≤ n. If s is the average column degree of the input matrix, this algorithm computes a column basis with a cost of Õ(nmω-1s) field operations in K. Here the soft-O notation is Big-O with log factors removed while ω is the exponent of matrix multiplication. Note that the average column degree s is bounded by the commonly used matrix degree that is also the maximum column degree of the input matrix.

Scalings form a class of group actions that have both theoretical and practical importance. A scaling is accurately described by an integer matrix. Tools from linear algebra are exploited to compute a minimal generating set of rational invariants, trivial rewriting and rational sections for such a group action. The primary tools used are Hermite normal forms and their unimodular multipliers. With the same line of ideas, a complete solution to the scaling symmetry reduction of a polynomial system is also presented.

In this paper we present a deterministic algorithm for the computation of a minimal nullspace basis of an m x n input matrix of univariate polynomials over a field K with m ≤ n. This algorithm computes a minimal nullspace basis of a degree d input matrix with a cost of O~ (nω ⌈md/n⌉) field operations in K. Here the soft-O notation is Big-O with log factors removed while ω is the exponent of matrix multiplication. The same algorithm also works in the more general situation on computing a shifted minimal nullspace basis, with a given degree shift [equation] whose entries bound the corresponding column degrees of the input matrix. In this case if ρ is the sum of the m largest entries of s, then a s-minimal right nullspace basis can be computed with a cost of O~(nωρ/m) field operations.

Recognizing hand drawn mathematical matrices on tablet computers has proven to be a particularly challenging task. While individual expression recognition can be simplified by assuming the entire content is a single semantic construct, a single math expression, a matrix is composed of multiple expressions arranged in rows and columns. These expressions must first be segmented into matrix elements, and then each individual matrix element expression must be recognized. In this work, we show how a simple algorithm on in-air (i.e. non-inking) strokes can be used to analyze the drawing order of a matrix. Once the drawing order is recognized, we show how outlier analysis on in-air packets gives rapid, reliable segmentation of matrix elements.

In this paper we give a new, fast algorithm for solving the simultaneous Padé approximation problem. The algorithm is fraction-free and is suitable for computation in domains where growth of coefficients in intermediate computations are a central concern. The algorithm gives significant improvement on previous fraction-free methods, in particular when solved via the use of vector Hermite-Padé approximation using the FFFG order basis algorithm previously done by the authors. The improvements are both in terms of bit complexity and in reduced size of the intermediate quantities.

In this paper we give an efficient algorithm for computation of order basis of a matrix of power series. For a problem with an m x n input matrix over a field K, m ≤ n, and order σ, our algorithm uses O(MM(n, ⊂O~(nω⌈mσ/n⌉) field operations in B.K, where the soft-O notation O~ is Big O with log factors omitted and MM(n,d) denotes the cost of multiplying two polynomial matrices with dimension n and degree d. The algorithm extends earlier work of Storjohann, whose method can be used to find a subset of an order basis that is within a specified degree bound δ using O~(MM(n,δ)) field operations for δ≥⌈ mσ/n⌉.

Many on-line (interactive) mathematics recognition systems allow the creation of typeset equations, normally in LaTeX, but they do not support mathematical problem solving. In this paper, we present MathBrush, a system that allows users to draw math input using a pen-input device on a tablet computer, recognizes the math expression, and then supports mathematical transformation and problem solving using back-end Computer Algebra Systems (CAS). We describe the architecture of the MathBrush system, which includes modules that support symbol recognition, semantic analysis, the transfer of recognized expressions to back-end CAS, and interface techniques for interacting with CAS output. We also identify unique challenges associated with recognition for math problem solving, such as the need for deeper semantic analysis than is required by LATEX, and the need to deal with ambiguities in user input. Our experiences serve to inform researchers seeking to design interactive mathematics recognition systems geared toward mathematical problem solving.

In this paper, we examine the problem of computing the greatest common divisor (GCD) of univariate polynomials represented in different bases. When the polynomials are represented in Newton basis or a basis of orthogonal polynomials, we show that the well-known Sylvester matrix can be generalized. We give fraction-free and modular algorithms to directly compute the GCD in the alternate basis. These algorithms are suitable for computation in domains where growth of coefficients in intermediate computations are a central concern. In the cases of Newton basis and bases using certain orthogonal polynomials, we also show that the standard subresultant algorithm can be applied easily. If the degrees of the input polynomials is at most n and the degree of the GCD is at least n/2, our algorithms outperform the corresponding algorithms using the standard power basis.

We consider the problem of sparse interpolation of an approximate multivariate black-box polynomial in floating-point arithmetic. That is, both the inputs and outputs of the black-box polynomial have some error, and all numbers are represented in standard, fixed-precision, floating point arithmetic. By interpolating the black box evaluated at random primitive roots of unity, we give efficient and numerically robust solutions. We note the similarity between the exact Ben-Or/Tiwari sparse interpolation algorithm and the classical Prony's method for interpolating a sum of exponential functions, and exploit the generalized eigenvalue reformulation of Prony's method. We analyze the numerical stability of our algorithms and the sensitivity of the solutions, as well as the expected conditioning achieved through randomization. Finally, we demonstrate the effectiveness of our techniques in practice through numerical experiments and applications.

We consider the problem of finding closed form solutions of linear differential equations having coefficients which are elliptic functions. For second order equations we show how to solve such an ode in terms of doubly periodic functions of the second kind. The method depends on two procedures, the first using a second symmetric power of an ode along with a decision procedure for determining when such equations have elliptic function solutions while the second involves the computation of exponential solutions.

An orthogonal Ore ring is an abstraction of common properties of linear partial differential, shift and q-shift operators. Using orthogonal Ore rings, we present an algorithm for finding hyperexponential solutions of a system of linear differential, shift and q-shift operators, or any mixture thereof, whose solution space is finite-dimensional. The algorithm is applicable to factoring modules over an orthogonal Ore ring when the modules are also finite-dimensional vector spaces over the field of rational functions.

We present a new algorithm for row reduction of a matrix of skew polynomials. The algorithm can be used for finding full rank decompositions and other rank revealing transformations of matrices of skew polynomials. The algorithm is intended for computation in exact arithmetic domains where the growth of coefficients in intermediate computations is a central concern. This coefficient growth is controlled by using fraction-free methods. This allows us to obtain a polynomial-time algorithm: for an m x s matrix of input skew polynomials of degree N with coefficients whose lengths are bounded by K the algorithm has a worst case complexity of O(m5s4N4K2) bit operations.

We present a new algorithm for determining all factorizations of a polynomial f in the domain ZN[x], a non-unique factorization domain, given in terms of parameters. From the prime factorization of N, the problem is reduced to factorization in Zph[x] where p is a prime and k ⪈ 1. If pk does not divide the discriminant of f and one factorization is given, our algorithm determines all factorizations with complexity Ο(n3 M(k log p)) where n denotes the degree of the input polynomial and M(t) denotes the complexity of multiplication of two t-bit numbers. Our algorithm improves on the method of von zur Gathen and Hartlieb, which has complexity Ο(n7 k(klog p + log n2). The improvement is achieved by processing all factors at the same time instead of one at a time and by computing the kernels and determinants of matrices over Zpk in an efficient manner.

This work demonstrates that it is possible to obtain accurate values of callable bonds using a fully numerical approach, provided that the PDE is discretized appropriately. To facilitate comparisons with results reported by Buttler and Waldvogel (1996), we consider models with a single factor: the instantaneous risk free interest rate. We emphasize, however, that it is straightforward to extend the numerical methods described to cases where the Green's function cannot be determined analytically as well as to cases with time-dependent parameters (typically used to match current term structures of interest rates/interest rate volatilities), or multi-factor interest rate models.

In t,his paper we st,ucly the problen ~ of transforrniug, via in-vertible colu1tln opcrat.ious ~ it matrix polyioruial into a varicty Of.shiftcd forms. Esarnplcs of forms c:overed in out frmwa-ork include a colunm rctluccd form: il triangular fornlz R I%!rInite IlOrInd fOrll1 or it Popov IlorIllal fOrIll alollg wit,11 their shifted courltcrpart,s. I3y obt.aiuiug tlcgrvc bounds for uuiniodiilar niiill,iplicrs of shifted Popor fornis we are able t,o c11lbct1 tlic probleni of conqmtiug il normal forni into 0Iic of deterniiJIing a sliift.ctl forni Of R InininIal pOl~IlOIIliill hiISiS fur all XWXiiltWl IIliLtris polyioruial. Shifted niiuind polynomial lXPX?S cm be conlpu1.d via sigma bases [2! 31 iIlld ill POpOv forui Vii1 Mahler s~sbenls [il. Tl ic 1. d t,t, cr Iut:t.lIotl gives a fractiorl-frw algorithm for computing niatris riormal forms. Key words: Popov Form. Herr&c N~~r~n~d Fornl 1

tire present a fraction-free approach to the computation of matrix Pad& systems.The method relies on determining a modified Schur complement for the coefficient matrices of the linear systems of equations that are associated to matrix Pad& approximation problems.By using this modified Schur complement for these matrices we are able to obtain a fast hybrid fraction-free algorithm for their computation.The algorithm is genera] and requires no extra assumptions on its input.

When a computer algebra system has an assumption facility, it is possible to distinguish between integration problems with respect to a real variable, and those with respect to a complex variable. Here, a class of integration problems is defined in which the integrand consists of compositions of continuous functions and signum functions, and integration is with respect to a real variable. Algorithms are given for evaluating such integrals. 1 Introduction In recent years, `assume' or `declare' facilities have been implemented in most of the available computer algebra systems (CAS). As well, such facilities have been gaining wider acceptance within the user community. The presence of these facilities has altered the way CAS behave, and many established areas of symbolic computation need to be reconsidered. The topic of this paper is an example of the impact on one traditional field of computer algebra, namely, symbolic integration. Because the early versions of many present-day CAS cou...

In a recent paper, Cabay, Jones, and Labahn develop a fast, iterative, lookahead algorithm for numerically computing Pade ´-Hermite systems and simultaneous Pade ´systems along a diagonal of the associated Pade ´tables.Included in their work is a detailed error analysis showing that the algorithm is weakly stable.In this article, we describe a Fortran implementation, VECTOR_PADE, of this algorithm together with a number of numerical experiments.These experiments show that the theoretical error bounds obtained by Cabay, Jones, and Labahn reflect the general behavior of the actual error, but that in practice these bounds are large overestimates.

This paper presents a new algorithm for computing the Hermite normal form H of an A c Z "m of rank m together with a unimodular pre-multiplier matrix U such that UA = H.Our algorithm requires O-(m '-lnM(mlog[lAl\)) bit oper@ions to produce both H and U. Here, IIAII = max,j lAij [, M(t) bit operations are sufficient to multiply two (t] -bit integers, and 0 is the exponent for matrix multiplication over rings: two m x m matrices over a ring R can be multiplied in O(me ) ring operations from R, The previously fastest algorithm of Hafner & McCurley requires 0-(m2nM (m log I IAI [)) bit operations to produce H, but does not produce a unimodular matrix U which satisfies UA = H.Previous methods require on the order of 0-(n3M(m log I[Al 1)) bit operations to produce a U -our algorithm improves on this significantly in both a theoretical and practical sense.

We present a Las Vegas probabilistic algorithm for reducing the computation of Hermite normal forms of rectangular polynomial matrices.In particular, the problem of computing the Hermite normal form of a rectangular m x n matrix (with m > n) reduces to that of computing the Hermite normal form of a matrix of size (n + 1) x n having entries of similar coefficient size and degree.The main cost of the reduction is the same as the cost of fraction-free Gaussian elimination of an m x n polynomial matrix.As an application, the reduction allows for the efficient computation of one-sided GCD'S of two matrix polynomials along with the solution of the matrix diophantine equation associated to such a GCD.

A set of new formulae for the inverse of a block Hankel (or block Toeplitz) matrix is given. The formulae are expressed in terms of certain matrix Padé forms, which approximate a matrix power series associated with the block Hankel matrix. By using Frobenius-type identities between certain matrix Padé forms, the inversion formulae are shown to generalize the formulae of Gohberg–Heinig and, in the scalar case, the formulae of Gohberg–Semencul and Gohberg–Krupnik. The new formulae have the significant advantage of requiring only that the block Hankel matrix itself be nonsingular. The other formulae require, in addition, that certain submatrices be nonsingular. Since effective algorithms for computing the required matrix Padé forms are available, the formulae are practical. Indeed, some of the algorithms allow for the efficient calculation of the inverse not only of the given block Hankel matrix, but also of any nonsingular block principal minor.

We present a new fast algorithm for the calculation of a Padé-Hermite form for a vector of power series. When the vector of power series is normal, the algorithm is shown to calculate a Padé-Hermite form of type (n0, ··· , nk) in O(k·(n02 + ··· + nk2)) operations. This complexity is the same as that of other fast algorithms for computing Padé-Hermite approximants. However, unlike other algorithms, the new algorithm also succeeds in the non-normal case, usually with only a moderate increase in cost.

For matrix power series with coefficients over a field, the notion of a matrix power series remainder sequence and its corresponding cofactor sequence are introduced and developed. An algorithm for constructing these sequences is presented. It is shown that the cofactor sequence yields directly a sequence of Pade fractions for a matrix power series represented as a quotient $B(z)^{ - 1} A(z)$. When $B(z)^{ - 1} A(z)$ is normal, the complexity of the algorithm for computing a Pade fraction of type $(m,n)$ is $O(p^3 (m + n)^2 )$, where p is the order of the matrices $A(z)$ and $B(z)$. For a power series that are abnormal for a given $(m,n)$, Padé fractions may not exist. However, it is shown that a generalized notion of Padé fraction, the Padé form, which is introduced in this paper, does always exist and can be computed by the algorithm. In the abnormal case, the algorithm can reach a complexity of $O(p^3 (m + n)^3 )$, depending on the nature of the abnormalities. In the special case of a scalar power series, however, the algorithm complexity is $O(m + n)^2 $, even in the abnormal case.