## Copyright (C) 1996, 1997 John W. Eaton
##
## This file is part of Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} {} polyder (@var{p})
## @deftypefnx {Function File} {} polyder (@var{p},@var{q})
## @deftypefnx {Function File} {[@var{num},@var{den}] = } polyder (@var{p},@var{q})
## Return the coefficients of the derivative of polynomials,
## polynomial quotients or polynomial products.
## If one argument is given, the derivative of the polynom @var{p} is 
## calculated.
## polyder(@var{p},@var{q}) calculates the derivative of the polynomial
## product. [@var{num},@var{den}] = polyder (@var{p},@var{q}) calculates
## the numerator @var{num} and the denominator @var{den} of the derivative
## of @var{p}/@var{q}.
## @seealso{polyderiv}
## @end deftypefn

## Author: jwe

## Modified: Kai Habel <kai.habel@gmx.de> (April 2000)
## Added support for polynomial products and quotients

function [...] = polyder (p,q)
	
	if nargin == 1
		if !is_vector(p)
			error('polyder(p): p must be a vector');
		endif
		#polynom
		lp=length(p);
		if lp==1
			vr_val(0);
			return;
		endif
		tmp = p (1:(lp-1)) .* [(lp-1):-1:1];
		vr_val(tmp);
	elseif nargin == 2
		if !(is_vector(p)&is_vector(q))
			error('polyder(p,q): p and q must be vectors');
		endif

		lp=length(p);
		lq=length(q);

		p1=0;
		if lp>1
			p1 = p (1:(lp-1)) .* [(lp-1):-1:1];
		endif
	
		q1=0;
		if lq>1
			q1 = q (1:(lq-1)) .* [(lq-1):-1:1];
		endif

		pq1=conv(p,q1);
		qp1=conv(q,p1);
		if (length(pq1)>length(qp1))
			qp1=prepad(qp1,length(pq1));
		elseif (length(pq1)<length(qp1))
			pq1=prepad(pq1,length(qp1));
		endif

		if nargout <= 1
			#product
			tmp=pq1+qp1;

			#remove heading zeros
			if any(tmp)
				tmp([ 1:find(tmp!=0)-1 ])=[];
			else
				tmp=0;
			endif
			vr_val(tmp);	
		elseif nargout == 2
			#quotient
			num=qp1-pq1;
			den=conv(q,q);

			#remove heading zeros
			if any(num)
				num([ 1:find(num!=0)-1 ])=[];
			else
				num=0;
			endif

			if any(den)
				den([ 1:find(den!=0)-1 ])=[];
			else
				den=0;
			endif

			vr_val(num);
			vr_val(den);
		endif
	else
		error('polyder: number of arguments must be 1 or 2');
	endif
endfunction
