Flattening Karatsuba's recursion tree into a single summation

TL;DR

This paper introduces an interleaved splitting scheme for Karatsuba's recursion tree, simplifying access to nodes and enabling the flattening of the entire tree into a single convolution formula with fewer multiplications, applicable to infinite power series.

Contribution

It presents a novel interleaved splitting method that flattens Karatsuba's recursion tree into a single summation, reducing elementary multiplications and extending applicability to infinite series.

Findings

Flattened recursion tree into a single formula.

Reduced number of elementary multiplications.

Applicable to infinite power series.

Abstract

The recursion tree resulting from Karatsuba's formula is built here by using an interleaved splitting scheme rather than the traditional left/right one. This allows an easier access to the nodes of the tree and $2n-1$ of them are initially flattened all at once into a single recursive formula. The whole tree is then flattened further into a convolution formula involving less elementary multiplications than the usual Cauchy product. Unlike the traditional splitting scheme, the interleaved approach may also be applied to infinite power series, and corresponding formulas are also given.

Equations

$$AB=\underset{\text{branch 0}}{\underbrace{\left(X+1\right)A_0{B}_0}}+\underset{\text{branch 1}}{\underbrace{X\left(X+1\right)A_1{B}_1}}-\underset{\text{branch 2}}{\underbrace{X\left(A_1-A_0\right)\left(B_1-B_0\right)}}$$

$$\begin{array}[]{l}a_{7}x^{7}+a_{6}x^{6}+a_{5}x^{5}+a_{4}x^{4}+a_{3}x^{3}+a_{2}x^{2}+a_{1}x+a_{0}\\[8.0pt] \quad=\left(a_{7}x^{6}+a_{5}x^{4}+a_{3}x^{2}+a_{1}\right)x+\left(a_{6}x^{6}+a_{4}x^{4}+a_{2}x^{2}+a_{0}\right)\\[8.0pt] \quad=\left(\left(a_{7}x^{4}+a_{3}\right)x^{2}+\left(a_{5}x^{4}+a_{1}\right)\right)x+\left(\left(a_{6}x^{4}+a_{2}\right)x^{2}+\left(a_{4}x^{4}+a_{0}\right)\right)\end{array}$$

$$\frac{1-x^{2^m}}{1-x}{a}_r{b}_r{x}^r\text{with}r=(d_m{d}_{m-1}\dots d_3{d}_2{d}_1{)}_2$$

$$(1+x)(1+x^2)(1+x^4)\cdots =1+x+x^2+x^3+x^4+x^5+\dots x^{2^m-1}\text{.}$$

$$(x^1{)}^{d_1}(x^2{)}^{d_2}(x^4{)}^{d_3}\dots (x^{2^{m-1}}{)}^{d_m}$$

$$\frac{1-x^n}{1-x}\left(A\left(x\right)\odot B\left(x\right)\right)$$

$$\frac{1-x^n}{1-x^{2^k}}{x}^{2^{k-1}+j}$$

$$x^{2^{k-1}+j}\left(\frac{P(x)}{x^{2^{k-1}+j}}\times \frac{Q(x)}{x^{2^{k-1}+j}}\right)=\frac{P(x)\times Q(x)}{x^{2^{k-1}+j}}$$

$$\begin{array}[]{l}-\displaystyle\sum_{k=1}^{\log_{2}\left(n\right)}\;\displaystyle\sum_{j=0}^{2^{k-1}-1}\;\displaystyle\frac{1-x^{2^{k-1}}}{\left(1-x\right)x^{2^{k-1}+j}}\\[18.0pt] \qquad\quad\displaystyle\left(\displaystyle\frac{1-x^{n}}{1-x^{2^{k}}}x^{2^{k-1}+j}\odot\left(1-x^{2^{k-1}}\right)A\left(x\right)\right)\\[18.0pt] \qquad\times\left(\displaystyle\frac{1-x^{n}}{1-x^{2^{k}}}x^{2^{k-1}+j}\odot\left(1-x^{2^{k-1}}\right)B\left(x\right)\right)\end{array}$$

$$\begin{array}[]{l}-\displaystyle\sum_{m=1}^{n-1}\displaystyle\frac{1-x^{2^{\lfloor\log_{2}(m)\rfloor}}}{\left(1-x\right)x^{m}}\\[18.0pt] \qquad\quad\displaystyle\left(\displaystyle\frac{\left(1-x^{n}\right)x^{m}}{1-x^{2^{\lfloor\log_{2}(m)\rfloor+1}}}\odot\left(1-x^{2^{\lfloor\log_{2}(m)\rfloor}}\right)A\left(x\right)\right)\\[18.0pt] \qquad\times\left(\displaystyle\frac{\left(1-x^{n}\right)x^{m}}{1-x^{2^{\lfloor\log_{2}(m)\rfloor+1}}}\odot\left(1-x^{2^{\lfloor\log_{2}(m)\rfloor}}\right)B\left(x\right)\right)\end{array}$$

