Numbers in Base X Part 1: Rising and Falling Division

When I was first introduced to polynomial arithmetic, my teacher said it's very similar to regular arithmetic. In fact, because you don't need to worry about carrying, performing addition, subtraction, and multiplication, is significantly easier. Are polynomials really that similar to numbers in base $x$? I felt compelled to find out, unaware of the fact that I was about to go down an unfathomably deep rabbit hole. This article will introduce you to various unusual approaches to dividing both polynomials and integers, and the algorithms we'll develop will serve as a foundation for many future explorations.

Preliminaries

In school, we learn that the position of each digit in the standard representation of an integer corresponds to the power of the associated base, usually $10$, that multiplies it. For example: $123$ is really just $1 \cdot 10^2 + 2 \cdot 10^1 + 1 \cdot 10^0$. I'll omit $10^0$, which is just $1$, and the power of $10^1$ from here on for the sake of brevity. Let's subtract $29$ from $123$ using the expanded representation:

$$\begin{align*} 123 - 29 &= 1 \cdot 10^2 + 2 \cdot 10 + 3 - 2 \cdot 10 - 9\\ & = 1 \cdot 10^2 +(2-2) 10 +(3-9) \\ & = 1 \cdot 10^2 + 0 \cdot 10 - 6 \end{align*}$$

If we had $x$ instead of $10$ in this example, we'd be done with the subtraction. However, we can only use digits $0$ to $9$ in the base $10$ representation, and we have to deal with the $-6$ by borrowing from the next highest power:

$$\begin{align*} 123 - 29 &= 1 \cdot 10^2 + 0 \cdot 10 - 6 \\ & = (9 \cdot 10 + 1 \cdot 10) - 6 \\ & = 9 \cdot 10 + 10 - 6 \\ & = 9 \cdot 10 + 4 \\ & = 94 \end{align*}$$

Multiplication is also straightforward, provided we remember its distributive property, which tells us we can multiply each term inside one bracket by the terms in the other:

$$\begin{align*} 123 \cdot 29 &= (1 \cdot 10^2 + 2 \cdot 10 + 3) (2 \cdot 10 + 9)\\ & = 1 \cdot 10^2 (2 \cdot 10 + 9) + 2 \cdot 10 (2 \cdot 10 + 9) + 3 (2 \cdot 10 + 9)\\ & = 2 \cdot 10^3 + 9 \cdot 10^2 + 4 \cdot 10^2 + 18 \cdot 10 + 6 \cdot 10 + 27 \\ & = 2 \cdot 10^3 + 13 \cdot 10^2 + 24 \cdot 10 + 27 \end{align*}$$

Once again, we're required to normalise the representation, meaning every multiple of a power of $10$ has to be less than $9$. This is easily accomplished if we rewrite the "digits" as multiples of $10$ with a remainder:

$$\begin{align*} 123 \cdot 29 & = 2 \cdot 10^3 + 13 \cdot 10^2 + 24 \cdot 10 + 27 \\ & = 2 \cdot 10^3 + (10 + 3) 10^2 + (2 \cdot 10 + 4) 10 + (2 \cdot 10 + 7) \end{align*}$$

Multiply through and carefully regroup to get the final result:

$$\begin{align*} 123 \cdot 29 &= 2 \cdot 10^3 + (10 + 3) 10^2 + (2 \cdot 10 + 4) 10 + (2 \cdot 10 + 7) \\ & = (2+1)10^3 + (3+2)10^2 + (4 + 2)10 + 7 \\ & = 3 \cdot 10^3 + 5 \cdot 10^2 + 6 \cdot 10 + 7 \\ & = 3567 \end{align*}$$

I would've done this much faster using the notation and algorithms taught in elementary school, but it's not too difficult conceptually.

Division is where my younger self got stuck. Like many students, I had memorised the long division algorithm without understanding why it works, and was also vaguely aware that I could, in principle, subtract repeatedly. The latter option seemed way too inelegant, not to mention tedious, and the only way I could adapt standard long division to the expanded representation was to normalise after each multiplication and subtraction step. This really bothered me because I could perform all the other operations while treating the intermediary results as polynomials, only normalising at the end. What made the problem worse is that the polynomial division algorithm my teacher taught me required the use of rational numbers, and those are nowhere to be seen when dividing integers! I gave up on this approach to division eventually, but the urge to compare integers and polynomials would frequently resurface over the next 20 years. Enough with the navel-gazing, we've got quite a bit of work ahead of us before we get to play!

Polynomial Division

So how do we divide polynomials? The algorithms I was taught were pretty much identical to those presented on Wikipedia, but I'd like to offer an alternative perspective that, in my opinion, provides a lot more insight into what's happening.

Multi-Pass Division

Returning to our running example, this time in base $x$ and with reversed digits, let's try to divide $f(x):3x^2 + 2x + 1$ by $g(x):2x + 9$. The strategy will be to rewrite $x^2$ and $x$ as multiples of the divisor, which requires us to first subtract all the lower degree (power of $x$) terms of the divisor from itself to isolate the leading term (the highest degree). Once we have a single term, we can clear out its coefficient (the leading coefficient) and multiply by an appropriate power of $x$. Symbolically, the process looks like this:

$$\begin{align*} f(x): 3x^2 + 2x + 1 &= 3(\frac{x}{2}(2x+9-9))+2x+1 \\ &= \frac{3x}{2}(g(x)-9) + 2x + 1 \\ &= \frac{3x}{2}g(x) + (2 - \frac{27}{2})x + 1 \\ &= \frac{3x}{2}g(x) -\frac{23}{2}x + 1 \end{align*}$$

Take a moment to fully understand the symbolic manipulations because the entire algorithm is just a repeated application of this idea until we have a bunch of multiples of the divisor, and a remainder that has a lower degree. We still have one more term to rewrite:

$$\begin{align*} f(x) &= \frac{3x}{2}g(x) -\frac{23}{2}(\frac{1}{2}(2x+9-9)) + 1 \\ &= \frac{3x}{2}g(x) -\frac{23}{4}(g(x)-9) + 1 \\ &= \frac{3x}{2}g(x) -\frac{23}{4}g(x) + 1 + \frac{207}{4} \\ &= \frac{3x}{2}g(x) -\frac{23}{4}g(x) + \frac{211}{4} \end{align*}$$

