share/bernstein/bernstein_utilities.mac

   1 /*
   2 ;;  Author Barton Willis
   3 ;;  University of Nebraska at Kearney
   4 ;;  Copyright (C) 2011 Barton Willis
   5
   6 ;;  This program is free software; you can redistribute it and/or modify
   7 ;;  it under the terms of the GNU General Public License as published by
   8 ;;  the Free Software Foundation; either version 2 of the License, or
   9 ;;  (at your option) any later version.
  10
  11 ;;  This program is distributed in the hope that it will be useful,
  12 ;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 ;;  GNU General Public License for more details.
  15 */
  16
  17 bernstein_expand(p,v) := block([d, dd, bbinomial,xcoeff, dm,q],
  18   bbinomial : lambda([a,b], xreduce("*", map(lambda([s,w], binomial(s,w)), a, b))),
  19   xcoeff : lambda([p,x,d], (p : ratexpand(p), map(lambda([z,n], p : ratcoeff(p,z,n)), x, d), p)),
  20   p : ratexpand(p),
  21   dm : map(lambda([s], hipow(p, s)),v),
  22   if not polynomialp(p,v, lambda([s], lfreeof(v,s))) then error("Input to bernstein_factor must be a polynomial"),
  23   d : map(lambda([s], map(lambda([w], hipow(s,w)), v)), if ?mplusp(p) then args(p) else [p]),
  24   dd : map(lambda([s], setify(makelist(k,k,0,s))),dm),
  25   dd : listify(apply('cartesian_product,dd)),
  26   q :  map(lambda([I], multibernstein_poly(I,dm,v) * lsum(bbinomial(I,J) * xcoeff(p, v, J)/ bbinomial(dm,J),J, d)),dd),
  27   xreduce("+",q))$