Essays/Polynomial+Division

Ken Iverson's APL Style essay talks briefly about polynomial division. This "essay" is mostly a translation of his treatment of polynomial division into J.

Here, instead of using J's +//.@(*/) to multiply polynomials, he develops a different approach:

   t=: [ +/ .*  b +/ .* ]
   b=: ([ =/ ] -/~ _1 i.@+ +&#)&i.&#

Here, t multiplies polynomials, and b creates an element selection array.

And, he suggests this example as illustrative:

   ]e=.(c=.1 2 1) t (d=.1 3 3 1)
1 5 10 10 5 1

(Here, e would represent the polynomial (1) + (5*y) + (10*y^2) + (10*y^3) + (5^y^4) + (1*y^5). Similarly 0 0 1 would represent the polynomial y^2 (or, if you prefer: (0) + (0*y) + (1*y^2)).)

Anyways, he goes on to implement polynomial division, based on the idea that {{ x +/ .* x b y }} builds a matrix.

First, he expresses bq which builds that element selection array array from c and e (instead of c and d -- d would be the unknown that we are trying to find):

      bq=: {{ (i.#x) =/ (i.#y) -/ i.1+ y-&# x }}

Or, equivalently:

   bq=: ([ =/ ] -/ 1 i.@+ -&#~)&i.&#

And, here, he suggests that we can think of inverting the matrix c +/ .* c bq e as representing polynomial division.

But, for the general case, we want to invert a square matrix. So, he proposes two different approaches:

   dbho=: {{ (d{.x)%.(2#d=. <./$M){.M=.y +/ .*y bq x }}
   dblo=: {{ (d{.x)%.(2#d=.-<./$M){.M=.y +/ .*y bq x }}

Here, dbho yields a high order remainder and dblo yields a low order remainder. To illustrate, he proposes these examples:

   E=. 1 5 10 10 7 4
   C=. 1 2 1
   ]Q=. E dbho C
1 3 3 1
   ]R=. E-C t Q
0 0 0 0 2 3

And, for a low order remainder:

   e=.4 7 10 10 5 1
   c=. 1 2 1
   ]q=. e dblo c
1 3 3 1
   ]r=. e-c t q
3 2 0 0 0 0

---

Another approach to polynomial division uses the "long division" or "synthetic division"

pDiv=: {{
  q=. $j=. 2 + x -&# y
  'x y'=. x,:y
  while. j=. j-1 do.
    q=. q, r=. x %&{. y
    x=. 1 |.!.0 x - y*r
  end.q
}}

   1 5 10 10 5 1 pDiv 1 2 1
1 3 3 1
   1 5 10 10 5 1 pDiv 1 3 3 1
1 2 1
   1 5 10 10 7 4 pDiv 1 3 3 1
1 2 1

---

The FFT essay indicates that polynomial multiplication corresponds with J's * on the fourier transformation of the polynomial. Thus, if there's no remainder to contend with, polynomial division corresponds with J's % on the fourier transformation of the polynomial.

In other words:

   %&.(fft :.ifft)/1 5 10 10 5 1,:1 2 1
1 3 3 1 0 0
   %&.(fft :.ifft)/1 5 10 10 5 1,:1 3 3 1
1 2 1 7.40149e_17 0 1.4803e_16

(Some minor cleanup may be necessary on these results... but when working with large polynomials, the speedup from this approach is significant.) [[Category:Math R.3.2]