Multiplying polynomials with ball coefficients

Let $A = \sum_i a_i x^i$ and $B = \sum_j b_j x^j$ be polynomials whose coefficients are balls $a_i = [m_i \pm r_i]$ , $b_j = [m'_j \pm r'_j]$ . A ball product $C$ must contain $\sum_{i+j=k} \alpha_i \beta_j$ for all $\alpha_i \in a_i$ , $\beta_j \in b_j$ and every $k$ . The classical algorithm evaluates each such sum as one ball dot product (Dot products of balls), which costs $O(n^2)$ multiplications for length $n$ but gives nearly optimal enclosures. Fast algorithms (Karatsuba, FFT) cannot simply be run in ball arithmetic: they rely on cancellations that ball arithmetic does not see, so the radii of the output grow like the sum of the absolute values of all intermediate quantities, which for polynomials with coefficients of very different magnitudes destroys the accuracy entirely ([vdH2008], [Hig2002]). FLINT uses the following algorithm instead: the midpoints are multiplied exactly with fast integer polynomial arithmetic, after a scaling and a splitting into blocks that keep the integers small, and the radii are bounded by a separate product of nonnegative polynomials of low precision, following van der Hoeven [vdH2008] and Johansson [Joh2017].

Error bound for the product#

Write $A = A_m + E$ with $A_m = \sum m_i x^i$ and $E = \sum \varepsilon_i x^i$ , $|\varepsilon_i| \le r_i$ , and likewise $B = B_m + E'$ . For a polynomial $P$ , let $|P|$ denote the polynomial whose coefficients are the absolute values of those of $P$ , and write $P \preceq Q$ if every coefficient of $P$ is at most the corresponding coefficient of $Q$ .

Lemma 1. With the notation above, $A_r = \sum r_i x^i$ and $B_r = \sum r'_j x^j$ ,

|AB - A_m B_m| \preceq |A_m| B_r + A_r (|B_m| + B_r), \qquad |A^2 - A_m^2| \preceq A_r (2|A_m| + A_r).

Proof. $AB - A_m B_m = A_m E' + E B_m + E E' = A_m E' + E(B_m + E')$ , and the coefficients of a product of polynomials satisfy $|PQ| \preceq |P|\,|Q|$ by the triangle inequality, with $|E| \preceq A_r$ , $|E'| \preceq B_r$ and $|B_m + E'| \preceq |B_m| + B_r$ . For the square, $A^2 - A_m^2 = 2 A_m E + E^2 = E(2A_m + E)$ .

So it suffices to compute the product $A_m B_m$ of the midpoints with rigorously bounded rounding error, and upper bounds for one or two products of polynomials with nonnegative coefficients. These products need only a few bits of relative accuracy, because they bound errors.

Scaling and blocks#

The coefficients of a power series often decay or grow geometrically: the coefficients of $e^x$ are $1/k!$ , those of $1/(1 - 100x)$ are $100^k$ . Representing all of them as integers with a common exponent would need $O(n \log n)$ or $O(n)$ bits per coefficient. The substitution $x \to 2^{-c} x$ , which multiplies the coefficient of $x^k$ by $2^{-ck}$ , is exact in binary floating-point arithmetic and commutes with multiplication: if $\tilde A(x) = A(2^{-c} x)$ , then $\tilde A(x) \tilde B(x) = (AB)(2^{-c} x)$ , and the coefficients of $AB$ are those of $\tilde A \tilde B$ multiplied by $2^{ck}$ .

FLINT chooses $c$ as the nearest integer to the average slope of the exponents between the first and the last nonzero midpoint of the two operands: if the first and last nonzero midpoints of $A$ have indices $a_1 < a_2$ and exponents $e_1, e_2$ , and those of $B$ indices $b_1 < b_2$ and exponents $f_1, f_2$ , then

c = \left\lfloor \frac{(e_2 - e_1) + (f_2 - f_1)}{(a_2 - a_1) + (b_2 - b_1)} + \frac12 \right\rfloor.

This is cheap and gives slowly varying scaled coefficients for power series with similar decay rates; it is not optimal in general.

