(**** LECTURE 12: Polynomial Arithmetic ****) (* $Id: polynom.ML,v 1.2 1998/09/08 14:42:44 lcp Exp lcp $ *) (*Polynomials in one variable (univariate polynomials) are lists of (exponent, coefficient) pairs. Exponents are decreasing and coeeficients are non-zero For exact results, coefficients should be rational rather than real. $Id: polynom.ML,v 1.2 1998/09/08 14:42:44 lcp Exp lcp $ **) (*Sum of two univariate polynomials*) fun polysum [] us = us : (int*real)list | polysum ts [] = ts | polysum ((m,a)::ts) ((n,b)::us) = if m>n then (m,a) :: polysum ts ((n,b)::us) else if n>m then (n,b) :: polysum us ((m,a)::ts) else (*m=n*) if a+b=0.0 then polysum ts us else (m, a+b) :: polysum ts us; (*Product of two terms*) fun termprod (m,a) (n,b) = (m+n, a*b) : (int*real); (*Product of two univariate polynomials*) fun polyprod [] us = [] | polyprod ((m,a)::ts) us = polysum (map (termprod(m,a)) us) (polyprod ts us); (*General exponentiation of the function f: defined only for k>0*) fun funpower f x k = if k=1 then x else if k mod 2 = 0 then funpower f (f x x) (k div 2) else f x (funpower f (f x x) (k div 2)); (*Raising a polynomial to a positive integer power*) val polypower = funpower polyprod; (*Convert an unsorted term list to a polynomial; like merge sort*) fun makepoly [] = [] | makepoly [(m,a)] = if a=0.0 then [] else [(m,a)] | makepoly ts = let val k = length ts div 2 in polysum (makepoly (take(ts,k))) (makepoly (drop(ts,k))) end; (*Make a polynomial from a list of integer coefficients for all terms. e.g. transforms [1,2,3] to x^2 + 3 *) fun coeffpoly [] = [] | coeffpoly(0::ns) = coeffpoly ns | coeffpoly(n::ns) = (length ns, real n) :: coeffpoly ns; val xplus1 = coeffpoly [1,1] and xminus1 = coeffpoly [1,~1]; (** Display of a polynomial as a string **) fun coeffshow a = if a=1.0 then "" else makestring a; fun termshow (0,a) = makestring a | termshow (1,a) = coeffshow a ^ "x" | termshow (m,a) = coeffshow a ^ "x^" ^ makestring m; fun polyshowing [] = "" | polyshowing((m,a)::ts) = if a<0.0 then " - " ^ termshow(m,~a) ^ polyshowing ts else " + " ^ termshow(m, a) ^ polyshowing ts; fun polyshow [] = "" | polyshow((m,a)::ts) = if a<0.0 then " - " ^ termshow(m,~a) ^ polyshowing ts else termshow(m, a) ^ polyshowing ts; (** Timings of the multiplication function for large powers **) polypower xplus1 5; polypower xplus1 10; polypower xplus1 32; polypower xplus1 100; (*Non GC 0.0742 GC 0.0 both 0.0742*) polypower xplus1 200; (*Non GC 0.2890 GC 0.0429 both 0.3320*) polypower xplus1 300; (*Non GC 0.6562 GC 0.5194 both 1.1757*) (** All timings were on orwell. They are worse than quadratic. val p400= polypower xplus1 400; (*Non GC 1.2186 GC 1.0155 both 2.2342*) polypower xplus1 800; (*Non GC 5.3004 GC 12.229 both 17.530*) polyprod p400 p400; (*Non GC 3.5427 GC 9.8079 both 13.3507*) **) (*Product of two univariate polynomials, by balanced merging*) fun polyprod2 [] us = [] | polyprod2 [(m,a)] us = map (termprod(m,a)) us | polyprod2 ts us = let val k = length ts div 2 in polysum (polyprod2 (take(ts,k)) us) (polyprod2 (drop(ts,k)) us) end; val polypower2 = funpower polyprod2; polypower2 xplus1 100; (*Non GC 0.0742 GC 0.0 both 0.0742*) polypower2 xplus1 200; (*Non GC 0.2851 GC 0.0 both 0.2851*) polypower2 xplus1 300; (*Non GC 0.6366 GC 0.0234 both 0.6601*) (** All timings were on orwell. They are slightly worse than quadratic. polypower2 xplus1 400; (*Non GC 1.1796 GC 0.0742 both 1.2538*) polypower2 xplus1 800; (*Non GC 4.9371 GC 0.4609 both 5.3980*) polyprod2 p400 p400; (*Non GC 3.4138 GC 0.2929 both 3.7067*) **) (** Division: quotient and remainder **) (*Division by zero -- empty polynomial -- raises exception Match*) fun polyquorem ts ((n,b)::us) = let fun quo [] qs = (rev qs, []) | quo ((m,a)::ts) qs = if m