Fifty Shades of J/Chapter 45

Table of Contents ... Glossary ... Previous Chapter ... Next Chapter

A partial solution to a partial problem

Principal Topics

p. (polynomial), %. (matrix inverse), /. (oblique), ~ (reflex), -. (less), polynomial quotients, polynomial multiplication, polynomial factors, roots of equations

Introduction

A seemingly innocuous post on the J Forum asked if anyone had a general routine for resolving partial fractions. Given that the heavy power-horses of p. and %. are already present in J, it seemed to require just small extensions of these to solve the problem, and it may be that this is indeed the case. Nevertheless I found this to be quite a tricky exercise, and as the title above suggests, the path to a general partial fraction algorithm as given here is not quite complete. However, some of the J features which turn up on the way are interesting in their own right, notably the use of explicit rank.

First, the basic problem is relatively simple, namely that of rewriting a quotient of polynomials such as $f(x)/g(x)$ in which $f(x)$ is of lower order than $g(x)$, in a form such as

${\frac {k1}{a_{1}-x}}+{\frac {k_{2}}{a_{2}-x}}+\ldots$

where $a_{1},a_{2},\ldots$ are the roots of $g(x)$. Initially, assume that these are distinct.

Start by defining a polynomial as a list of coefficients in ascending power order which is the meaning implicit in the right argument of p., so that the algebraic function $f(x)=2x^{2}+5x-3$ is represented by the polynomial _3 5 2. The value returned by monadic p. is

   p. _3 5 2
┌─┬──────┐
│2│_3 0.5│
└─┴──────┘

that is the roots _3 0.5 along with the highest order coefficient 2 which helps distinguish $f(x)$ from say $g(x)=6x^{2}+15x-9$, which has the same roots.

To model partial fractions, an initial decision has to be made between using boxes or lists. For example, the partial fraction ${\frac {1+x}{1+3x+2x^{2}}}$ could be represented either by a list of boxes

1 1;1 3 2

or by a list of polynomials

1 1 0
1 3 2

My general principle is to use lists wherever possible unless data-inherent heterogeneity makes lists too burdensome. Since the contents of boxes are sealed by definition, the only available operations without opening are the relatively simple ones of joining and shaping. With lists, operations specifically appropriate to the data types are fully available, subject to the constraints of structural rectangularity which may force the insertion of fill characters. These can sometimes be benign, as in the case of polynomials where adding a couple of zeros on the right of, say, the polynomial 1 0 2 merely adds two non-contributory terms $0x^{3}+0x^{4}$ to the function $1+2x^{2}$.

Using this representation of polynomials, a fraction such as $f(x)/g(x)$ is a 2-list of polynomials, for example ${\frac {4-x}{1+2x^{2}}}$ is the 2-list

4 _1 0
1 0 2

Further, the model is readily extendible by representing a number of such polynomial fractions as a list of 2-lists of lists, for example ${\frac {2}{4+x}}-{\frac {6}{1+5x}}$ is

    ]t0=:  >"1 (2;4 1) ,: (_6;1 5)
 2   0
 4   1

_6   0 
 1   5

Basic partial fractions to composite partial fraction

Next, the dyadic verb pmult (polynomial multiply) multiplies two polynomials and is often cited as an illustration of the adverb oblique /.

   pmult=:+//.@(*/)     NB. (dyad) polynomial multiplication
   4 _1 pmult 1 0 2
4 _1 8 _2

(Assume in what follows that defined verbs are monadic unless there is a specific comment to the contrary.)

