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
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 in which is of lower order than , in a form such as
where are the roots of . 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 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 from say , 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 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 to the function .
Using this representation of polynomials, a fraction such as is a 2-list of polynomials, for example 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 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 :
]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 is not unique – it could be the “obvious” factorisation of or it could be , or 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 or equivalently .
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
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 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 where the appropriate partial fraction structure is . 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, .
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 , 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
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 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
- .
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. )