$$\begin{array}[]{l}f(x)\ast g(x)=\displaystyle\frac{f\left(x\right)\odot g\left(x\right)}{1-x}-\displaystyle\sum_{m=1}^{\infty}\;\displaystyle\frac{1-x^{2^{\lfloor\log_{2}m\rfloor}}}{\left(1-x\right)x^{m}}\\[18.0pt] \qquad\qquad\qquad\qquad\qquad\begin{array}[]{ll}&\displaystyle\left(\displaystyle\frac{x^{m}}{1-x^{2^{\lfloor\log_{2}m\rfloor+1}}}\odot\left(1-x^{2^{\lfloor\log_{2}m\rfloor}}\right)f\left(x\right)\right)\\[18.0pt] \ast&\displaystyle\left(\displaystyle\frac{x^{m}}{1-x^{2^{\lfloor\log_{2}m\rfloor+1}}}\odot\left(1-x^{2^{\lfloor\log_{2}m\rfloor}}\right)g\left(x\right)\right)\end{array}\end{array}$$

$$S_{d}=\left\{m\,\Big{|}\,1\leqslant m\leqslant d,\left(\left(d-m\right)\,\textrm{mod}\,\,2^{\lfloor\log_{2}m\rfloor+1}\right)<2^{\lfloor\log_{2}m\rfloor}\right\}$$

$$|S_d|={\text{A268289}}_d=d-2^k\tau \left(\frac{d+1}{2^k}-1\right)$$

$$u_{n+2^s}=u_n+\left(n+1\right)\left(\lfloor {\log }_2(n)\rfloor -s+2\right)+2^s-2^{\lfloor {\log }_2(n)\rfloor +1}$$

$$\begin{array}[]{l}a_{7}x^{7}+a_{6}x^{6}+a_{5}x^{5}+a_{4}x^{4}+a_{3}x^{3}+a_{2}x^{2}+a_{1}x+a_{0}\\[8.0pt] \quad=\left(a_{7}x^{6}+a_{5}x^{4}+a_{3}x^{2}+a_{1}\right)x+\left(a_{6}x^{6}+a_{4}x^{4}+a_{2}x^{2}+a_{0}\right)\\[8.0pt] \quad=\left(\left(a_{7}x^{4}+a_{3}\right)x^{2}+\left(a_{5}x^{4}+a_{1}\right)\right)x+\left(\left(a_{6}x^{4}+a_{2}\right)x^{2}+\left(a_{4}x^{4}+a_{0}\right)\right)\end{array}$$

$$\left(a_7-a_3\right),\left(a_5-a_1\right),\left(a_6-a_2\right)\text{and}\left(a_4-a_0\right)$$

$$\left(1+x\right)\left(1+x^2\right)\left(1+x^4\right)\left(1+x^8\right)\dots$$

$$\begin{array}[]{r}A\times B\,=\!\!\!\!\!\!\!{\displaystyle\sum_{\begin{array}[]{c}\text{\footnotesize$f\in\mathbb{Z}[X],$}\\ \text{\footnotesize$f\,\big{|}\,\sum_{k=0}^{n-1}X^{k}$}\end{array}}}\!\!\!\!\!\!\displaystyle f\left(f\,\dot{f}^{\textrm{\tiny\hskip 1.0ptW}}\odot\dot{f}A\odot\dot{f}B\right)\end{array}$$

$$A\times B=\sum _{k=0}^{n-1}{\bar{f}}_k\left({\sigma }_k{\bar{f}}_k{X}^k\odot {\dot{f}}_{k}A\odot {\dot{f}}_{k}B\right)$$

$$a(x)\ast b(x)=\sum _{k=0}^{\infty }{f}_k(x)\left({\sigma }_k{f}_k(x)x^k\odot {\dot{f}}_k(x)a(x)\odot {\dot{f}}_k(x)b(x)\right)$$

$$\begin{array}[]{l}\displaystyle g(x)\ast h(x)\,=\,\sum_{m=0}^{\infty}x^{m}\sum_{k=0}^{m}\sigma_{k}\\[12.0pt] \displaystyle\qquad\qquad\qquad\sum_{j=k}^{m}\tau_{k}(m,j)\,\tau_{k}(j,k)\left(\sum_{t=j-k}^{j}\upsilon_{k}(j,t)\,\sigma_{t}a_{t}\right)\left(\sum_{t=j-k}^{j}\upsilon_{k}(j,t)\,\sigma_{t}b_{t}\right)\end{array}$$

$$\begin{array}[]{l}\displaystyle g(x)\ast h(x)\,=\,\sum_{m=0}^{\infty}x^{m}\sum_{k=0}^{m}\sigma_{k}\\[12.0pt] \displaystyle\qquad\qquad\qquad\qquad\sum_{j=k}^{m}T(k+m-j,k)\,T(j,k)\\[12.0pt] \displaystyle\qquad\qquad\qquad\qquad\qquad\left(\sum_{t=j-k}^{j}T(k,j-t)\,\sigma_{t}a_{t}\right)\left(\sum_{t=j-k}^{j}T(k,j-t)\,\sigma_{t}b_{t}\right)\end{array}$$

Full text

1 Introduction

The fast multiplication algorithm discovered by Anatoly Karatsuba in 1960 (and published two years later) is known to be the oldest algorithm faster than the “grade school” method (when involved numbers or polynomials are large enough); while newer algorithms are still faster for sufficiently large numbers or polynomials, it is still widely used today for multiplicating medium-sized numbers or polynomials.

Due to its recursive divide-and-conquer approach, implementing this algorithm with no care about various issues (mostly related to storage of the temporary data) will lead to poor and often slow programs. Furthermore, the triple recursion involved by the algorithm, along with propagating changes in the computed data due to consecutive subtractions, makes implementing it in an iterative style more challenging.

The purpose of this paper is to deeply rewrite Karatsuba’s formula in such a way that an iterative implementation would at first glance naturally arise. The last example of high-level pseudocode given in the paper thus relies on a single simple loop. The current paper will rather focus on identifying and writing down the new formula rather than on computational issues.

While keeping the very same number of elementary products, namely $n^{\log_{2}3}$ multiplications, where $n$ is the degree of the polynomials, the formulas that will be presented in section 5 do not necessarily keep also the number of additions and subtractions as low as in the state-of-the-art implementations of the algorithm — deeper study of the question has not been made however.

2 Karatsuba’s recursion tree

Let $A$ and $B$ be two polynomials in the same indeterminate $x$; the divide-and-conquer paradigm to be used in the next sections requires splitting both $A$ and $B$ with the help of a third polynomial $X$ (most likely a monomial) in the same indeterminate $x$; the exact purpose of $X$ will be discussed in the current section as well as the following one.

While the ideas discussed here may be applied to any variant of Karatsuba’s initial formula, we take the following one as a starting point and group all terms as factors around each of its three distinct branches:

[TABLE]

with $A=A_{1}X+A_{0}$ and $B=B_{1}X+B_{0}$. The formula is intended to be applied recursively until elementary products are encountered, and at each step the five terms $A_{0}$, $A_{1}$, $B_{0}$, $B_{1}$ and $X$ must be redefined according to the exact level of the recursion. Each node in the recursion tree will be labelled here according to a radix-3 labelling system — by aggregating the arbitrary reference numbers, which are specified in (1), of the successive branches leading to it. Computing the product $AB$ now amounts to summing all contributions corresponding to the leaf nodes of a tree.

When a node is reached (by starting from the root node) without walking on any branch-2 nodes, we call it here a direct node; it will otherwise be called indirect. Obviously, direct nodes are reached by walking along paths whose label does not contain any digit $2$. An indirect node whose parent is a direct node will be called primary indirect. The following diagram illustrates that by showing direct and primary indirect nodes in a ternary recursion tree: hatched nodes are the primary indirect ones while all others are direct nodes (furthermore, shaded nodes are direct leaf nodes). Non-primary indirect nodes are discarded here, since they will be taken into account during the recursive process applied to each primary indirect node.

012012012012012012012

Obviously, computing the whole product can also be seen as summing the contributions of all direct leaf nodes and primary indirect nodes (which are either shaded or hatched on the diagram above). This strategy will be deeper studied in the section 4.

3 The art of splitting polynomials

Karatsuba’s algorithm seems to be most of the time implemented or studied by using the same splitting scheme which involves taking apart terms of lower and higher degree (this may intuitively be seen as a left/right approach); it will be referred to here as the “traditional splitting scheme”.

Several other splitting schemes are fully compliant with formula (1), namely any scheme taking apart groups of some power-of-$2$ sequential terms. The simplest one will be considered from now on: taking apart terms of even and odd rank. Such choice will have two main benefits: identifying the exact labelling number of a node is now easier and applying the algorithm to infinite power series will also be possible. This splitting scheme will be referred to here as the “interleaved splitting scheme”. Applying it recursively is illustrated below:

[TABLE]

New iterations of the interleaved splitting scheme lead to increasingly-sparse polynomials: each next term has initially a non-null coefficient, then each second term, then each fourth term, etc. If the term $X$ in the formula (1) is initially some $x$, then it will become $x^{2}$ at the second iteration, then $x^{4}$, etc.

Lemma 3.1**.**

Let $K$ be a direct leaf node, built according to the interleaved splitting scheme, and reached by following the $(d_{1}d_{2}d_{3}\dots d_{m})_{3}$-labelled path (where all $d_{k}$ are digits from $\{0,1\}$ since a direct node has no digit $2$ in its path and with $m=\log_{2}n$). The contribution of this node $K$ to the whole summation is

[TABLE]

which means that an initial radix-3 string is merely read later as a binary string with no further conversion (other than reverting the order of the digits).

Proof.

The left factor is the same for all direct leaf nodes; it does not rely on the path to the node $K$ but only on the number $m$ of iterations. Its purpose is to accumulate all iterated $X+1$ factors from the formula (1), namely