This was the last step because the remaining term, $211/4$, has a lower degree than the divisor, so we ran out of powers of $x$ to rewrite. We can factor out the $g(x)$ to get our final result:

$$\begin{align*} f(x) = (\frac{3x}{2} -\frac{23}{4})g(x) + \frac{211}{4} \end{align*}$$

We effectively rewrote $f(x)$ as a multiple of $g(x)$, the quotient, with an added remainder polynomial of lower degree. Conceptually, this algorithm splits the divisor into the leading term and the lower degree terms, and, at each step, we slide the lower degree terms to the right, and subtract them from the corresponding numerator coefficients. Before the subtraction, each divisor coefficient must be multiplied by the coefficient of the rewritten power of $x$, and divided by the divisor's leading coefficient. We can bring the resulting numerator coefficient under a common denominator, and, to avoid having to keep track of the various powers of the leading divisor coefficient, we can pre-emptively multiply all numerator coefficients, below the modified terms, by the divisor.

Here's what this algorithm looks like when translated into python code:

from fractions import Fraction # rational numbers support

def polyDivMultiPass(num: list, div: list) -> tuple[list, list]:
    divLen = len(div)
    assert divLen > 0, 'An empty divisor list implies division by 0'
    divDeg = divLen - 1 # the divisor's degree

    leadingZeroTerms = None
    for i in range(divDeg, -1, -1):
        if div[i] != 0:
            leadingZeroTerms = divDeg - i
            divLen -= leadingZeroTerms
            divDeg -= leadingZeroTerms
            break
    assert leadingZeroTerms is not None, 'All divisor coefficients were zero'

    numDeg = len(num) - 1 # numerator degree
    # Not exiting early won't impact correctness, but it can
    # convert integers to rationals unnecessarily.
    if numDeg - divDeg < 0:
        return [], num

    divLc = div[divDeg] # the divisor's leading coefficient
    divLcAccum = divLc
    divLcAccumPrev = 1

    output = num.copy() # copy the numerator to preserve the input arguments

    for i in range(numDeg, divDeg - 1, -1): # decreasing sequence
        prevCoeff = output[i]

        for k in range(1, divLen): # skip the leading divisor coefficient
            # The value of i-k will produce a decreasing sequence of indices,
            # corresponding to the lower degree divisor terms that flow out 
            # out of the rewritten term. Multiplying by divLc is done to
            # bring the values under a common denominator.
            output[i - k] = divLc * output[i - k] - prevCoeff * div[divDeg - k]

        # Multiply all remaining coefficient by the leading divisor coefficient
        for k in range(i - divLen, -1, -1):
            output[k] *= divLc

        output[i] = Fraction(prevCoeff, divLcAccum)
        # Store the previous power of the leading divisor coefficient so we can
        # divide the remainders by it. Eliminating this store requires a 
        # restructuring of the loop that will make the code harder to read.
        divLcAccumPrev = divLcAccum
        divLcAccum *= divLc

    # We didn't get to the remainders in the loop above, so make sure to
    # divide them as well.
    for i in range(divDeg):
        output[i] = Fraction(output[i], divLcAccumPrev)

    return output[divDeg:], output[:divDeg] # split into quotient and remainder

def printQuotRem(quotRem: tuple[list, list]):
    print('Quotient: [',', '.join([str(p) for p in quotRem[0]]),']'
          ' Remainder: [',', '.join([str(p) for p in quotRem[1]]),']')

I added "multi-pass" to the name of the function because we're doing multiple passes over the output list before we compute all coefficients. The input polynomials are specified as coefficient lists in increasing power order, which might look odd if you're not a programmer, but it makes adding higher powers of $x$ more efficient. I've also added a little helper function to print out the output of polyDivMultiPass. Let's divide $6 x^5+5 x^4+4 x^3+3 x^2+2 x+1$ by $7 + 8 x + 9 x^2$ using our code:

printQuotRem(polyDivMultiPass([1, 2, 3, 4, 5, 6], [7, 8, 9]))

Quotient: [ 872/2187, -10/243, -1/27, 2/3 ] Remainder: [ -3917/2187, -1972/2187 ]

Translating the result back to mathematical notation:

$$\frac{2 x^3}{3}-\frac{x^2}{27}-\frac{10 x}{243}+\frac{872}{2187}, -\frac{1972 x}{2187}-\frac{3917}{2187}$$

If you install the SymPy python library (enter pip install sympy in a command prompt or terminal), you can validate the output against its polynomial division algorithm:

from sympy import *

x = var('x')
n = poly(1 + 2 * x + 3 * x**2 + 4 * x**3 + 5 * x**4 + 6 * x**5)
d = poly(7 + 8 * x + 9 * x**2)
quotient, remainder = div(n, d)
print(quotient.all_coeffs(), remainder.all_coeffs())

[2/3, -1/27, -10/243, 872/2187] [-1972/2187, -3917/2187]

If you're a programmer, try to spot the redundant computations and see if you can come up with a better implementation.

A Change in Perspective

Doing multiple passes over the coefficient list is a straightforward adaptation of how I presented the division algorithm, but we can also look at the computations from the perspective of a single output coefficient! Let's say you wanted to compute the 3rd quotient coefficient directly, and, to make the analysis simpler, we can assume the divisor polynomial is monic, meaning its leading coefficient is $1$. What information do we need to compute it? Well, if the divisor has degree $1$, this quotient coefficient will obviously depend on the constant divisor term (the degree $0$ one) from the rewrite of the previous (2nd) quotient coefficient. In other words, we have to multiply the constant divisor coefficient by the previously computed quotient, and then subtract it from the 3rd highest numerator coefficient. Similarly, if the divisor was of degree $2$, the 3rd quotient coefficient would depend on the previous $2$ quotient coefficients, but, this time, the first (highest) coefficient needs to be multiplied by the constant divisor coefficient, and the 2nd by the degree $1$ divisor coefficient. In general, a quotient coefficient will depend on, at most, $n$ previous coefficients, where $n$ is the degree of the divisor minus $1$.

