--- /dev/null
+## Copyright (C) 2003 Willem J. Atsma <watsma@users.sf.net>
+##
+## This program is free software; you can redistribute it and/or
+## modify it under the terms of the GNU General Public
+## License as published by the Free Software Foundation;
+## either version 2, or (at your option) any later version.
+##
+## This software is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied
+## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+## PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public
+## License along with this software; see the file COPYING. If not,
+## see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{c} =} sym2poly (@var{p}, @var{x})
+## Returns the coefficients of the symbolic polynomial expression @var{p}
+## as a vector. If there is only one free variable in @var{p} the
+## coefficient vector @var{c} is a plain numeric vector. If there is more
+## than one free variable in @var{p}, a second argument @var{x} specifies the
+## free variable and the function returns a cell vector of symbolic expressions.
+## The coefficients correspond to decreasing exponent of the free variable.
+##
+## Example:
+## @example
+## symbols
+## x = sym("x");
+## y = sym("y");
+## c = sym2poly (x^2+3*x-4); # c = [1,3,-4]
+## c = sym2poly (x^2+y*x,x); # c = @{2,y,0@}
+## @end example
+##
+## If @var{p} is not a polynomial the result has no warranty.
+##
+## @seealso{poly2sym,polyval,roots}
+## @end deftypefn
+
+## Created: 18 April 2003
+## Changed: 25 April 2003
+## Removed the use of differentiate to get to coefficients - round-off
+## errors cause problems. Now using newly created sumterms().
+## Changed: 6 May 2003
+## Removed the attempt to use ldegree(), degree() and coeff() - results
+## with these are inconsistent.
+
+function c = sym2poly(p,x)
+
+ BADPOLY_COEFF_LIMIT = 500;
+
+ if is_vpa(p)
+ ## polynomial is one vpa number
+ c = to_double(p);
+ if length(c)!=1
+ error("Argument is not a polynomial.");
+ endif
+ return
+ endif
+
+ if !is_ex(p)
+ error("Argument has to be a symbolic expression.")
+ endif
+
+ pvars = findsymbols(p);
+ if isempty(pvars)
+ ## It is possible that we get an expression without any symbols.
+ c = to_double(p);
+ return;
+ endif
+ nvars = length(pvars);
+
+ if nvars>1 && exist("x")!=1
+ error("Symbolic expression has more than 1 free variable; no variable specified.")
+ elseif exist("x")!=1
+ x = pvars{1};
+ endif
+
+ p = expand(p);
+
+ ## GiNaC has commands to access coefficients directly, but in octave this often
+ ## does not work, because for example x^2 typed in octave results in a
+ ## non-integer power in GiNaC: x^2.0 .
+
+ [num,den] = numden(p);
+ tmp = findsymbols(den);
+ for i=1:length(tmp)
+ if tmp{i}==x
+ error("Symbolic expression is a ratio of polynomials.")
+ endif
+ endfor
+
+ p = expand(p);
+ p_terms = sumterms(p);
+ ## if this is well behaved, I can find the coefficients by dividing with x
+ c_ex = cell;
+ for i=1:length(p_terms)
+ tmp = p_terms{i};
+ for j=1:BADPOLY_COEFF_LIMIT
+ if disp(differentiate(tmp,x))=="0"
+ break;
+ endif
+ tmp = tmp/x;
+ endfor
+ if j==BADPOLY_COEFF_LIMIT
+ printf("Please examine your code or adjust this function.\n");
+ printf("This error may occur because the passed expression is not a polynomial.\n");
+ error("Reached the set limit (%d) for the number of coefficients.",BADPOLY_COEFF_LIMIT)
+ endif
+ if (length(c_ex)<j) || isempty(c_ex{j})
+ c_ex(j)=tmp;
+ else
+ c_ex(j) = c_ex{j}+tmp;
+ endif
+ endfor
+ order = length(c_ex)-1;
+
+ all_numeric = true;
+ for i=1:(order+1)
+ if isempty(c_ex{i})
+ c_ex(i) = vpa(0);
+ endif
+ cvar=findsymbols(c_ex{i});
+ ncvar = length(cvar);
+ if ncvar
+ all_numeric=false;
+ for j=1:ncvar
+ if disp(cvar{j})==disp(x)
+ printf("Possibly this error occurs because two symbols with the same name\n");
+ printf("are different to GiNaC. Make sure the free variable exists as a\n");
+ printf("symbol in your workspace.\n");
+ error("The symbolic expression is not a polynomial.")
+ endif
+ endfor
+ endif
+ endfor
+
+ c_ex = reverse(c_ex);
+
+ if all_numeric
+ for i=1:(order+1)
+ c(1,i)=to_double(c_ex{i});
+ endfor
+ else
+ c = c_ex;
+ endif
+
+endfunction
+
+%!shared
+% symbols
+% x=sym("x"); y=sym("y");
+%!test
+% assert(sym2poly (x^2+3*x-4), [1,3,-4]);
+%!test
+% assert (sym2poly (x^2+y*x,x), {2,y,0})
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff