equ2vect(a, b, c, r)

function [x, a, b, c]= equ2vect(a, b, c, r)
% equ2vect  quation relle du second degr vectorise
%  a b c sont des vecteurs de mme taille ou des scalaires
%
%  a(i)x^2 +b(i)x+c(i) est rsolue dans R pour tous les i
%  si l'quation est du premier degr, la seconde solution est NaN
%  si la chane r commence par un r ou est absente, seules les solutions
%  relles sont calcules  moins que l'un des coefficients des quations
%  ne soit complexe.
%  si seules les racines relles sont calcules et que l'quation n'en
%  possde pas les deux solutions sont NaN
%
%  En sortie a, b, c sont optionnels et sont des vecteurs colonnes de mme
%  taille.

if nargin <= 3
  r = 'r';
end; 
n = max([length(a); length(b); length(c)]);
delta = ones(n,1).*NaN;
x = [delta, delta];

a = a(:);
if length(a) == 1
  a = a(ones(size(delta)),:);
end;
b = b(:);
if length(b) == 1
  b = b( ones(size(delta)),:);
end;
c = c(:);
if length(c) == 1
  c = c(ones(size(delta)),:);
end;

if length(a) ~= length(b) | length(a) ~= length(c)
  error('a, b et c doivent tre scalaires ou de mme dimension');
end;

r = (upper(r(1)) == 'R') & ...
    all(isreal(a)) & all(isreal(b)) & all(isreal(c));

test = (a == 0);                                            % premier degr
ind = find(test);
x(ind,1) = -c(ind)./b(ind);
x(ind,2) = Inf*x(ind,1);

test = ~test;                                               % second degr
delta = NaN.*ones(n, 1);
delta(test) = b(test).*b(test)-4.*a(test).*c(test);

if r
  test = (delta >= 0);                                      % delta positif
end;  
delta(test) = sqrt(delta(test));
ind = find(test);
x(ind,1) = [0.5.*(-b(ind)-(sign(real(b(ind)))+(b(ind)==0)).*delta(ind))./a(ind)];

ind = find(test & x(:,1) ~=0);
x(ind,2) = [(c(ind,1)./a(ind,1))./x(ind,1)];
ind = find(test & x(:,1) ==0);
x(ind,2) = zeros(size(ind));