Feel free to use this idea to compute the quotient coefficients for the examples above if you're struggling to visualise what's going on. If it helps, what I see in my head is something similar to a linear recurrence, where an array of the negated lower degree divisor coefficients, expressed in increasing power order, are sliding against the numerator coefficients, in decreasing order; each step generates an entry in a quotient array, which is parallel to the numerator one. Where the quotient array and the sliding divisor overlap, we multiply the overlapping coefficient, sum all the products, and subtract the result from the corresponding numerator coefficient. If there are no overlapping entries in the arrays, you can assume the missing values are zero. Let's divide $[1, 2, 3, 4, 5, 6]$ by $[7, 8, 1]$ as an example:

$$\underline{5} - 8 \cdot 6 = -43$$

$$\underline{4} + 8 \cdot 43 - 7 \cdot 6 = 306$$

$$\underline{3} - 8 \cdot 306 + 7 \cdot 43 = -2144$$

The values inside the boxes are the computed coefficients, and the underlined numbers are the numerators we need to subtract from in each step. Unfortunately, this approach doesn't work for the remainder coefficients, besides the highest degree one, because we're not doing any rewriting for them. The solution turns out to be very straightforward if we apply the same analysis we used for the quotient. The constant degree term of the remainder would only be affected by the last rewrite step, so we need to subtract the product of the constant quotient and divisor coefficients. The degree $1$ remainder coefficient will depend on the last $2$ quotients, but this time the constant quotient coefficient is multiplied by the degree $1$ divisor coefficient, and the degree $1$ quotient coefficient by the constant divisor coefficient. At this point you might recognise that the steps are similar to the quotient computation, except we're sliding the divisor against the quotient array in reverse, and are not using the previous output in each step:

$$\underline{1} + 2144 \cdot 7 = \boxed{15009}$$

$$\underline{2} - 306 \cdot 7 + 2144 \cdot 8 = \boxed{15012}$$