[TABLE]

The fraction is a shorthand notation for this expression111The $1/(1-x)$ part is an usual generating function for $1+x+x^{2}+x^{3}+x^{4}+\dots$ and multplicating it by $(1-x^{2^{m}})$ allows to truncate the series to an arbitrary degree..

The $x^{r}$ part, which is the fourth factor in the formula from the lemma, comes by following the given path made of selected branches 0 and branches 1 in different levels of the tree. Accumulating the iterated $X$ factors from the branch 1 in the formula (1) again is done as:

[TABLE]

where the digit $d_{k}$ is used for indicating whether the $x^{2^{k-1}}$ factor is accumulated or not — which is only the case in a branch 1, allowing to use the digit as an exponent. Of course, the whole product is equal to $x^{r}$.

Finally, the $a_{r}$ and $b_{r}$ factors are identified by induction: in the formula (1) $A$ and $B$ are increasingly-sparse polynomials having both $a_{k}$ and $b_{k}$ as their constant term for some arbitrary direct node of the tree, where $k$ also is the degree of all previously accumulated $X$ factors. This is obviously true at the root level of the tree (with $k=0$), and this property remains true at each new level of the tree whatever the selected branch is222This property can also be noticed in the illustration of the interleaved splitting scheme at the beginning of the current section for a polynomial of degree $7$, formally applying this splitting scheme being deeply related to applying the formula (1) in regards to the degree of computed terms. In this example, the constant term of each sparse polynomial in parentheses at any level is $a_{k}$ when the external factor of this polynomial is $x^{k}$.. Thus, the selected terms in $A$ and in $B$ are the coefficients of the terms of degree $r$ in each polynomial. ∎

4 Partially flattening the recursion tree

In this section, we compute all direct leaf nodes at once while gathering separately all primary indirect nodes (thus taking care of all shaded or hatched nodes in the diagram from the section 2).

The Lemma 3.1 helps achieving the first part of this goal: gathering all $a_{r}b_{r}x^{r}$ terms is done by computing the termwise product of $A$ and $B$ — which will be written down as $A(x)\odot B(x)$ — for both polynomials and series. Furthermore, the following conventions will be used from now on:

• •

explicit symbols for multiplication and convolution ($\times$ and $\ast$ respectively) for formulae which are applied recursively;

• •

an implicit notation if one of the factors is a polynomial — or the generating function of an integer sequence — containing only [math] and $1$ coefficients.

Let $A$ and $B$ be two polynomials of degree $n-1$ (for clarity, $n$ being a power of $2$). Once Karatsuba’s formula has been applied repeatedly $\log\left(n\right)$ times, matching coefficients in $A$ and $B$ are multiplied together for all direct leaf nodes of the recursion tree. We can group the $n$ elementary multiplications as:

[TABLE]

where, from a computational point of view, the left factor does not involve a true multiplication but rather $n$ shift/add operations — which match the required additions in Karatsuba’s formula.

We now have to gather all primary indirect nodes, as stated at the end of the section 2 and at the beginning of the current section. There are $2^{k}$ of them on the level $k$ of the tree ($k=0$ corresponding to the first level under the root node) as can be noticed on the diagram from the section 2, and we will collect them in the so-called level-order333A tree is traversed in level-order when each node on a level is visited before going to a lower level. Of course, we only care here about primary indirect nodes..

By using two nested summation symbols, we can “iterate” on all levels (there are $\log_{2}(n)$ of them) and then on all $2^{k}$ primary indirect nodes in the $k^{\textrm{th}}$ level.

We now have to remember that $A$ and $B$ are sparse polynomials, splitting them according to the interleaved splitting scheme and then subtracting the lowest part from the highest is equivalent to the two consecutive steps:

• •

multiplicating $A(x)$ and $B(x)$ by $1-x^{2^{k-1}}$;

• •

keeping only terms of degree $x^{0\cdot 2^{k}}$, $x^{1\cdot 2^{k}}$, $x^{2\cdot 2^{k}}$, $x^{3\cdot 2^{k}}$, etc.

The second task involves filtering a polynomial in order to keep specific terms and cancel all others; this is easily achieved by performing the termwise product of such polynomial with a relevant mask being another polynomial whose coefficients are in $\{0,1\}$. Building this mask is achieved by using some convenient tools from the theory of generating functions: it has to be remembered first that $1/(1-x^{2^{k}})$ expands to $1+x^{2^{k}}+x^{2\cdot 2^{k}}+x^{3\cdot 2^{k}}+\dots$, which can be truncated and shifted as required with:

[TABLE]

where $x^{2^{k-1}+j}$ gives the relevant “offset” for the $j^{\textrm{th}}$ primary indirect node at the $k^{\textrm{th}}$ level of the tree. Again, the binary encoding of $2^{k-1}+j$ is closely related to the radix-3 string labelling the path for reaching the parent of a given primary indirect node: we merely track all accumulated $X$ factors by iterating with the formula (1) before reaching the primary indirect node; more precisely, we want to match the term of lowest degree still available in $A$ and in $B$ (since many terms have been discarded through the splitting process).

