function [x,d,V0,V,verifm]=resoud_system(A,b,epsi,choent)

% résolution d'un système linéaire avec structure de l'espace des solutions
%
% *********************************************************
%  [x,d,V0,V,verifm]=resoud_system(A,b,epsi,choent) :
%     Détermine l'espace affine des solutions avec aucun élément, un seul
%     élément ou une base affine d'un moins un élément
%
%  * Variables d'entrée
%     * A : matrice carrée numérique ou symbolique (n,n)
%     * b : vecteur colonne
%  * Variable d'entrée optionnelle :
%     * epsi : nombre positif de test d'égalité à zéro (pris égal à
%              eps*1e4 par défaut)
%     * choent : entier égal à 0 par défaut
%        * 0 : pas de manipulation d'entier
%        * 1 : manipulation d'entiers sur opérations sur les lignes
%  * Variables de sortie :
%     * x : vecteur colonne solution si elle est unique
%           vide sinon
%     * d : entier égal à 1 si la solution est unique, 0 sinon
%     * V0 : vecteur de R^n
%     * V :  matrice rectangulaire (n,p)
%            tels que si l'espace affine des solution est de dimension
%            supérieure à 1, il est de la forme
%            V0+sum_{i=1}^p x_i V(:,i) où  x_i est un réel.
%            V0+V est donc l'ensemble de p vecteurs indépendants de
%            l'espace des solutions
%  * Variable de sortie optionnelles :
% %   * verifm : réel positif ou nul dont la nullité traduit que le système
%       linéaire Ax=b est bien vérifié dans le cas où il y a plus d'une
%       solution. Vide sinon.
%
%  EXEMPLE :
%  a)
%     A=[1 2;3 4];B=[1 1].';[x,d,V0,V,verifm]=resoud_system(A,B)
%  b)
%     A=[1 2;2 4];B=[1 1].';[x,d,V0,V,verifm]=resoud_system(A,B)
%  c)
%     A=[1 2;2 4];B=[1 2].';[x,d,V0,V,verifm]=resoud_system(A,B)
%  d)
%     B=[1 2 3 5];
%     A=[9 7 1 8];
%     A=[B;B+A;A;2*B+3*A];
%     B=B.';
%     [x,d,V0,V,verifm]=resoud_system(A,B)
%  e)
%     b=[1 2 3 5];a=[9  7 1 8];
%     A=[b;b+a;a;2*b+3*a];
%     B=[sum(b);sum(a+b);sum(a);sum(2*b+3*a)];
%     [x,d,V0,V,verifm]=resoud_system(A,B)
%  f)
%     [x,d,V0,V,verifm]=resoud_system(ones(5,5),zeros(5,1))
%  g)
%     l1=[1 2 3];
%     l2=[2 4 1];
%     l3=[1 2 0];
%     A=[l1;l2;l3];
%     b=A*[3 2 1].';
%     [L,U,P]=lu(A)
%     [Ab,bb,permut]=pivot_gauss(A,b)
%
%
%
% ************ Fonctions auxiliaires utilisées ************
%
%    pivot_gauss
%
% *********************************************************
%
% 2014 by          Jérôme BASTIEN
%                  Université Claude Bernard Lyon I, UFRSTAPS, Laboratoire CRIS, Villeurbanne
%                  E-Mail : jerome.bastien@univ-lyon1.fr
%

nin=nargin;
nout=nargout;
if nin<=2|isempty(epsi)
    epsi=eps*1e4;
end
if nin<=3|isempty(choent)
    choent=0;
end
[n,no]=size(A);
if ~isnumeric(A)&isnumeric(b)
    b=sym(b);
end
if isnumeric(A)&~isnumeric(b)
    A=sym(A);
end
testsy=~isnumeric(A);
if choent==1
    xx=[A(:)-floor(A(:));b-floor(b)];
    if testsy
        [dodomax,dodoind]=max(abs(double(xx)));
        testlo=xx(dodoind)~=sym(0);
    else
        testlo=max(abs(xx))>epsi;
    end
    if testlo
        error('des coefficients non entiers');
    end
end
[Ab,bb,permut]=pivot_gauss(A,b,epsi,0,choent);
V=sum(abs(Ab),2);
if testsy
    N=max(find(diag(Ab)~=sym(0)));
else
    N=max(find(V>epsi));
end
if isempty(N)
    N=0;
end
if N==n
    x=A\b;
    er=A*x-b;
    d=1;
    V0=[];
    V=[];
    if nout>=5
        if testsy
            [mm,indmax]=max(abs(double(er)));
            verifm=er(indmax);
        else
            verifm=max(abs((er)));
        end
    end
else
    d=0;
    x=[];
    y2=bb(N+1:n);
    if testsy
        [mm,indmax]=max(abs(double(y2)));
        test=y2(indmax)~=sym(0);
    else
        test=max(abs(y2))>epsi;
    end
    if test
        V0=[];
        V=[];
        if nout>=5
            verifm=[];
        end
    else
        if N==0
            V0=zeros(n,1);
            V=eye(n);
            if nout>=5
                verifm=0;
            end
            if testsy
                V0=sym(V0);
                V=sym(V);
                verifm=sym(verifm);
            end
        else
            C=Ab(1:N,1:N);
            D=Ab(1:N,N+1:n);
            if testsy
                V0=[C\bb(1:N);sym(zeros(n-N,1))];
                V=[-C\D;sym(eye(n-N,n-N))];
            else
                V0=[C\bb(1:N);zeros(n-N,1)];
                V=[-C\D;eye(n-N,n-N)];
            end
            V0(permut,:)=V0;
            V(permut,:)=V;
            if nout>=5
                er=A*(V0(:,ones(1,n-N))+V)-b(:,ones(1,n-N));
                er=er(:);
                if testsy
                    [mm,indmax]=max(abs(double(er)));
                    verifm=er(indmax);
                else
                    verifm=max(abs(er));
                end
            end
        end
    end
end