Next, the scaled midpoints of each operand are split into maximal runs of consecutive indices, called blocks, such that within a block all nonzero midpoints are integer multiples of a common power of two $2^{e}$ and have absolute value below $2^{e + H}$ , where the height $H$ is at most $3p + 512$ for the working precision $p$ . A block $t$ starting at index $s_t$ thus defines an integer polynomial $P_t$ with $\tilde A_m = \sum_t 2^{e_t} x^{s_t} P_t(x)$ , exactly, and similarly $\tilde B_m = \sum_u 2^{e'_u} x^{s'_u} Q_u(x)$ . Then

\tilde A_m \tilde B_m = \sum_{t, u} 2^{e_t + e'_u} x^{s_t + s'_u} P_t Q_u,

where each $P_t Q_u$ is a product of integer polynomials, computed exactly by _fmpz_poly_mulmid (Kronecker substitution, Schönhage–Strassen or the small-prime FFT, see Multiplying polynomials over Z). Each output coefficient is the sum of the contributions of all block pairs; each contribution is rounded to $p$ bits when it is added, with the rounding error added to the radius. Well-scaled operands form a single block each, and then the midpoint product is one exact integer polynomial product followed by one rounding per coefficient.

The products of nonnegative polynomials in Lemma 1 are computed in the same way. Their coefficients are mag_t numbers with 30-bit mantissas (Representation of radii and midpoints); after the same scaling they are split into blocks of height at most $3 \cdot 30 + 512 = 602$ bits. For a pair of blocks of which the shorter has fewer than 1000 coefficients, the product is computed in double-precision floating-point arithmetic, each coefficient being converted exactly to a double with a fixed exponent offset (the block height guarantees that neither the operands nor the products overflow or underflow); otherwise it is computed exactly with integers. The results are converted to upper bounds with mag_set_d_2exp_fmpz and mag_set_fmpz_2exp_fmpz, which round upwards.

Lemma 2 (Floating-point sums of nonnegative products). Let $s = \sum_{i=1}^{L} x_i y_i$ with nonnegative doubles $x_i, y_i$ , $L \le 1000$ , and let $\tilde s$ be the value computed by multiplying and adding in any order in IEEE double precision with rounding to nearest, without overflow or underflow. Then $s \le (1 + 10^{-9}) \tilde s$ .

Proof. Each product and each partial sum is rounded once, with a relative error factor $1 + \delta$ , $|\delta| \le u = 2^{-53}$ . Since all quantities are nonnegative, each term $x_i y_i$ enters $\tilde s$ multiplied by a product of at most $L$ such factors, each at least $1 - u$ , so $\tilde s \ge (1 - u)^L s$ . Hence $s \le (1 - u)^{-L} \tilde s \le (1 + 1.2 \cdot 10^{-13}) \tilde s$ for $L \le 1000$ . (A fused multiply-add only removes roundings.)

The implementation multiplies the double results by $1 + 10^{-9}$ , a comfortable margin over Lemma 2.

Algorithm 1: Block multiplication

Input: $A$ , $B$ with ball coefficients, a window $[\ell, h)$ of output indices, precision $p$ . Output: balls $c_\ell, \ldots, c_{h-1}$ .

Drop the coefficients of each operand that only contribute below $x^\ell$ or above $x^{h-1}$ , and the trailing zero coefficients.
If some coefficient is not finite, use the classical algorithm.
Compute the scaling exponent $c$ .
Radii: compute upper bounds for the coefficients $\ell, \ldots, h-1$ of $|\tilde A_m| \tilde B_r + \tilde A_r(|\tilde B_m| + \tilde B_r)$ (or of $\tilde A_r(2|\tilde A_m| + \tilde A_r)$ for a square), omitting the terms with $A_r = 0$ or $B_r = 0$ , by block products as described above, and store them as the radii of $c_k$ .
Midpoints: split $\tilde A_m$ , $\tilde B_m$ into blocks of height at most $3p + 512$ ; for each pair of blocks, compute the exact middle product of the integer polynomials and add its coefficients, times the power of two, to the balls $c_k$ with rounding to $p$ bits.
Multiply $c_k$ by $2^{ck}$ .

Theorem 3. For every choice of $\alpha_i \in a_i$ and $\beta_j \in b_j$ , the coefficient $\sum_{i+j=k} \alpha_i \beta_j$ of the product lies in the ball $c_k$ computed by Algorithm 1, for $\ell \le k < h$ .