In order to follow very closely the formula (1), each termwise product should immediately be divided by $x^{2^{k-1}+j}$ (before multiplicating the two newly-built polynomials) because we actually want each polynomial to have a constant term while separately accumulating the required $X$ factors; but the division can occur later since

[TABLE]

where the three factors match those in the branch 2 in the formula (1), except for one thing: while all $X$ factors have been correctly accumulated here, the $X+1$ ones (from the branches [math] and $1$ in the previous steps) are still missing — but there is exactly one of them at each level of the tree, as previously, and we can accumulate them with the same method than in the formula (2).

Thus, gathering all primary indirect nodes finally gives:

[TABLE]

where $k$ is the level of each considered row of nodes, and $j$ the rank of each direct subtracting node on the level $k$.

Since $2^{k-1}+j$ merely iterates over $1,2,3,\dots,n-1$, it is easy to use a single summation for directly iterating over all the considered subtrees:

[TABLE]

Summing both parts (2) and (3) results in $A\left(x\right)\times B\left(x\right)$ and is more or less equivalent to Karatsuba’s algorithm from a computational point of view — as long as some variable substitution is done before each recursive call in order to map sparse polynomials to new polynomials of smaller degree.

Using the previously described interleaved splitting scheme now allows to apply Karatsuba’s recursive formula to infinite power series (which is not the case with the traditional splitting scheme).

In the formulas (2) and (3), all $1-x^{n}$ numerators in the generating functions are intended to truncate periodical sequences of unitary and null coefficients to the required length. Extending these formulas to infinite power series is then easily achieved by removing such numerators:

[TABLE]

As a separate question, we may wonder for which values the index of summation $m$ in the formulas (3) and (4) will contribute to the computation of a term of degree $d$ in the final result. Otherwise said, we want to gather the primary indirect nodes involved in a given resulting term — these primary indirect nodes are labelled according to the previously specified enumeration $1,2,3,\dots,n-1$ by following the level-order (on the diagram of the section 2, the seven hatched nodes will be labelled $1,2,3,\dots,7$). The set $S_{d}$ of all such indices $m$ is:

[TABLE]

which comes directly from the formula (3): the $m^{\textrm{th}}$ primary indirect node have some $x^{m}P(x)$ contribution — with some $P(x)=c_{0}+c_{1}x^{2^{k+1}}+c_{2}x^{2\cdot 2^{k+1}}+\dots$ (at the $k^{\textrm{th}}$ level of the tree, and of course $k=\lfloor\log_{2}m\rfloor$). Because of all $X+1$ accumulated factors, each computed coefficient will also be shifted $2^{k}$ times “to the right”. The definition (5) gathers all primary indirect nodes such that one term in the involved sparse polynomial has a degree “close” enough to $d$.

The cardinality $|S_{d}|$ of such sets of indices is empirically found to be the sequence A268289 in the On-Line Encyclopedia of Integer Sequences, namely the cumulated differences between the number of digits 1 and the number of digits 0 in the binary expansions of integers up to $d$ . Another explicit expression for $|S_{d}|$ resorting to the $\tau$ Takagi function can also be given:

[TABLE]

with $d$ some non-negative444By convention ${\texttt{A268289}}_{0}=0$, the formal definition of the sequence being slightly different than the plain english one above which would imply the wrong statement ${\texttt{A268289}}_{0}=-1$. integer and $k=\lfloor\log_{2}(d)\rfloor$.

While inserting an extended proof of the previous identity would be far beyond the scope of the current paper, a quick hint will help building such a proof: when terms of the three sequences for all indices up to $2^{s}-1$ are known, we build the following terms up to the index $2^{s+1}-1$ with the help of the same building rule for the three sequences:

[TABLE]

and since the three sequences share the same initial terms, we finally prove that they are identical.

5 Fully flattening the recursion tree

Having described in the previous section how to handle two branches of the tree at once by using the termwise multiplication formula, we now go one step further. The key ideas from the previous section obviously are:

• •

handling several leaf nodes as a whole by using termwise products of polynomials;

• •

performing all required operations (additions, subtractions, shifts) by accumulating factors finally expanding as polynomials with unitary and null coefficients;

• •

keeping or cancelling coefficients of polynomials by performing the termwise product of the latter with relevant masks.

We will extend the application of these ideas in order to handle all the $n^{\log_{2}3}$ leaf nodes by computing only $n$ termwise products. As a starting point, we consider again the interleaved splitting scheme which is illustrated in Section 3 by the following example:

[TABLE]

The last line of the example shows four sparse polynomials which would occur in four different paths, namely $(11)_{3}$, $(10)_{3}$, $(01)_{3}$ and $(00)_{3}$ (when reading these paths in the reversed order as binary encoded strings, we can identify the index of the constant term for all these polynomials). Instead of exploring again the branches 0 and 1, we now explore the branch 2 from each of these four nodes — reaching the lowest level of the tree with four elementary terms:

