# Fifty Shades of J/Chapter 45

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

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

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

- k1/(a1-x) + k2/(a2-x) + …

where a1 ,a2,.. 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 (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 (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 {2/(4+x)} – {6/(1+5x)} is

t0 2 0 4 1 _6 0 1 5

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

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 (11x + 8)/(3 - 2x - 8x^{2}) :

u1 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 nth. 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 {-4/(-6 - 8x)} – {1.5/(-0.5 + x)} or equivalently {2/(3 + 4x)} + {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

- (23+55x+8x
^{2})/(3+13x-18x^{2}-40x^{3})

to resolve, define:

u2 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 `u0` within constant factors :

pf cp/t0 10 0 20 5 _1.2 0 0.2 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.

For complex roots, take as an example 1/(1+x^{5}) represented by

u3=.2 6$1 0 0 0 0 0 1 0 0 0 0 1 u3 1 0 0 0 0 0 1 0 0 0 0 1

pf operates as for real roots

pf u3 _0.161803j_0.117557 0 _0.809017j_0.587785 1 _0.161803j0.117557 0 _0.809017j0.587785 1 0.2 0 1 1 0.0618034j_0.190211 0 0.309017j_0.951057 1 0.0618034j0.190211 0 0.309017j0.951057 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 :

cp/2{.pf u3 NB. a quadratic partial fraction 0.4 _0.323607 0 1 _1.61803 1 cp/_2{.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.

Next, multiple roots, say the resolution of (6+8x+3x^{2})/(1+x)^{3} where the appropriate partial fraction structure is k_{1}/(1+x)^{3} + k_{2}/(1+x)^{2} + 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, (6+8x+3x^{2})/(1+x)^{3} = 1/(1+x)^{3} +2(1+x)^{2} +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 (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

- (4 + 18x + 9x
^{2})/(2 + 3x)^{3}= {-4/(2 + 3x)^{3}} + {2/(2 + 3x)^{2}} + {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 (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 m1 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

- {4/(2 + 3x)
^{3}} - {2/(2 + 3x)^{2}} - {1/(2 + 3x)} + {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 rtof=: ,.&1@-@(1&of) NB. convert roots to factors .. facs=: (0&of)mhead rtof NB. .. and adjust for multiplier of=: >@{pfd NB. pick from denom 0=mult, 1=roots 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 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. )