(* :Name: Algebra`Horner` *)

(* :Title: Horner Form Factorization Of Polynomials. *)

(* :Author: Mark Sofroniou *)

(* :Summary:
This package adds the function Horner for factoring univariate and
multivariate polynomials in efficient computational form.
*)

(* :Context: Algebra`Horner` *)

(* :Package Version: 1.1 *)

(* :Copyright: Copyright 1995-1999, Wolfram Research, Inc.
*)

(* :History:
 Tidied up the interface, February 1997.
 Rewritten for release with version 3.0, November 1995.
 (Original Version by Mark Sofroniou, September 1993.) *)

(* :Keywords:
 Horner, Nested, Polynomial. *)

(* :Source:
 D. E. Knuth, The Art of Computer Programming  Vol 2: Seminumerical Algorithms,
 Second Edition, Addison Wesley, New York, 1981. *)

(* :Mathematica Version: 3.0 *)

(* :Limitations:
 Stack recursion limits. *)

(* :Discussion:
Pre-conditioning of the coefficients can sometimes improve numerical stability
and further reduce the operation count.  Various other considerations arise such
as if the polynomial is not dense, if we require evaluations on an evenly spaced
grid, or if computations are performed in parallel.

The continued fraction form is usually preferable for rational functions
of polynomials.
*)

BeginPackage["Algebra`Horner`"]

Horner::usage = "Horner[p] gives the polynomial or rational function of
polynomials p in Horner or nested form with respect to the default variables
Variables[p].  Horner[p, v] gives p in Horner form with respect
to the variable or variable list v.  Horner[p1/p2, v1, v2] gives the
rational function of polynomials p1/p2 in Horner form using the variables
or variable lists v1 and v2 for p1 and p2 respectively."

Unprotect[Horner];

Begin["`Private`"]

(* Currently supported forms are polynomials and rational functions of
polynomials *)

expand[p_/q_]:= Expand[ p ]/Expand[ q ];
expand[p_]:= Expand[ p ];

(* Horner's rule. Top level function does error checking. *)

(* Default variables. *)

Horner[ p_/q_ ]:=
  Block[{$RecursionLimit=Infinity},
    Module[{ans, ep, eq},
      ep = Expand[p];
      eq = Expand[q];
      ans  /; (ans = Catch[ HornerRule[ ep, Variables[ep] ]/
                 HornerRule[ eq, Variables[eq] ] ];
                   ans=!=$Failed)
    ]
  ];

Horner[ p_ ]:=
  Block[{$RecursionLimit=Infinity},
    Module[{ans, ep},
      ep = expand[p];
      ans /; (ans = Catch[ HornerRule[ ep, Variables[ep] ] ];
                  ans=!=$Failed)
    ]
  ] /; !MatchQ[ p, HoldPattern[_/_]];


(* Specified variables. *)

VarQ[x_]:= VectorQ[x] || Head[x]=!=List;
ToList[v_List]:= v;
ToList[v_]:= {v};

Horner[ p_, v_?VarQ ]:=
  Block[{$RecursionLimit=Infinity},
    Module[{ans},
      ans /; (ans = Catch[ HornerRule[ expand[p], ToList[v] ] ];
                  ans=!=$Failed)
    ]
  ];

Horner[ p_/q_, vp_?VarQ, vq_?VarQ ]:=
  Block[{$RecursionLimit=Infinity},
    Module[{ans},
      ans /; (ans = Catch[ HornerRule[ Expand[p], ToList[vp] ]/
                  HornerRule[ Expand[q], ToList[vq] ] ];
                ans=!=$Failed)
    ]
  ];

(* Rational function of polynomials *)

HornerRule[ p_/q_, v_ ]:= HornerRule[p, v]/HornerRule[q, v];

(* End of recursion, or no variables found by Variables. *)

HornerRule[ p_, {} ]:= p;

(* Nothing to Horner *)

HornerRule[ 0, _ ] = 0;

(* This uses sparse polynomial coefficient code, which returns a list
 of pairs of the form {exponent, coefficient} in decreasing order of
 exponent. *)

sparseCoefficientList = Internal`NestedTermsList;

(* Horner's rule for univariate polynomials. Variable is paired with
 coefficients as {var^exp, coeff}. *)

HornerRule[ p_, {v_} ]:=
  Fold[
    ConstructSum,
    0,
    PowersAndCoeffs[ sparseCoefficientList[ p, v ], v ]
  ];

(* Horner's rule for multivariate polynomials as recursive univariate
 decomposition using indirect recursion. Variable is paired with coefficients
 as {var^exp,hcoeff} where hcoeff is the coefficient in Horner form. *)

HornerRule[ p_, {v_, rest__} ]:=
  Fold[
    ConstructSum,
    0,
    PowersAndCoeffs[
      sparseCoefficientList[ p, v ], v, HornerRule[ #, {rest} ]&
    ]
  ];

(* Check for invalid polynomial. *)

Horner::fail = "The expression `1` is not a valid polynomial in the variable `2`";

PowersAndCoeffs[ _?(SameQ[#, sparseCoefficientList]&)[ p_, v_ ], __ ]:=
  (Message[ Horner::fail, p, v ]; Throw[ $Failed ]);

(* Powers are ordered in descending order. Calculate offset powers as
 successive differences (not necessarily incremental). *)

PowersAndCoeffs[ eclist_, v_ ]:=
  Thread[ {OffsetPowers[ v, Part[#, 1] ], Part[#, 2]} ]& @ Thread[eclist];

(* This is the recursive part for multivariates. *)

PowersAndCoeffs[ eclist_, v_, RecurseHorner_ ]:=
  Thread[ {OffsetPowers[ v, Part[#, 1] ], Map[RecurseHorner, Part[#, 2] ]} ]& @
    Thread[eclist];

(* Calculate powers of variable as difference between successive exponents.
 Use Listable attribute of Power. *)

OffsetPowers[ v_, e_ ]:= v^(e - Join[ Drop[e, 1], {0} ]);

(* Accumulate Horner form of polynomial from running sum, variable^exponent
 and coefficient. *)

ConstructSum[ sum_, {ve_, c_} ]:= ve*(c + sum);


End[]; (* End `Private` context. *)

(* Protect exported symbols. *)

SetAttributes[Horner, {Protected, ReadProtected}];

EndPackage[]; (* End package context. *)