[TABLE]

which have to be multiplicated with corresponding $b$ coefficients. We want to compute the four multiplications as $\left(1-x^{4}\right)A(x)\odot\left(1-x^{4}\right)B(x)$. Extraneous coefficients are easily cancelled with a mask. While the mask should be applied before the termwise product from a computational point of view, we focus rather on building the most concise formula and the use of the mask will be postponed.

It has to be noticed that the $1-x^{2^{k-1}}$ factor is the same for all nodes from the $k^{\textrm{th}}$ level; of course such factors may also be accumulated when walking on several branches $2$ on a given path. Since all nodes must be visited by the recursion process, all possible selections of such factors have to be considered. The ternary initial tree built from the formula (1) now becomes a more classical binary tree: at each level we can choose between either accumulating a new $1+X$ factor (to be used after the termwise product) or accumulating a new $1-X$ factor (to be used before the termwise product). Indeed, walking on a branch [math] or $1$, implies accumulating the relevant $1+X$ factor according to the initial identity (1), while walking on the branch $2$ implies subtracting coefficients to other ones wich is performed here by accumulating the relevant $1-X$ factor. The $1-X$ factors have an arithmetical purpose and must be applied on the actual numerical values before the termwise multiplication (which is going to definitely discard some coefficients), while the $1+X$ factors have a shifting/adding purpose and must be applied after the termwise purpose when the unwanted coefficients have been discarded.

Accumulating various $1+X$ factors actually plays two different roles in the computation: one has been already described (shifting and adding some terms), we focus now on the other one. The mask to be applied for discarding extraneous and useless subtracted coefficients happens to be the very same polynomial made of accumulated $1+X$ factors as long as we shift and truncate it accordingly.

This is proved as follows; when computing the product

[TABLE]

and considering a specific factor $1+x^{2^{k}}$, we have to remember that the latter gives some control on each block of $2^{k+1}$ consecutive coefficients: either the left pattern of $2^{k}$ coefficients is duplicated on the right part or not; for instance, removing the single factor $1+x^{2}$ would lead to the sequence of coefficients $1,1,0,0,1,1,0,0,1,1,0,0,\dots$ When “tracking” the coefficients during the recursion process (according to the interleaved splitting scheme), the $1-x^{2^{k}}$ factor means subtracting the left part of such a block from the right part, and of course we must cancel half of the coefficients in each block.

We can now write down a formula for the whole tree by choosing some way to iterate over all selections of $1-X$ factors or over all selections of $1+X$ factors; this could be done, for instance, by iterating over all subsets of $\{1,2,\dots,\log_{2}n\}$. By noticing that the mask is being iterated over all divisors of an appropriate polynomial, we can utilise some extra conventions and notations to arrive at an elegant formulation:

[TABLE]

with $\dot{f}$ selecting all unselected $1+X^{2^{k}}$ factors in $f$ and negating — for each one — the coefficient of their non-constant term, and $\dot{f}^{\textrm{\tiny\hskip 1.0ptW}}$ the leading term (including the coefficient) of $\dot{f}$. The superscript character W stands for “weight”555This symbol is compact but not very common; it can be found however in an article by Shigeru Kuroda, Shestakov-Umirbaev reductions and Nagata’s conjecture on a polynomial automorphism (2007)..

The divisors of the $\sum_{k=0}^{n-1}X^{k}$ polynomial are of course all possible selections of factors in $(1+X)(1+X^{2})(1+X^{4})\dots$ The $\dot{f}^{\textrm{\tiny\hskip 1.0ptW}}$ term is used for shifting the mask at the beginning of the meaningful subtracted terms, and the coefficient of this term is either $-1$ or $1$ according to the number of subtractions (respectively odd or even) in the current path.

Of course the expected $n^{\log_{2}3}$ elementary multiplications are embedded in the previous formula, since the whole idea of the current section was to track them and to shift/add/subtract them according to the initial (1) formula; no supplementary multiplication was added anywhere.

An alternate version of the formula comes from the fact that $\dot{f}^{\textrm{\tiny\hskip 1.0ptW}}$ is either some $-x^{k}$ or $x^{k}$ with $k$ being the sum of powers of $2$ involved in the product of $1-x^{2^{j}}$ when building the $\dot{f}$ polynomial. We can thus iterate over $1,2,\dots,n-1$ as the degree of $\dot{f}^{\textrm{\tiny\hskip 1.0ptW}}$ and build other terms from it:

[TABLE]

