% Copyright 2015 by Jonathan D. Hauenstein

% Deflate Caprasse system from Section 5.2 of 
% "Certifying isolated singular points and their multiplicity structure" 
% by J.D. Hauenstein, B. Mourrain, and A. Szanto

function [Fext,Vext,Fint,Vint,Fiso,Viso,Fiso2,Viso2,Fi,Vi] = deflate52_sys()
    % Output: Fext - polynomial system from extrinsic slicing deflation
    %         Vext - variables arising from extrinsic slicing deflation
    %         Fnf - polynomial system from intrinsic slicing deflation
    %         Vnf - variables arising from intrinsic slicing deflation
    %         Fiso - polynomial system from isosingular deflation
    %         Viso - variables arising from isosingular deflation
    %         Fiso2 - polynomial system from isosingular deflation with full
    %                   rank submatrix
    %         Viso2 - variables arising from isosingular deflation with full
    %                   rank submatrix
    %         Fi - polynomial system from |i|=1 deflation
    %         Vi - variables arising from |i|=1 deflation
    
    % setup the original system
    syms x1 x2 x3 x4;
    X = [x1 x2 x3 x4];
    F = [x1^3*x3-4*x1^2*x2*x4-4*x1*x2^2*x3-2*x2^3*x4-4*x1^2-4*x1*x3+10*x2^2+10*x2*x4-2; 
    	 x1*x3^3-4*x1*x3*x4^2-4*x2*x3^2*x4-2*x2*x4^3-4*x1*x3+10*x2*x4-4*x3^2+10*x4^2-2; 
         2*x1*x2*x4+x2^2*x3-2*x1-x3; 
         x1*x4^2+2*x2*x3*x4-x1-2*x3];
    
    % jacobian of F at [2,-sqrt(-3),2,sqrt(-3)] has a 2-dim'l null space
    J = jacobian(F,X);
    
    % extrinsic slicing deflation
    syms l1 l2 l3 l4;
    L = [l1;l2;l3;l4];
    Fext = [F;J*L;randi([-10,10],2,4)*L-ones(2,1)];
    Vext = [X transpose(L)];

    % display summary of extrinsic slicing deflation
    disp(['Extrinsic slicing deflation']);
    disp(['  functions: ',num2str(length(Fext))]);
    disp(['  variables: ',num2str(length(Vext))]);
        
    % intrinsic slicing deflation
    L = [l1;l2];
    Fint = [F;J*(randi([-10,10],4,2)*L + randi([-10,10],4,1))];
    Vint = [X transpose(L)];
    
    % display summary of intrinsic slicing deflation
    disp(['Intrinsic slicing deflation']);
    disp(['  functions: ',num2str(length(Fint))]);
    disp(['  variables: ',num2str(length(Vint))]);    
    
    % isosingular deflation -- add all 3x3 minors
    Fiso = F;
    for j = 1:4 
        for k = 1:4
            Fiso = [Fiso;det(J([1:j-1,j+1:4],[1:k-1,k+1:4]))];
        end;
    end;
    Viso = X;

	% display summary of isosisngular deflation
    disp(['Isosingular deflation']);
    disp(['  functions: ',num2str(length(Fiso))]);
    disp(['  variables: ',num2str(length(Viso))]);    

	% isosingular deflation -- add all 3x3 minors containg [1:2,1:2]
    Fiso2 = F;
    for j = 3:4 
        for k = 3:4
            Fiso2 = [Fiso2;det(J([1,2,j],[1,2,k]))];
        end;
    end;
    Viso2 = X;

	% display summary of isosingular deflation
    disp(['Isosingular deflation with full rank submatrix']);
    disp(['  functions: ',num2str(length(Fiso2))]);
    disp(['  variables: ',num2str(length(Viso2))]);    

    % |i|=1 deflation
    A = J(1:2,1:2);
    B = J(1:2,3:4);
    N = det(A)*[-inv(A)*B;eye(2,2)]*randi([1,10],2,1); % random null vec    
    Fi = [F;J(3:4,:)*N];
    Vi = X;

	% display summary of |i|=1 deflation
    disp(['|i|=1 deflation']);
    disp(['  functions: ',num2str(length(Fi))]);
    disp(['  variables: ',num2str(length(Vi))]);        
    
return;

