J's #. implements horner's rule. FYI,
-- Raul On Wed, Dec 12, 2012 at 4:17 PM, Boyko Bantchev <[email protected]> wrote: > 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