with $\dot{f}_{k}=(1-X)^{d_{0}}(1-X^{2})^{d_{1}}(1-X^{4})^{d_{2}}\dots$ by referring to the binary digits of $k=(\dots d_{2}d_{1}d_{0})_{2}$, with also $\bar{f}_{k}=(1+X)^{1-d_{0}}(1+X^{2})^{1-d_{1}}(1+X^{4})^{1-d_{2}}\dots$ where the bar is intended to show that the product is truncated to the same “format” than $\dot{f}_{k}$, made of $\log_{2}n$ different factors, despite the infinite number of leading zeros in the binary encoding of $k$, and with $\sigma_{k}={\texttt{A106400}}_{k}$ in the On-Line Encyclopedia of Integer Sequences ($\sigma_{k}=-1$ if the binary weight of $k$ is odd and $\sigma_{k}=1$ otherwise). Several expressions for $\sigma_{k}$ are published on the page of the sequence A106400 — one involving an hypergeometric ${}_{2}F_{1}$ function.

Formula (8) can now easily be adapted to infinite power series as:

[TABLE]

by merely removing the bar from $\bar{f}_{k}$, now defining $f_{k}$ as an infinite product.

6 Implementing the new formula

Implementing the formula (7) with no care about the true purpose of each part does not lead to a very efficient code: some elementary multiplications would be done though they are going to be cancelled soon after, many operations involving null coefficients in very sparse polynomials could be avoided, etc. But a first attempt can be given as a proof of concept and we give below two pieces of pseudocode intended to be used with any computer algebra system handling the polynomial type; they do not focus on low-level implementation issues (how more or less sparse polynomials are internally represented in order to give the most efficient access to their coefficients). Termwise multiplication of polynomials should of course be already implemented.

Two polynomials $f$ and $f^{\prime}$ are used below for accumating shift/add and shift/subtract operations according to what was previously discussed.

The following code shows how subtracting and masking factors are accumulated while iterating on the branches of the tree; the key idea is to use the $1+X$ factors for selecting subtracted terms as well as for propagating them and the $1-X$ corresponding factors for performing the subtractions: {codebox} \Procname$\proc{Multiply}(A,B)$ \li$d\leftarrow\lceil\log_{2}(1+\max(\deg A,\deg B))\rceil$ \li$n\leftarrow 2^{d}$ \li$s\leftarrow 0$ \li\For$k\leftarrow 0$ \To$n-1$ \li\Do$f\leftarrow 1$ \li$f^{\prime}\leftarrow 1$ \li\For$j\leftarrow 0$ \To$d-1$ \li\Do\If$k\mathrel{\&}2^{j}\neq 0$        \Commenttest if bit $j$ of $k$ is set \li\Then$f^{\prime}\leftarrow(1-X^{2^{j}})\,f^{\prime}$       \Commentsubtracting factor \li\Else$f\leftarrow(1+X^{2^{j}})\,f$       \Commenttermwise mask \End\End\li$\id{degree}\leftarrow\deg f^{\prime}$          \Commentdegree of $f^{\prime}$ \li$\id{lt}\leftarrow\proc{coeff}(f^{\prime},X,\id{degree})\,X^{\id{degree}}$          \Commentleading term in $f^{\prime}$ \li$s\leftarrow s+f\,(f\id{lt}\odot\,f^{\prime}A\odot f^{\prime}B)$ \End\li\Returns

where it can be seen that the variable $f$ has two distinct purposes: the algebraic one and also a tracking purpose for identifying how many terms have to be kept in the termwise multiplication, as explained in the previous section.

The previous pseudocode however loses all benefits of traditional implementations of Karatsuba’s algorithm because the same subtractions are computed for distinct values of $k$. Fortunately iterating over the binary expansions of $k$ by using the reflected binary code (Gray code) instead of the standard radix-$2$ labelling system preserves the required number of subtractions.

Since iterating over such binary expansions (the Gray code) implies flipping a single bit between two consecutive integers, we are now visiting the nodes in a new unintuitive order in such a way that the accumulated factors can be reused from a node to another one by performing two single changes only: dividing the current product by some $1-X^{k}$ or by some $1+X^{k}$ (according to the exact location of the flipped bit) in order to “cancel” the branch being left, as well as multiplicating it by some $1+X^{k}$ or by some $1-X^{k}$ (according to the exact location of the flipped bit) in order to “enable” the new visited branch.

Of course, polynomial divisions are exact here since we merely divide some $(1+X)(1+X^{2})(1+X^{4})\dots$ polynomial by one of its previously “enabled” factors for cancelling it — the same for another $(1-X)(1-X^{2})(1-X^{4})\dots$ polynomial. Actually, both polynomial are always kept synchronized in the following way: if some $(1\pm X^{k})$ factor is “enabled” in one polynomial, it is disabled in the other one (because both branches are not visited together at the same level).

