octave_packages/symbolic-1.1.0/sym2poly.m

   1 ## Copyright (C) 2003 Willem J. Atsma <watsma@users.sf.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or
   4 ## modify it under the terms of the GNU General Public
   5 ## License as published by the Free Software Foundation;
   6 ## either version 2, or (at your option) any later version.
   7 ##
   8 ## This software is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied
  10 ## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11 ## PURPOSE.  See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public
  15 ## License along with this software; see the file COPYING.  If not,
  16 ## see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn {Function File} {@var{c} =} sym2poly (@var{p}, @var{x})
  20 ## Returns the coefficients of the symbolic polynomial expression @var{p}
  21 ## as a vector. If there is only one free variable in @var{p} the
  22 ## coefficient vector @var{c} is a plain numeric vector. If there is more
  23 ## than one free variable in @var{p}, a second argument @var{x} specifies the
  24 ## free variable and the function returns a cell vector of symbolic expressions.
  25 ## The coefficients correspond to decreasing exponent of the free variable.
  26 ##
  27 ## Example:
  28 ## @example
  29 ## symbols
  30 ## x = sym("x");
  31 ## y = sym("y");
  32 ## c = sym2poly (x^2+3*x-4);    # c = [1,3,-4]
  33 ## c = sym2poly (x^2+y*x,x);    # c = @{2,y,0@}
  34 ## @end example
  35 ##
  36 ## If @var{p} is not a polynomial the result has no warranty.
  37 ##
  38 ## @seealso{poly2sym,polyval,roots}
  39 ## @end deftypefn
  40
  41 ## Created: 18 April 2003
  42 ## Changed: 25 April 2003
  43 ##    Removed the use of differentiate to get to coefficients - round-off
  44 ##     errors cause problems. Now using newly created sumterms().
  45 ## Changed: 6 May 2003
  46 ##    Removed the attempt to use ldegree(), degree() and coeff() - results
  47 ##     with these are inconsistent.
  48
  49 function c = sym2poly(p,x)
  50
  51   BADPOLY_COEFF_LIMIT = 500;
  52
  53   if is_vpa(p)
  54     ## polynomial is one vpa number
  55     c = to_double(p);
  56     if length(c)!=1
  57       error("Argument is not a polynomial.");
  58     endif
  59     return
  60   endif
  61
  62   if !is_ex(p)
  63     error("Argument has to be a symbolic expression.")
  64   endif
  65
  66   pvars = findsymbols(p);
  67   if isempty(pvars)
  68     ## It is possible that we get an expression without any symbols.
  69     c = to_double(p);
  70     return;
  71   endif
  72   nvars = length(pvars);
  73
  74   if nvars>1 && exist("x")!=1
  75     error("Symbolic expression has more than 1 free variable; no variable specified.")
  76   elseif exist("x")!=1
  77     x = pvars{1};
  78   endif
  79
  80   p = expand(p);
  81
  82   ## GiNaC has commands to access coefficients directly, but in octave this often
  83   ## does not work, because for example x^2 typed in octave results in a
  84   ## non-integer power in GiNaC: x^2.0 .
  85
  86   [num,den] = numden(p);
  87   tmp = findsymbols(den);
  88   for i=1:length(tmp)
  89     if tmp{i}==x
  90       error("Symbolic expression is a ratio of polynomials.")
  91     endif
  92   endfor
  93
  94   p = expand(p);
  95   p_terms = sumterms(p);
  96   ## if this is well behaved, I can find the coefficients by dividing with x
  97   c_ex = cell;
  98   for i=1:length(p_terms)
  99     tmp = p_terms{i};
 100     for j=1:BADPOLY_COEFF_LIMIT
 101       if disp(differentiate(tmp,x))=="0"
 102         break;
 103       endif
 104       tmp = tmp/x;
 105     endfor
 106     if j==BADPOLY_COEFF_LIMIT
 107       printf("Please examine your code or adjust this function.\n");
 108       printf("This error may occur because the passed expression is not a polynomial.\n");
 109       error("Reached the set limit (%d) for the number of coefficients.",BADPOLY_COEFF_LIMIT)
 110     endif
 111     if (length(c_ex)<j) || isempty(c_ex{j})
 112       c_ex(j)=tmp;
 113     else
 114       c_ex(j) = c_ex{j}+tmp;
 115     endif
 116   endfor
 117   order = length(c_ex)-1;
 118
 119   all_numeric = true;
 120   for i=1:(order+1)
 121     if isempty(c_ex{i})
 122       c_ex(i) = vpa(0);
 123     endif
 124     cvar=findsymbols(c_ex{i});
 125     ncvar = length(cvar);
 126     if ncvar
 127       all_numeric=false;
 128       for j=1:ncvar
 129         if disp(cvar{j})==disp(x)
 130           printf("Possibly this error occurs because two symbols with the same name\n");
 131           printf("are different to GiNaC. Make sure the free variable exists as a\n");
 132           printf("symbol in your workspace.\n");
 133           error("The symbolic expression is not a polynomial.")
 134         endif
 135       endfor
 136     endif
 137   endfor
 138
 139   c_ex = reverse(c_ex);
 140
 141   if all_numeric
 142     for i=1:(order+1)
 143       c(1,i)=to_double(c_ex{i});
 144     endfor
 145   else
 146     c = c_ex;
 147   endif
 148
 149 endfunction
 150
 151 %!shared
 152 % symbols
 153 % x=sym("x"); y=sym("y");
 154 %!test
 155 % assert(sym2poly (x^2+3*x-4), [1,3,-4]);
 156 %!test
 157 % assert (sym2poly (x^2+y*x,x), {2,y,0})