Depending on your maths background, you might've also noticed that we're effectively performing a partial discreet convolution of the coefficient lists to produce the values we need to subtract from the numerator. You can verify this using the numpy library in python (pip install numpy if you don't have it):

import numpy
print(numpy.convolve([7, 8, 1], [-2144, 306, -43, 6]))

[-15008 -15010      3      4      5      6]

This isn't surprising given how we could, if we wanted to, compute the remainder as $\mathrm{numerator}-\mathrm{quotient} \cdot \mathrm{divisor}$, with the polynomial multiplication being a discrete convolution. The remainder algorithm above is simply an optimised version of this computation.

Modifying the monic algorithm to support non-monic divisors is also quite easy, because we can factor out the leading coefficient of the divisor, and bring it into the numerator:

$$\begin{align*} \frac{3x^2 + 2x + 1}{2x + 9} &= \frac{3x^2 + 2x + 1}{2 (x + \frac{9}{2})} \\ &= \frac{\frac{3}{2}x^2 + \frac{2}{2}x + \frac{1}{2}}{x + \frac{9}{2}} \end{align*}$$

The problem is now reduced to monic division with rational coefficients. We can use python's Fraction class to convert all coefficients to fractions before passing them to a monic division function (or implement that function so it uses rationals internally), but I find it more interesting to create an algorithm similar to polyDivMultiPass, which computes the output numerators directly, and only generates the appropriate rational coefficients at the very end.

Using the multi-pass algorithm for inspiration, it's obvious that each numerator should be multiplied by an increasing power of the leading divisor coefficient, which we can label $L$. We must also bring the previous coefficients under a common denominator, and, because each of their denominators is a consecutive power of $L$, we can simply multiply each one by an increasing power of $L$, starting from $0$ for the last rewritten term, and going as far back $L^{k-1}$, where $k$ is the degree of the divisor. Note how each coefficient of the divisor will always be multiplied by the same power of $L$, meaning we can precompute and cache their product.

The numerators in the remainder computation should all be multiplied by the same power of $L$, and the subtracted quotient numerators should be adjusted in the same way as the quotient computation, except the order is reversed. Unfortunately, precomputing the multiplier doesn't work for remainders, except for the leading coefficient. Bringing it all together, we get the following python function:

def polyDivRecurrence(num: list, div: list) -> tuple[list, list]:
    divLen = len(div)
    assert divLen > 0, 'An empty divisor list implies division by 0'
    divDeg = divLen - 1 # the divisor's degree

    leadingZeroTerms = None
    for i in range(divDeg, -1, -1):
        if div[i] != 0:
            leadingZeroTerms = divDeg - i
            divLen -= leadingZeroTerms
            divDeg -= leadingZeroTerms
            break
    assert leadingZeroTerms is not None, 'All divisor coefficients were zero'

    numDeg = len(num) - 1 # numerator degree
    diffDeg = numDeg - divDeg
    if diffDeg < 0:
        return [], num

    divLc = div[divDeg] # the divisor's leading coefficient
    # Precompute the required powers of the leading divisor coefficient and the
    # scaled lower degree divisor coefficients to keep it consistent with the 
    # pen-and-paper notation. Consumes more memory but avoids redundant computations.
    divLcAccum = divLc
    divLcPowers = [1, divLc]
    for i in range(diffDeg):
        divLcAccum *= divLc
        divLcPowers.append(divLcAccum)

    divDegMinus1 = divDeg - 1
    divScaled = [div[divDegMinus1]] # the divisor must be non-empty as a precondition
    for i in range(divDeg - 2, -1, -1): # avoiding list comprehensions for clarity
        divScaled.append(divLcPowers[divDegMinus1 - i] * div[i])

    quotientLen = diffDeg + 1
    # pre-allocate the lists, because we'll fill them backwards
    quotientNums = [0] * quotientLen
    quotients = [0] * quotientLen

    for i in range(numDeg, divDeg - 1, -1):
        iRev = numDeg - i # incrementing index; could `enumerate` the range too
        # Fill the quotient backwards to avoid having to reverse the list; the
        # index is off by 1 to avoid unnecessary subtractions in the next loop.
        quotientIndex = quotientLen - iRev 

        quotNum = num[i] * divLcPowers[iRev] # scale the associated numerator coefficient
        for k in range(0, min(iRev, divDeg)):
            quotNum -= divScaled[k] * quotientNums[quotientIndex + k]

        quotientIndex -= 1 # the quotient index is off by one here, fix it
        quotientNums[quotientIndex] = quotNum
        quotients[quotientIndex] = Fraction(quotNum, divLcPowers[iRev + 1])

    remainders = []
    for i in range(divDeg):
        remNum = num[i] * divLcAccum
        for j in range(i + 1):
            # The values of quotientNums[j] * divLcPowers[j] could be precomputed if the
            # remainder degree tends to be very large. Would use more working memory.
            remNum -= quotientNums[j] * divLcPowers[j] * div[i-j]
        remainders.append(Fraction(remNum, divLcAccum))

    return quotients, remainders

Pen and Paper

If you don't have access to a computer, I find that the notation inspired by this alternate algorithm is much more convenient when dividing large polynomials. The pen-and-paper example should also help clarify what the code does. Let's divide $6 x^5+5 x^4+4 x^3+3 x^2+2 x+1$ by $9 x^2+8 x+7$ as an example. The same "slide the divisor against the quotient coefficients" idea still applies, but we need to scale by the appropriate powers of $L$ in each step. Start by computing all powers of $L=9$ needed for the division, and list them in increasing power order. The degrees of the numerator and divisor are $5$ and $2$ respectively, meaning we'll need $5 - 2 + 1 = 4$ powers of $9$:

9, 81, 729, 6561

If you have access to a calculator, you can type $9 * 9$, and then press the equals key $3$ times to get all the powers of $9$. Next, list all the negated lower degree coefficients of the divisor in decreasing power order, then multiply each by the power of 9 above and one entry to left - the first value will be multiplied by 1 implicitly. This means we get $-8 \cdot 1$ and $-7 \cdot 9 = -63$. You can do it all in one step, placing the results on the next line, but this time list them in increasing power order to avoid getting disoriented once we start calculating the quotient coefficients. To wrap up the preparations, you should bring the highest degree numerator coefficient down to the third line, resulting in the following:

9, 81, 729, 6561
-63, -8
6

Our first quotient coefficient is $6/9 = 2/3$, and to compute the next one we evaluate one step of our recurrence relation by multiplying the next highest term of the numerator, which is $5$, by the lowest power of $9$ in the first row, and then adding the result of the linear recurrence formed by the 2nd and 3rd row. To avoid having to write a bunch of leading zeros in the first row, you can implicitly multiply by $0$ by ignoring the $-63$ entry. This gives us the next value in the third row: $5 \cdot 9 - 8 \cdot 6 = -3$:

9, 81, 729, 6561
-63, -8
6, -3

This makes $-3/81 = -1/27$ our next quotient coefficient. Continue by bringing down the $4$ from the numerator, multiply it by the next power of $9$, which is $81$, and add the result of the linear recurrence: $4 \cdot 81 + 8 \cdot 3 - 63 \cdot 6 = -30$:

9, 81, 729, 6561
-63, -8
6, -3, -30

The third quotient coefficient is $-30/729$, or $-10/243$. Identifying the inputs for the linear recurrence is quite easy, you start from the last computed value in the third row, and pair it with the corresponding value from the second row. For the last quotient coefficient, that's $8 \cdot 30 + 3 \cdot 63$, which we add to $3 \cdot 729$ to get the last quotient numerator, $2616$:

9, 81, 729, 6561
-63, -8 
6, -3, -30, 2616

Don't forget to divide by 6561 to get the final quotient coefficient: $2616/6561 = 872/2187$. I like to underline the denominator value from the first row after each coefficient computation so I don't lose track of where I am.

Computing the leading coefficient of the remainder can be done using the same process as the quotient coefficients. Bring down the 2 from the numerator and add the recurrence: $2 \cdot 6561 - 8 \cdot 2616 + 63 \cdot 30 = -5916$. I like to place the result in line with the quotient, separated by a vertical bar:

9, 81, 729, 6561
-63, -8
6, -3, -30, 2616    |    -5916 

Divide by 6561 to get $-5916/6561 = -1972/2187$.

I visualise the remainder computation as a table where the first row comes from the low degree coefficients of the quotient, as many of them as there are lower degree terms in the divisor. The second row is comprised of powers of the leading divisor coefficient in decreasing power order, and the third contains the negated lower degree divisor coefficients in increasing power order.

-30, 2616
  9,    1
 -7,   -8

Computing a coefficient of the remainder is now a matter of multiplying all values in the same column and summing the results: $30 \cdot 9 \cdot 7 - 2616 \cdot 1 \cdot 8 = -19038$ This value is then added to the corresponding scaled numerator coefficient: $2 \cdot 6561 - 19038 = -5916$. It just so happens that this computation coincides with what we did for the quotient. The next coefficient is computed by shifting the bottom row to the right, discarding the coefficient that went out of the table and filling the new empty spot with 0:

-30, 2616
  9,    1
  0,   -7

Notice how the first two rows don't chage. This means the values can be pre-computed! Carry out the same procedure by bringing down the last numerator coefficient: $1 \cdot 6561 - 2616 \cdot 1 \cdot 7 = -11751$:

9, 81, 729, 6561
-63, -8
6, -3, -30, 2616    |    -5916, -11751 

The constant term of the remainder is therefore $-11751/6561 = -3917/2187$. If the divisor degree was higher, we would've kept shifting the bottom row to the right until we filled it with zeros. With a little bit of practice, you can do this without explicitly constructing the table, though I prefer to write the pre-multiplied quotients in the 4th row, because I find it less error-prone when the remainder has a large degree:

9, 81, 729, 6561
-63, -8
6, -3, -30, 2616    |    -5916 , -11751
-270, 2616

Your phone's calculator can make quick work of each step in the computation, or you can use a separate sheet of paper to keep the bookkeeping notation above nice and tidy. I find this notation significantly easier to work with compared to synthetic division, though, nowadays, I mostly let the computer do the work for me.

If you're planning to teach this notation to someone else, I implore you to not skip the intuition behind its construction, as is often the case in school maths classes. It would be fairly difficult for students to reverse engineer the logic behind such a compact notation, and it encourages mindless memorisation.

Rising Division

One aspect of polynomial division that always bothered me was the focus on the leading term. It makes sense in the context of integer long division, but polynomials are different in one important aspect: $x$ can be less than $1$, and when it gets very close to $0$, arbitrarily close, if you will, the contribution of the higher powers of $x$ become vanishingly small. In fact, close to $0$, the smaller the power of $x$, the more it contributes to the overall value! The idea behind this vague notion of "close to $0$" is that, for any polynomial, you can pick a value for $x$, at and below which the absolute size of every term strictly decreases as the exponent of $x$ increases. I've assumed positive $x$ for simplicity, but the same idea applies for negative values, and, when you include those, it's more appropriate to think of a range centred on zero instead.

So, how do we divide polynomials under the assumption that the constant term is actually the leading one? Simple, we apply the same rewriting process as regular division, but starting from the low degree terms instead. Let's see how that works for our first polynomial division example, where we divided $f(x): 3x^2 + 2x + 1$ by $g(x): 2x + 9$:

$$\begin{align*} f(x): 3x^2 + 2x + 1 &= 3x^2+2x+1(\frac{1}{9}(9 + 2x - 2x)) \\ &= 3x^2+2x+(\frac{1}{9}(g(x) - 2x)) \\ &= 3x^2 + (2 - \frac{2}{9})x + \frac{1}{9}g(x) \\ &= 3x^2 + \frac{16}{9}(\frac{x}{9}(9 + 2x - 2x)) + \frac{1}{9}g(x) \\ &= 3x^2 + \frac{16}{9}(\frac{x}{9}g(x) - \frac{2}{9}x^2)) + \frac{1}{9}g(x) \\ &= (3-\frac{32}{81})x^2 + \frac{16x}{81}g(x) + \frac{1}{9}g(x) \\ &= \frac{211x^2}{81} + (\frac{16x}{81} + \frac{1}{9})g(x) \end{align*}$$

