Recent posts on quadratics and root finding, especially in a
demonstrational manner, provoked me to post the following. It's
rather long and demonstrates techniques of another programming
language – a smooth mix of symbolic and numeric computations. The
reason I am posting is that it seems rather useful to me but I don't
believe to have seen this kind of treatment of the subject in the J
literature or the forums. Perhaps someone will be provoked, in turn,
to search for a way of doing something similar in J.
Define a quadratic polynomial:
q0 = 2*x^2+(-3)*x+1;
It can be simplified using the Horner's rule, ax²+bx+c ≡ (ax+b)x+c:
q1 = reduce q0 with a*x^2+b*x+c = (a*x+b)*x+c end;
It is worth defining the rule separately, as a transforming function
applicable to any quadratic polynomial and not just q0:
qhorner = reduce with a*x^2+b*x+c = (a*x+b)*x+c end;
Now the definition of q1 is:
q1 = qhorner q0;
Let's check what we have so far (comments after // show the output):
q0; // 2*x^2+(-3)*x+1
q1; // (2*x+-3)*x+1
(Both these are symbolic expressions with x a variable.)
As promised, qhorner can be applied to other quadratics, including
ones whose coefficients are expressions:
q2 = qhorner ((R+S)*x^2+(T+4)*x+W);
q2; // ((R+S)*x+(T+4))*x+W
Now let's define a rule that substitutes whatever we would wish
for the 'x' in an expression (quadratic or other):
subst e t = reduce e with x = t end;
subst can be used to evaluate an expression to a number, but that
is only one of its uses. Check how it works:
subst q0 5; // 36.0
subst q1 5; // 36
subst q2 5; // ((R+S)*5+(T+4))*5+W
Expressions also can be substituted for x:
q01 = subst q0 (5*u+7);
q11 = subst q1 (5*u+7);
q21 = subst q2 (5*u+7);
Checking these definitions:
q01; // 2*(5*u+7)^2+(-3)*(5*u+7)+1
q11; // (2*(5*u+7)+-3)*(5*u+7)+1
q21; // ((R+S)*(5*u+7)+(T+4))*(5*u+7)+W
Any of these can be simplified (or the opposite) by substituting for
the 'u'. Here is how this works with u = -2:
reduce q01 with u = -2 end; // 28.0
reduce q11 with u = -2 end; // 28
reduce q21 with u = -2 end; // ((R+S)*(-3)+(T+4))*(-3)+W
And now let's devise a new transforming function, which finds the
roots of a quadratic:
roots = reduce with a*x^2+b*x+c = a*(x-x1)*(x-x2)
when d = sqrt (b^2-4*a*c);
x1 = (-b-d)/(2*a);
x2 = (-b+d)/(2*a) end end;
Using roots we obtain:
roots q0; // 2*(x-0.5)*(x-1.0)
Again, we can evaluate by substituting for the 'x':
subst q0 1; // 0.0
subst q0 5; // 36.0
Time for generalization! We shall use the Horner's rule to transform
and evaluate polynomials of arbitrary degree, not just quadratic ones.
First, define the core operation:
horner v c = v*x+c;
Next, a function that, given a list of coefficients, creates a
polynomial in x (foldl1 is similar to J's '/' but works from left
to right):
poly cs = foldl1 horner cs;
Worth noting is that it is due to the horner function that the
polynomials created by poly are defined in x.
Now, e.g.,
p0 = poly [a,b,c];
p0; // (a*x+b)*x+c
is the 'hornered' version of a general quadratic polynomial, and
p1 = poly [2,-3,1];
p1; // (2*x+-3)*x+1
is the same as q1 above.
But we can indeed define polynomials of higher degrees, e.g.:
p2 = poly [1,-2,3,5];
p2; // ((1*x+-2)*x+3)*x+5
representing a 'hornered' version of x³-2x²+3x+5.
Let's do some evaluations, again using subst from above:
subst p0 5; // (a*5+b)*5+c
subst p1 5; // 36
subst p2 5; // 95
We may wish to tabulate a polynomial over lists of arguments. Here is
a helper function to simplify this:
tab p args = map (subst p) args;
Thus:
tab p1 (-5..5); // [66,45,28,15,6,1,0,3,10,21,36]
tab p2 (-5..5); // [-185,-103,-49,-17,-1,5,7,11,23,49,95]
Actually, we can do more than what poly does. With the following
definition:
polyhorn cs = init p, last p
when p = scanl1 horner cs end;
the Horner's rule is fully applied, because the partial results from
applying it along the list of coefficients cs are collected in a list
rather than thrown out as with poly (which is due to replacing foldl1
with scanl1). Thus, e.g.:
polyhorn [a,b,c]; // [a,a*x+b],(a*x+b)*x+c
or
polyhorn [a,b,c,d]; // [a,a*x+b,(a*x+b)*x+c],((a*x+b)*x+c)*x+d
In general, if cs is a list of coefficients representing a polynomial
p and v is any numeric value, then the expression
subst (polyhorn cs) v
is a pair of values: a list of coefficients representing a polynomial
q, and a number r=p(v), such that p(x) = (x-v).q(x)+r holds. In other
words, we can now divide a polynomial (p) by a linear binomial (x-v)
obtaining a quotient (q) and a remainder r, the latter being the value
of p at v. For example,
subst (polyhorn [1,-2,3,5]) 5; // [1,3,18],95
shows that x³-2x²+3x+5 = (x-5)*(x²+3x+18)+95, and that, in
particular, (x³-2x²+3x+5)(5) is 95.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