Proof. Step 1 only removes terms that do not occur in the requested coefficients. The scaling in steps 3 and 6 is exact. In step 5, the block representation of the midpoints is exact and the block products are exact, so the only errors are the roundings when adding the contributions, which ball addition adds to the radii. The midpoint of the result after step 5 therefore differs from the coefficient of $\tilde A_m \tilde B_m$ by at most the radius accumulated in step 5. By Lemma 1, the coefficient of the exact product differs from that of $\tilde A_m \tilde B_m$ by at most the coefficient of the bound polynomial, for which step 4 computed an upper bound: the conversions of midpoints to mag_t in step 4 round upwards, the integer products are exact, the double products are bounded by Lemma 2, and the final conversions round upwards.

The block algorithm is about as accurate as the classical one when the scaled coefficients of each operand vary by much less than $2^{3p}$ , as is typical for power series of analytic functions: each coefficient of the product is then rounded once, to $p$ bits relative to the largest term in its sum. For badly scaled operands it degenerates gracefully to many small blocks, and in the worst case (every coefficient a block of its own) to the classical algorithm with more overhead.

Selection and cost#

_arb_poly_mullow uses the classical algorithm if the output or one of the operands has length at most 7, or if the output length or twice the length of the shorter operand is at most a cutoff that is 110 for $p \le 128$ and $10^4/(\ln p)^3$ clamped to $[8, 60]$ above (increased by 25% for squares at $p \ge 256$ and again at $p \ge 4096$ ). At high precision, the integer products are long, and the block algorithm wins earlier.

For operands of length $n$ that form single blocks, the midpoint product is one integer polynomial product with $O(p + \log n)$ -bit coefficients, that is, $O(M(n(p + \log n)))$ bit operations, and the radius products cost $O(n^2)$ double operations for $n < 1000$ (fast in practice, since they are simple dot products) and $O(M(n \log n))$ beyond. The classical algorithm costs $O(n^2 M(p))$ .

Complex polynomials#

For complex polynomials, _acb_poly_mullow and _acb_poly_mulmid first try _acb_poly_mulmid_fmpzi, when _acb_poly_mulmid_use_fmpzi predicts that it is faster (roughly, when the number of coefficient products per input coefficient times $p^{0.57}$ is at least 240, or 330 for squares). It scales both operands as above, using the larger of the real and imaginary midpoints for the slope, writes all real and imaginary midpoints of each operand as Gaussian integers times a single power of two, and fails if this needs more than $3p + 512$ bits. The midpoints of the product are then one product of Gaussian integer polynomials (_fmpzi_poly_mulmid), and the radii of the real and imaginary parts are bounded separately: for $(a_r + i a_i)(b_r + i b_i)$ with radii $\rho$ ,

\operatorname{rad}(\operatorname{Re}) \le |a_r| \rho_{b_r} + \rho_{a_r}(|b_r| + \rho_{b_r}) + |a_i| \rho_{b_i} + \rho_{a_i}(|b_i| + \rho_{b_i})

and symmetrically for the imaginary part, which is Lemma 1 applied to $a_r b_r - a_i b_i$ and $a_r b_i + a_i b_r$ ; the eight (four for a square) products of nonnegative polynomials are computed by _arb_poly_addmulmid_mag. If the operands are not well-scaled, the transposed algorithm computes the real and imaginary parts of the product with four real products _arb_poly_mulmid, each of which uses its own blocks. The Gauss variant with three real products, $(a + bi)(c + di) = (ac - bd) + ((a + b)(c + d) - ac - bd) i$ , is faster but loses accuracy when $|a|$ and $|b|$ (or $|c|$ and $|d|$ ) differ greatly, since the error of the large products then swamps the small part; it is available but not used by default.

Use in FLINT#

_arb_poly_mullow_block and _arb_poly_mulmid_block implement Algorithm 1; _arb_poly_mullow_classical and _arb_poly_mulmid_classical the classical algorithm with arb_dot; _arb_poly_mullow, _arb_poly_mulmid and _arb_poly_mul choose between them. _arb_poly_addmulmid_mag exposes the radius products.
_acb_poly_mulmid_fmpzi, _acb_poly_mullow_transpose, _acb_poly_mullow_transpose_gauss and _acb_poly_mullow_classical are the complex algorithms, selected by _acb_poly_mullow.
Over the generic ball contexts of gr_poly.h, _gr_poly_mullow and _gr_poly_mulmid call these functions, so that every generic polynomial algorithm over arb and acb coefficients uses them.
Matrix multiplication uses the same idea of exact integer products of scaled blocks (Ball matrix multiplication).