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]