The following version, though not optimized by itself from an implementation point of view (because it still relies on high-level polynomial types), gives the prototype of a more optimized iterative version of Karatsuba’s algorithm: {codebox} \Procname$\proc{Multiply2}(A,B)$ \li$d\leftarrow\lceil\log_{2}(1+\max(\deg A,\deg B))\rceil$ \li$n\leftarrow 2^{d}$ \li$f\leftarrow(1-X^{n})/(1-X)$            \Commentinitial mask (all bits set) \li$f^{\prime}\leftarrow 1$ \li$g\leftarrow 0$            \CommentGray-code counterpart of $k$ \li$s\leftarrow f\,(f\odot A\odot B)$            \Commentcase $g=k=0$ \li\For$k\leftarrow 1$ \To$n-1$ \li\Do$j\leftarrow\lfloor\log_{2}(k\;\textrm{xor}\;k-1)\rfloor$          \Commentleast significant set bit in $k$ \li\If$g\mathrel{\&}2^{j}=0$         \Commenttest if bit $j$ has to be set in $g$ \li\Then$f^{\prime}\leftarrow(1-X^{2^{j}})\,f^{\prime}$        \Commentsubtracting factor \li$f\leftarrow f/(1+X^{2^{j}})$        \Commenttermwise mask \li\Else$f^{\prime}\leftarrow f^{\prime}/(1-X^{2^{j}})$        \Commentsubtracting factor \li$f\leftarrow(1+X^{2^{j}})\,f$        \Commenttermwise mask \End\li$g\leftarrow g\;\textrm{xor}\;2^{j}$         \Commentupdate $g$ (Gray-code of $k$) \li$\id{degree}\leftarrow\deg f^{\prime}$         \Commentdegree of $f^{\prime}$ \li$\id{lt}\leftarrow\proc{coeff}(f^{\prime},X,\id{degree})\,X^{\id{degree}}$         \Commentleading term in $f^{\prime}$ \li$s\leftarrow s+f\,(f\id{lt}\odot\,f^{\prime}A\odot f^{\prime}B)$ \End\li\Returns

A lower-level implementation of this pseudocode should avoid actually storing the $f^{\prime}$ polynomial in a separate buffer and computing both $f^{\prime}A$ and $f^{\prime}B$ products — the idea being rather to directly store $f^{\prime}A$ and $f^{\prime}B$, and merely update them at each step of the loop.

Furthermore, efficiently implementing the previous pseudocode should take care of the subtracting and adding steps: since polynomials become very sparse for some values of $k$, very few terms should be manipulated at these points. Two main directions should be explored for that purpose: using linked lists for representing polynomials or tracking the remaining non-null coefficients by using an elaborated system of strides666This is one of the most important concepts behind the famous Numpy module for Python; strides allow to build views on parts of an existing array without actually copying it.. Elementary multiplications should of course be aware of the mask to be applied in order to avoid useless computation.

7 Computing arbitrary coefficients

Adapting the formula (9) to a more or less Cauchy-like one is actually achievable — though very inefficient from a computational point of view since the same multiplications will be performed again and again rather than propagated through the shift/add process. It may however have some interest in further theoretical investigations.

For that purpose, we do not need the $\odot$ termwise operator any longer since we have no interest in computing several coefficients at once, but we now need the $\&$ bitwise multiplication operator (the bitwise “and” operator) since the new formula will highly rely on testing whether such or such $1\pm X$ factor is selected or not.

Let $g(x)=a_{0}+a_{1}x+a_{2}x^{2}+\dots$ and $h(x)=b_{0}+b_{1}x+b_{2}x^{2}+\dots$, then

[TABLE]

with $\tau_{k}(m,j)=[k\mathbin{\&}(m-j)=0]$ and $\upsilon_{k}(j,t)=[k\mathbin{\&}(j-t)=j-t]$, both defined by using Iverson bracket, and $\sigma_{k}$ as previously defined in the formulas (8) and (9).

The $\tau_{k}$ function is used for testing whether a coefficient from a given degree will actually be shifted and added as a contribution to another given degree; the $\upsilon_{k}$ function is used for testing whether a coefficient of a given degree will be shifted and subtracted when accumulating all $1-X$ factors.

After having noticed that $\tau_{k}$ is related to the sequence A047999 and that $\upsilon_{k}$ is related to the sequence A106344, we decide to clean up the previous formula by relying on A047999 only, namely on the Sierpiński triangle, finally getting:

[TABLE]

with $T(n,k)$ being A047999 — defined as $T(n,k)=[k\mathbin{\&}(n-k)=0]$. While the resulting nested summations may look rather heavy, a quick glance at the famous graphical representation of Sierpiński triangle will show that most of the involved terms are null.

8 Conclusion

While the section 6 highly relies on pseudocode parts, identifying and publishing the formulas (7) to (9) as new theoretical convolution formulas actually was the true purpose of the current paper. While the number of elementary multiplications remains identical to the expected $n^{\log_{2}3}$ one, actually applying the formula probably involves more additions and subtractions than what would be the case by following the conventional recursive approach — their exact amount being however not investigated here.

The final formula, in the previous section, shows that tracking individual coefficients through the whole recursion tree is achievable; as a sophisticated convolution involving the Sierpiński triangle, it may be seen as a starting point for further investigations focused towards combinatorics.

Conflict of Interest: The authors declare that they have no conflict of interest.

The current article is accessible on http://export.arxiv.org/pdf/1902.08982 .