An easy way to verify this computation is to reverse the list of coefficients passed to our standard division function, then reverse the output again:

def polyDivRise(n: list, d: list) -> tuple[list, list]:
    quotient, remainder = polyDivMultiPass(list(reversed(n)), list(reversed(d)))
    quotient.reverse()
    remainder.reverse()
    return quotient, remainder

printQuotRem(polyDivRise([3, 2, 1], [2, 9]))

Quotient: [ 1/9, 16/81 ] Remainder: [ 211/81 ]

Even though it's inefficient, polyDivRise is good enough for our current needs, and I like the fact that it explicitly demonstrates its connection to the standard algorithm. Is this type of polynomial division actually useful or merely a curiosity? In my opinion, the most important application of rising division can be found in calculus and analysis, because it allows us to divide truncations of two infinite series. I'll show you an example later on!

Naming

When I first stumbled upon this idea I wanted to find the name mathematicians have assigned to it, but my search yielded no results. I eventually stumbled upon "division by increasing power order", mentioned in the documentation for the Giac computer algebra system, which prompted me to search French Wikipedia, where I finally found a proper description of this type of division. The name was too verbose for my taste, so I decided to call it rising polynomial division instead, because we rewrite starting from the lowest degree term; this line of reasoning suggested an alternative name for standard polynomial division - falling polynomial division. Even though it's effectively the same algorithm, I'll use these names from now on to differentiate between the two rewriting directions.

Expansion

Another polynomial division rule that always bothered me was the fact that we had to stop after rewriting the numerator term whose degree matched the degree of the denominator's leading term. The algorithm for computing decimal expansions of rational numbers was taught in primary school, so, by this point, I was quite comfortable with the idea of an infinite process. As a reminder, the decimal expansion of, say, $1/7$ is:

$$\frac{1}{7} = 0.\overline{142857}$$

The line above the fractional digits is called an overline and is used to indicate that those digits will repeat indefinitely. For practical reason, the repeating part is often omitted because even relatively small denominators might require the computation of hundreds of digits before the repeating pattern is established. The algorithm is fairly straightforward, and involves repeated multiplication by $1$, or $10/10$ to be specific:

$$\begin{align*} \frac{1}{7} &= \frac{1 \cdot 10}{7 \cdot 10} = \frac{7 + 3}{7 \cdot 10} = \frac{1}{10} + \frac{3}{7 \cdot 10} \\ &= \frac{1}{10} + \frac{3 \cdot 10}{7 \cdot 10^2} = \frac{1}{10} + \frac{4 \cdot 7 + 2}{7 \cdot 10^2} \\ &= \frac{1}{10} + \frac{4}{10^2} + \frac{2}{7 \cdot 10^2} \end{align*}$$

Computing the first $2$ digits like this is straightforward but cumbersome, which is why schools usually teach a more compact notation that takes advantage of the fact that at each step we multiply the previous remainder by $10$ and divide by the denominator to decompose it into a quotient and remainder. The quotient is output as a digit, and the remainder is carried into the next step:

$$\begin{align*} 1 &= \underline 0 \cdot 7 + \boxed 1 \\ \boxed 1 \cdot 10 &= \underline 1 \cdot 7 + \boxed 3 \\ \boxed 3 \cdot 10 &= \underline 4 \cdot 7 + \boxed 2 \\ \boxed 2 \cdot 10 &= \underline 2 \cdot 7 + \boxed 6 \\ \boxed 6 \cdot 10 &= \underline 8 \cdot 7 + \boxed 4 \\ \boxed 4 \cdot 10 &= \underline 5 \cdot 7 + \boxed 5 \\ \boxed 5 \cdot 10 &= \underline 7 \cdot 7 + \boxed 1 \end{align*}$$

The numbers with boxes around them are the remainders that we carry forward into next line, and the underlined numbers are the output digits of the expansion. You can omit the divisor, or replace it with an easy to write symbol, once you're comfortable with the procedure. In the decimal notation, each digit after the dot is implicitly multiplied by a negative power of $10$ that corresponds to the index of the digit relative to the dot. For example, $0.1428$ is a compact representation of the following:

$$0.1428 = \frac{1}{10} + \frac{4}{10^2} + \frac{2}{10^3} + \frac{8}{10^4}$$

This is very similar to the implicit positive powers of 10 for regular integers. Given how each digit is supposed to be less than $10$, how can we be sure that quotient after each division will produce a valid digit? Fortunately, the proof is quite simple and it relies on the observation that, by definition, each remainder is strictly smaller than the divisor. Suppose the divisor is $D$ and n-th remainder is $R_n$. We can form an equality based on the definition of a remainder $R_n < D$, then divide both sides by $D$, followed by a multiplication by $10$:

$$\begin{align*} R_n &< D \\ \frac{R_n}{D} &< 1 \\ \frac{R_n \cdot 10}{D} &< 10 \end{align*}$$