A useful first step in developing a partial fraction algorithm is to develop an inverse verb, that is one which combines basic (that is 2 by 2) partial fractions into a single composite partial fraction.

   cp=:(+/@(pmult"1|.)),:pmult&(1&{)    NB. combine basic p.frac’ns
   cp/t0
_22  4 0
  4 21 5

Distinct real roots

The roots and multiplier of the denominator of a partial fraction are obtained by

   pfd=:p.@(1&{)      NB. partial fraction denominator
   of=:>@{pfd         NB. (dyad) pick from p. 0=multiplier, 1=roots

Thus if u1 represents the partial fraction ${\frac {11x+8}{3-2x-8x^{2}}}$ :

   ]u1=:11 8 0,:3 _2 _8
11  8  0
 3 _2 _8
   0 of u1        NB. multiplier
_8
   1 of u1        NB. roots of denominator
_0.75 0.5

Although having p. return a multiplier (which is already part of the input data) as well as roots adds complexity to its result, this approach is well justified by considering that the factorisation of an expression such as $3-5x-2x^{2}$ is not unique – it could be the “obvious” factorisation of $(1-2x)(3+x)$ or it could be $(0.5-x)(6+2x)$, or $(2-4x)(1.5+0.5x)$ and so on. Thus in factorising a polynomial into linear factors (which is always possible because every n^th order polynomial has exactly n roots, allowing for possible complex values), the multiplier delivered by p. must be applied arbitrarily to one of the factors. The choice made here is to apply it to the first. A general verb for multiplying the head only of a list is most clearly expressed as an explicit function:

   mhead=:dyad :'(x*{.y),}.y'    NB. x*head of y,tail unchanged

following which roots are changed into factors by rtof, which reverses signs and catenates 1s, and adjusted by the multiplier with facs:

   rtof=:,.&1@-@(1&of)         NB. turn roots into factors ...
   facs=:(0&of)mhead rtof      NB. ... and adjust for multiplier

If there are just two roots, everything is in place to find the factors:

      facs u1                  NB. factors
  _6 _8
_0.5  1

... and then the partial fraction coefficients:

   11 8 %.|:facs u1            NB. coefficients
_1.5 _4

so that the resolution is ${\frac {-4}{-6-8x}}-{\frac {1.5}{-0.5+x}}$ or equivalently ${\frac {2}{3+4x}}+{\frac {3}{1-2x}}$.

More generally, the numerator coefficients are equated to the combined coefficients of the denominator polynomial factors following multiplication with one factor at a time omitted. This gives two nice examples of the use of explicit rank, one to exclude each item in turn in the list of factors, then pmult is used at rank 2 to multiply these factors together. For example, given the fraction

${\frac {23+55x+8x^{2}}{3+13x-18x^{2}-40x^{3}}}$

to resolve, define:

   ]u2=:>23 55 8;3 13 _18 _40    NB. partial fraction
23 55   8   0
 3 13 _18 _40

   facs u2            NB. basic polynomial factors of u2
 _30 _40
_0.5   1
 0.2   1

Now use less which for x-.y includes all items of x except those which are cells of y:

   minors=:-."2 1~
   minors facs u2
_0.5   1
 0.2   1

 _30 _40
 0.2   1

 _30 _40
_0.5   1

and ‘polynomial multiply’ within each pair of factors

   mat=:pmult/"2@minors@facs NB. lin eqns for pf coeffs
   mat u2
_0.1 _0.3   1
  _6  _38 _40
  15  _10 _40

At this point the power of matrix divide comes into play and is consolidated in a verb:

   pfc=:}:@(0&{) %.|:@mat    NB. partial fraction coefficients
   pfc u2
_20 _1.5 0.8

which gives the partial fraction coefficients which correspond to facs u2 in order.

Finally the verb form brings coefficients and factors together in the basic partial fraction representation:

   form=:(,:~,)"1 0     NB. (dyad) merge factors and coefficients
   pf=:facs form pfc    NB. partial fractions
   pf u2
 _20   0
 _30 _40

_1.5   0
_0.5   1

 0.8   0
 0.2   1

It is easy to confirm that cp/pf u2 is identical to u2 and similarly for u1. As further confirmation using the first example pf cp/t0 is the same as t0 within constant factors:

   pf cp/t0
  10 0
  20 5

_1.2 0
 0.2 1
   (5 0.2 * t0) -: pf cp/t0
1

The problem of distinct real denominator roots has thus been fully dealt with, which leaves two matters outstanding, namely complex roots, and repeated roots.

Complex roots

For complex roots, take as an example ${\frac {1}{1+x^{5}}}$ represented by

   ]u3=:2 6$1 0 0 0 0 0 1 0 0 0 0 1
1 0 0 0 0 0
1 0 0 0 0 1

pf operates as for real roots

cleancomp =: (**@|)&.+. NB. remove complex rounding errors
   cleancomp pf u3
_0.161803j_0.117557 0
_0.809017j_0.587785 1

 _0.161803j0.117557 0
 _0.809017j0.587785 1

0.0618034j_0.190211 0
 0.309017j_0.951057 1

 0.0618034j0.190211 0
  0.309017j0.951057 1

                0.2 0
                  1 1

If the implementation dependent assumption is made that imaginary complex pairs occur together, there is no problem in combining basic fractions into partial fractions with real coefficients as in :

   cleancomp cp/ 0 1{pf u3          NB. a quadratic partial fraction
0.4 _0.323607 0
  1  _1.61803 1
   cleancomp cp/ 2 3{pf u3         NB. another quadratic partial fraction
0.4 0.123607 0
  1 0.618034 1

It is a straightforward matter to write a verb which post-processes the results of pf by combining each complex fraction with its neighbour.

Multiple roots

Next, multiple roots, say the resolution of ${\frac {6+8x+3x^{2}}{(1+x)^{3}}}$ where the appropriate partial fraction structure is ${\frac {k_{1}}{(1+x)^{3}}}+{\frac {k_{2}}{(1+x)^{2}}}+{\frac {k_{3}}{(1+x)}}$. This is a relatively easy application of the binomial coefficients and matrix divide:

   bc=.!/~@i.           NB. binomial coefficients
   bc 3
1 1 1
0 1 2
0 0 1

and the relevant coefficients are found by

   6 8 3 %.bc 3
1 2 3

that is, ${\frac {6+8x+3x^{2}}{(1+x)^{3}}}={\frac {1}{(1+x)^{3}}}+{\frac {2}{(1+x)^{2}}}+{\frac {3}{(1+x)}}$.

More generally the binomial coefficients for the polynomial a b are found by

   bco=.dyad :0        NB. (dyad) coeffs. of (polynomial y)^x
r=.(,1),:y [ i=.2
while. i<x do.
  r=.r,({:r) pmult y [ i=.i+1
end.
)

so for a partial fraction ${\frac {4+18x+9x^{2}}{(2+3x)^{3}}}$, the first step is the matrix

   3 bco 2 3
1  0 0
2  3 0
4 12 9

followed by matrix divide to obtain the coefficients :

   4 18 9%.|:3 bco 2 3
_4 2 1

that is

${\frac {4+18x+9x^{2}}{(2+3x)^{3}}}={\frac {-4}{(2+3x)^{3}}}+{\frac {2}{(2+3x)^{2}}}+{\frac {1}{(2+3x)}}$

Notice that bco depends only on pmult and not on any of the factorisation verbs. If the denominator is multiplied by a further factor, say resolve ${\frac {4+14x+27x^{2}+18x^{3}}{(x+1)(2+3x)^{3}}}$ into partial fractions, each of the polynomials listed in 3 bco 2 3 must be multiplied by the new factor (hence pmult"1), and a third order set of binomial coefficients added:

   ]m1=:((3 bco 2 3)pmult"1(1 1)),{:4 bco 2 3
1  1  0  0
2  5  3  0
4 16 21  9
8 36 54 27

   4 14 27 18 %.|:m1
4 _2 _1 1

that is the resolution is

${\frac {4}{(2+3x)^{3}}}-{\frac {2}{(2+3x)^{2}}}-{\frac {1}{(2+3x)}}+{\frac {1}{1+x}}$ .

At the start I said that the journey was not complete, but at least a staging post has been reached from which it is just a matter of conscientious programming to achieve a general partial fraction algorithm.

Code Summary

pmult =: +//.@(*/)                     NB. (dyad) polynomial multiplication
cp    =: (+/@(pmult"1|.)),:pmult&(1&{) NB. combine partial fractions
pfd   =: p.@(1&{)                      NB. partial frtn denominator
of    =: >@{pfd                        NB. pick from denom 0=mult, 1=roots
rtof  =: ,.&1@-@(1&of)                 NB. convert roots to factors ..
facs  =: (0&of)mhead rtof              NB. .. and adjust for multiplier
mhead =: dyad :'(x*{.y),}.y'           NB. x*head of y,tail unchanged
minors=: -."2 1~                       NB. uses dyadic -. (less)
mat   =: pmult/"2@minors@facs          NB. lin. eqns for pfractn. coeffs
pfc   =: }:@(0&{) %. |:@mat            NB. partial fraction coefficients
form  =: (,:~,)"1 0                    NB. merge factors and coefficients
pf    =: facs form pfc                 NB. construct partial fraction
cleancomp=: (**@|)&.+.                 NB. remove complex rounding errors
bc    =: !/~@i.                        NB. binomial coefficients

bco   =: dyad :0                       NB. (dyad) coeffs. of (polynomial y)^x
r=. (,1),:y [ i=. 2
while. i<x do.
  r=. r,({:r)pmult y [ i=. i+1
end.
)

Script

File:Fsojc45.ijs