The left-hand side is how we compute the quotient and remainder, and we can always ensure the previous $R$ is less than the divisor by reducing the numerator with long division first to compute the integer part. This is a rather elementary proof by induction which, sadly, is usually not taught when decimal expansion is introduced. Note how it didn't matter that we multiplied both sides by 10. We could've chosen any other number and the inequality would still hold, which means we can compute expansions in any base! Here's $1/7$ in base $2$:

$$\begin{align*} 1 &= \underline 0 \cdot 7 + \boxed 1 \\ \boxed 1 \cdot 2 &= \underline 0 \cdot 7 + \boxed 2 \\ \boxed 2 \cdot 2 &= \underline 0 \cdot 7 + \boxed 4 \\ \boxed 4 \cdot 2 &= \underline 1 \cdot 7 + \boxed 1 \\ \end{align*}$$

which gives us $0.\overline{001}_2$. The little 2 subscript is the standard notation for specifying the base of the number. Turning this algorithm into code is very easy:

def intDivFalling(n: int, d: int, b: int = 10, termCount: int = 10) -> tuple[int, list]:
    expansion = []

    integerPart = n // d
    n %= d

    for _ in range(termCount):
        n *= b
        expansion.append(n // d)
        n %= d

    return integerPart, expansion

print(intDivFalling(1, 13))

(0, [0, 7, 6, 9, 2, 3, 0, 7, 6, 9])

Given how we can expand in any base $x$ through multiplication by $x/x$, why not do the same with polynomial falling division?

Rising and Falling Expansions

Let's apply our rational number expansion intuition to the falling division of $f(x): x$ by $g(x): x^2 - x - 1$:

$$\begin{align*} f(x): x &= \frac{x^2}{x} = \frac{(x^2 - x - 1) + x + 1}{x}=\frac{g(x) + x + 1}{x} \\ &= \frac{g(x)}{x} + \frac{(x+1)x}{x^2} = \frac{g(x)}{x} + \frac{x^2 + x}{x^2} \\ &= \frac{g(x)}{x} + \frac{g(x) + 2x + 1}{x^2} = \frac{g(x)}{x} + \frac{g(x)}{x^2} + \frac{2x+1}{x^2} \\ &= \frac{g(x)}{x} + \frac{g(x)}{x^2} + \frac{2x^2+x}{x^3} = \frac{g(x)}{x} + \frac{g(x)}{x^2} + \frac{2 \cdot g(x)}{x^3}+ \frac{3x+2}{x^3} \end{align*}$$

You could continue doing this by hand, but, if you're lazy like me, you might already be thinking of ways to reuse our polynomial division code. The good news is that you don't even need to touch the code! Repeatedly multiplying by $x/x$ and dividing is equivalent to multiplying the numerator by a power of $x$ that corresponds to the number of terms you want to compute, and then applying the standard division algorithm. All you have to do at the end is divide the result by the same power of $x$ to preserve the equality:

$$\begin{align*} f(x): x &= \frac{x^2 \cdot x^3}{x^4} = \frac{(g(x)+x+1)x^3}{x^4} = \frac{(g(x)x+x^2+x)x^2}{x^4} \\ &=\frac{(g(x)x + g(x)+2x+1)x^2}{x^4} = \frac{(g(x)x^2 + g(x)x+2x^2+x)x}{x^4} \\ &=\frac{(g(x)x^2 + g(x)x+2g(x) + 3x + 2)x}{x^4} \\ &=\frac{g(x)x^3 + g(x)x^2 +2g(x)x + 3x^2 + 2x}{x^4} \\ &=\frac{g(x)x^3 + g(x)x^2 +2g(x)x + 3g(x) + 5x+3}{x^4} \\ &=g(x)(\frac{1}{x}+\frac{1}{x^2} + \frac{2}{x^3} + \frac{3}{x^4}) + \frac{5x + 3}{x^4} \end{align*}$$

Multiplying by powers of $x$ is equivalent to prepending zeros in the coefficient list that represents polynomials:

n = [0] * 10 + [0, 1] # compute 10 terms by prepending 10 zeros, equal to a multiplication of x^10
d = [-1, -1, 1] # x^2 - x - 1
print(polyDiv(n, d))

[55, 34, 21, 13, 8, 5, 3, 2, 1, 1]

Applying the division by $x^{10}$ to undo the previous multiplication we get:

$$\frac{1}{x}+\frac{1}{x^2}+\frac{2}{x^3}+\frac{3}{x^4}+\frac{5}{x^5}+\frac{8}{x^6}+\frac{13}{x^7}+\frac{21}{x^8}+\frac{34}{x^9}+\frac{55}{x^{10}}$$

You might recognise the numerators as the Fibonacci numbers, and it shouldn't surprise you that they came out given the previous discussion about linear recurrences. The division turned into a pure linear recurrence because the leading divisor coefficient is $1$, and the lower degree coefficients were filled with zeros, so there was nothing to subtract from when computing the quotient coefficients.

Rising division is even simpler because we don't need to worry about adjusting the powers of $x$ at the end. Let's apply it to a special case of the geometric series formula:

$$\frac{1}{x-1}$$

This time we need to append zeros instead because the list will be reversed:

num = [1] + [0]*10
div = [1, -1]
print(polyDivRise(num, div))

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

Which corresponds to:

$$x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1$$

I promised to give you a calculus example earlier, so let's take a few terms from the $sin(x)$ and $cos(x)$ series expansions that we've seen before, and use them to compute the first few terms of the $tan(x)$ series:

$$tan(x)=\frac{sin(x)}{cos(x)} \approx\frac{x -\frac{x^3}{6}+\frac{x^5}{120}-\frac{x^7}{5040}}{1-\frac{x^2}{2}+\frac{x^4}{24}-\frac{x^6}{720}}$$

from fractions import Fraction as Q

num = [0, Q(1), 0, Q(-1, 6), 0, Q(1, 120), 0, Q(-1, 5040)] + [0]*6 # sin(x) series
div = [Q(1), 0, Q(-1, 2), 0, Q(1, 24), 0, Q(-1, 720)] # cos(x) series
printQuotRem(polyDivRise(num, div))

Quotient: [ 0, 1, 0, 1/3, 0, 2/15, 0, 17/315 ] Remainder: [ 0, 331/15120, 0, -13/6300, 0, 17/226800 ]

The quotient gives us the correct expansion up to the 7th power:

$$tan(x) \approx\ x + \frac{x^3}{3} + \frac{2 x^5}{15}+ \frac{17 x^7}{315}$$

The general formula for the tangent series coefficients is quite complicated, because it involves the Bernoulli numbers, and yet we can compute as many terms as we want with a simple algorithm! You can find more efficient formulas for the coefficients of various series by carefully studying the behaviour of rising division, and the patterns in each iteration of its computation. You can also use both rising and falling expansions to approximate the values of tricky integrals.

If you're familiar with calculus you might've realised that our algorithms are producing an infinite series approximation of a rational function. You might also suspect that the series will diverge around the zeros (poles) of the divisor, which, combined with a simple divergence testing method like the ratio test, can give give us one of the oldest polynomial root-finding algorithms - Bernoulli's Method. Optimisations based on this idea underlie many modern root-finding algorithms, which will be the subject of a future article.

We've only scratched the surface of the beautiful structure that underlies these simple algorithms, so feel free to play around with the code, evaluate it symbolically in terms of the roots of the divisor, and look for patterns.

Future Explorations

In future articles, we'll explore how expressing linear recurrences as matrix operations will allow us to create an important bridge between linear algebra and polynomials through the Companion matrix, whose eigenvalues can be used to approximate the roots of a polynomial. We'll also see how any term in a linear recurrence can be directly computed through a formula that depends on the roots of the polynomial division that corresponds to it. This also provides a natural connection between polynomials and a certain class of pseudo-random number generators.

Reversing polynomials expansion will turn out to be important for the efficient implementation of certain error correcting codes, the inversion of Toeplitz matrices, and even rational approximation.

Polynomial division leads to an algorithm for solving systems of polynomial equations, which, in turn, provides the tools necessary for a deep exploration of the theory of algebraic numbers with minimal background knowledge requirements. It turns out that the same ideas can be applied to the problem of partial fraction decomposition.

Division also turns out to be important in explaining what it means to change the basis of a polynomial, and how that relates to various exotic number bases. The subjects of interpolation and approximation theory will also make an appearance during this investigation.

Integers and rational numbers will closely accompany us at each step along the way, providing context and inspiration, particularly when it comes to data compression. This reminds me, before I wrap up, there's one last key tool missing from our division algorithms toolkit.

Rising Integer Division

After playing around with polynomial expansions for a while, I started wondering if there's an integer equivalent for rising division. Falling polynomial division was inspired by the classic integer expansion algorithm I learned in school, so maybe I can go in the opposite direction, and adapt the polynomial rewriting approach for integers! Let's try this idea with our old friend, $1/7$. Can I rewrite the $1$ as a sum of a digit, which is less than $10$, multiplied by $7$, and some remainder? We have a lot of choices if the remainder is negative, but given the fact that we'll continuously rewrite this remainder, it better be a multiple of the base, which will be $10$ in this example, so we can get an increasingly larger base exponent. These requirements correspond to the following symbolic constraints:

$$\begin{cases} 7 x + 10 y = 1 \\ 0 \le x < 10 \end{cases}$$

Is it even possible to find solutions to the equation? If $7$ and $10$ had a common factor $m$, greater than $1$, we could rewrite the equation as $m (a x + b y) = 1$, for some $a$ and $b$, which obviously has no integer solutions because $m > 1$. I picked $7$ and $10$ so they're relatively prime (coprime), meaning their greatest common factor is 1, but if we were trying to find an expansion for, say, $1/15$ in base $10$, we'd fail, because $15$ and $10$ have $5$ as a common factor. This is a rather serious limitation compared to the usual expansion, though we don't need to worry about it if we pick a prime number as the base, and factor out all multiples of the base from the denominator before we begin the expansion.

Putting this limitation aside, how do we find appropriate values for $x$ and $y$? The brute-force approach here would be to simply try out increasingly larger values for $x$ until the product has a lowest digit equal to $1$. It only takes $2$ attempts to arrive at $7 \cdot 3 = 21$, which means we can subtract $2 \cdot 10$ to get 1, giving us $x = 3$ and $y = -2$

A more elegant approach relies on the observation that we can rewrite $10$ as $7 + 3$, which allows us to reduce the size of the coefficients we're considering: $$1 = 7 x + (7 + 3) y = 7 (x + y ) + 3 y$$

We can repeat this process until one of the outer coefficients becomes $1$:

$$1 = (3 \cdot 2 + 1) (x + y ) + 3 y = (x + y) + 3(y + 2 (x + y))$$

One of the multipliers became 1, so at this point we have two obvious options for its corresponding expression: $- 2 + 3 \cdot 1$ or $4 - 3 \cdot 1$

The first gives us:

$$\begin{cases} -2 &= x + y \\ 1 &= y + 2 (x + y) = y + 2 (-2) \end{cases}$$

which we can solve through simple substitution to find $y = 5$ and $x=-7$. You can check for yourself that the second option results in $x = 13$ and $y = -9$. Neither of these satisfies the $0 \le x < 10$ inequality, but we can easily fix that by rewriting $-7$ as $3 - 10$; expand and re-collect the terms:

$$\begin{align*} (3 - 10) \cdot 7 + 5 \cdot 10 &= 3 \cdot 7 + (5 - 7) \cdot 10 \\ &= \boxed{3 \cdot 7 + (-2) \cdot 10 = 1} \end{align*}$$

Depending on your maths background, you might've recognised the rewriting trick as an application of the Extended Euclidean algorithm (EEA), and the equation $7 x + 10 y = 1$ as a special case of Bézout's identity. The EEA will be our constant companion throughout the upcoming series of articles, and we'll explore it in more detail in part 2!

With our $x$ and $y$ coefficients in hand, we can now rewrite $1/7$ as:

$$\frac{(3 \cdot 7 - 2 \cdot 10)}{7} = 3 + \frac{(-2 \cdot 10)}{7}$$

To proceed from here we need to repeat the same process for $-2$, which is similar to our remainder term in the regular expansion algorithm. Because we already have a way to express $1$, let's just multiply that equation on both sides by $-2$:

$$-2 = (-6) \cdot 7 + 4 \cdot 10$$

Negative digits aren't allowed in the usual decimal expansion, so let's express $-6$ as $4 - 10$, which allows us to combine the $-10$ term with the $4 \cdot 10$ one:

$$\begin{align*} -2 &= (4 - 10) 7 + 4 \cdot 10 \\ &= 4 \cdot 7 - 10 \cdot 7 + 4 \cdot 10 \\ &= 4 \cdot 7 + (4 - 7) 10 \\ &= 4 \cdot 7 + (-3) 10 \end{align*}$$

We've finally arrived at an appropriate expression for $-2$, which we can substitute back into our expansion:

$$\begin{align*} 3 + \frac{-2 \cdot 10}{7} &= 3 + \frac{(4 \cdot 7 + (-3)10) 10}{7} \\ &= 3 + 4 \cdot 10 + \frac{-3 \cdot 10^2}{7} \end{align*}$$

From here on we just repeat the same process, but with $-3$ instead:

$$\begin{align*} -3 &= -9 \cdot 7 + 6 \cdot 10 \\ &= (1 - 10) 7 + 6 \cdot 10 \\ &= 1 \cdot 7 + (6 - 7) 10 \\ &= 1 \cdot 7 + (-1) 10 \end{align*}$$

Rewrite the rational term of the expansion once again:

$$\begin{align*} \frac{-3 \cdot 10^2}{7} &= \frac{(1 \cdot 7 + (-1) 10) 10^2}{7} \\ &= 1 \cdot 10^2 + \frac{-1 \cdot 10^3}{7} \end{align*}$$

Our next digit is $1$ and rewriting $-1$ is equally straightforward:

$$\begin{align*} -1 &= -3 \cdot 7 + 2 \cdot 10 \\ &= (7 - 10) 7 + 2 \cdot 10 \\ &= 7 \cdot 7 + (2 - 7) 10 \\ &= 7 \cdot 7 + (-5) 10 \end{align*}$$

At this point you might be starting to see the pattern. The power of $10$ increases by $1$ and the corresponding coefficient (digit) is $7$. We can express $-5 \cdot 3$ as $5 - 2 \cdot 10$, which obviously makes the digit after it $5$. while the next remainder is $5 \cdot 2 - 2 \cdot 7 = -4$. Our expansion so far looks like this:

$$3 + 4 \cdot 10 + 1 \cdot 10^2 + 7 \cdot 10^3 + 5 \cdot 10^4 + (-4) 10^5 / 7$$

We can follow a similar logic for $-4 \cdot 3 = (8 - 2*10)$, giving us digit $8$ and remainder $4 \cdot 2 - 2 \cdot 7 = -6$. Repeat one more time: $-6 \cdot 3 = 2 - 2 \cdot 10$, which gives us digit $2$ and remainder $6 \cdot 2 - 2 \cdot 7 = -2$. Perhaps unsurprisingly, we eventually end up with a remainder that we've seen before, and in this case it goes all the way back when we computed the digit $4$. This means that the digits $4, 1, 7, 5, 8, 2$ will continue repeating. Here's what we've computed so far: $$3 + 4 \cdot 10 + 1 \cdot 10^2 + 7 \cdot 10^3 + 5 \cdot 10^4 + 8 \cdot 10^5 + 2 \cdot 10^6 + (-2)10^7 / 7$$

This looks an awful lot like a regular integer, except the powers continue to increase, and the repeating part is similar to repeating expansions. Writing out the powers of 10 is becoming cumbersome, so let's repurpose regular integer notation, and combine it with the overline used for repeated digits in regular expansions:

$$\overline{285714}3$$

Compare this against the usual base $10$ expansion we computed previously:

$$0.\overline{142857}$$

Do you see the pattern? Starting from the 2nd digit and going backwards we get 4, and 1, then we loop to the end and continue: $7, 5, 8, 2$. This a cyclic permutation (rotation) of the repeated digits.

Compact Notation

It would be nice to have a compact notation similar to the one we used for falling integer expansions:

$$\begin{align*} 1 &= \underline{3} \cdot 7 + \boxed{-2} 10 \\ \boxed{-2} &= \underline{4} \cdot 7 + \boxed{-3} 10 \\ \boxed{-3} &= \underline{1} \cdot 7 + \boxed{-1} 10 \\ \boxed{-1} &= \underline{7} \cdot 7 + \boxed{-5} 10 \\ \boxed{-5} &= \underline{5} \cdot 7 + \boxed{-4} 10 \\ \boxed{-4} &= \underline{8} \cdot 7 + \boxed{-6} 10 \\ \boxed{-6} &= \underline{2} \cdot 7 + \boxed{-2} 10 \end{align*}$$

The equality in each row is easy to verify, but the downside is that you have to perform all the important calculations somewhere else. If the divisor is very large, you can substitute any easy to write symbol for it. Also, feel free to omit the $10$ too, if you feel comfortable with the calculations. An alternative, more computation-friendly, notation can look something like this:

$$\begin{align*} 1 &= \underline{3} \cdot 7 + 10 \cdot \underline{\overline{-2}} \\ \boxed{-6} &= \underline{4} - 1 \cdot 10 \rightarrow (\underline{\overline{-2}} - \underline{4} \cdot 7) / 10 = \underline{\overline{-3}} \\ \boxed{-9} &= \underline{1} - 1 \cdot 10 \rightarrow (\underline{\overline{-3}} - \underline{1} \cdot 7) / 10 = \underline{\overline{-1}} \\ \boxed{-3} &= \underline{7} - 1 \cdot 10 \rightarrow (\underline{\overline{-1}}- \underline{7} \cdot 7) / 10 = \underline{\overline{-5}} \\ \boxed{-15} &= \underline{5} - 2 \cdot 10 \rightarrow (\underline{\overline{-5}} - \underline{5} \cdot 7) / 10 = \underline{\overline{-4}} \\ \boxed{-12} &= \underline{8} - 2 \cdot 10 \rightarrow (\underline{\overline{-4}} - \underline{8} \cdot 7) / 10 = \underline{\overline{-6}} \\ \boxed{-18} &= \underline{2} - 2 \cdot 10 \rightarrow (\underline{\overline{-6}} - \underline{2} \cdot 7) / 10 = \underline{\overline{-2}} \end{align*}$$