% Copyright 2015 by Jonathan D. Hauenstein

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

function [FA,VA,FB,VB,FC,VC,FD,VD] = deflate53_3()
    % Output: FA - polynomial system from method A
    %         VA - variables arising from method A
    %         FB - polynomial system from method B
    %         VB - variables arising from method B
    %         FC - polynomial system from method C
    %         VC - variables arising from method C
    %         FD - polynomial system from method D
    %         VD - variables arising from method D

    [FA,VA] = deflate53_3_A();
    
    % display summary of method A
    disp(['Method A']);
    disp(['  functions: ',num2str(length(FA))]);
    disp(['  variables: ',num2str(length(VA))]);    
    
    [FB,VB] = deflate53_3_B();
    
    % display summary of method B
    disp(['Method B']);
    disp(['  functions: ',num2str(length(FB))]);
    disp(['  variables: ',num2str(length(VB))]);    
    
    [FC,VC] = deflate53_3_C();
    
    % display summary of method C
    disp(['Method C']);
    disp(['  functions: ',num2str(length(FC))]);
    disp(['  variables: ',num2str(length(VC))]);    

	[FD,VD] = deflate53_3_D();
    
    % display summary of method D
    disp(['Method D']);
    disp(['  functions: ',num2str(length(FD))]);
    disp(['  variables: ',num2str(length(VD))]);    

return;

function [FA,VA] = deflate53_3_A()
    % deflate using method A -- intrinsic slicing deflation

    syms x y;
    X = [x y];

    F = [14*x + 33*y - 3*sqrt(5)*(x^2 + 4*x*y + 4*y^2 + 2) + sqrt(7) + x^3 + 6*x^2*y + 12*x*y^2 + 8*y^3;
         41*x - 18*y - sqrt(5) + 8*x^3 - 12*x^2*y + 6*x*y^2 - y^3 + 3*sqrt(7)*(4*x*y - 4*x^2 - y^2 - 2)];

    % first iteration -- Jacobian has 1 dim'l null space at Pt = [5^(1/2)/5 + (2*7^(1/2))/5, (2*5^(1/2))/5 - 7^(1/2)/5];
    J = jacobian(F,X);

    % deflate
    syms u;
    X = [X u];
    Fd = [F;J*(randi([-10,10],2,1)*u + randi([-10,10],2,1))];

    % second iteration -- Jacobian still has 1 dim'l null space
    Jd = jacobian(Fd,X);
    
    % deflate
    syms v1 v2;
    X = [X v1 v2];
    Fdd = [Fd;Jd*(randi([-10,10],3,2)*[v1;v2] + randi([-10,10],3,1))];

    % third iteration -- Jacobian still has 1 dim'l null space
    Jdd = jacobian(Fdd,X);
    
    % deflate
    syms w1 w2 w3 w4;
    X = [X w1 w2 w3 w4];
    Fddd = [Fdd;Jdd*(randi([-10,10],5,4)*[w1;w2;w3;w4] + randi([-10,10],5,1))];

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);

    % deflate
    syms r1 r2 r3 r4 r5 r6 r7 r8;
    X = [X r1 r2 r3 r4 r5 r6 r7 r8];
    Fdddd = [Fddd;Jddd*(randi([-10,10],9,8)*[r1;r2;r3;r4;r5;r6;r7;r8] + randi([-10,10],9,1))];    
    
    % jacobian is full rank -- find distinct nonzero polynomials
    FA = unique(Fdddd);
    FA = FA(subs(FA,X,rand(size(X))) ~= 0);
    VA = X;
    
return;

function [FB,VB] = deflate53_3_B()
    % deflate using method B -- isosingular defl. with full rank submatrix

    syms x y;
    X = [x y];

    F = [14*x + 33*y - 3*sqrt(5)*(x^2 + 4*x*y + 4*y^2 + 2) + sqrt(7) + x^3 + 6*x^2*y + 12*x*y^2 + 8*y^3;
         41*x - 18*y - sqrt(5) + 8*x^3 - 12*x^2*y + 6*x*y^2 - y^3 + 3*sqrt(7)*(4*x*y - 4*x^2 - y^2 - 2)];

    % first iteration -- Jacobian has 1 dim'l null space at Pt = [5^(1/2)/5 + (2*7^(1/2))/5, (2*5^(1/2))/5 - 7^(1/2)/5];
    J = jacobian(F,X);

    % deflate -- add det(J)
    Fd = [F;det(J)];
    
    % second iteration -- Jacobian still has 1 dim'l null space
    Jd = jacobian(Fd,X);
    
    % deflate -- add 2x2 minors containing (1,1)
    Fdd = Fd;
    for j = 2:3
        Fdd = [Fdd;det(Jd([1,j],:))];
    end;

    % third iteration -- Jacobian still has 1 dim'l null space
    Jdd = jacobian(Fdd,X);

    % deflate -- add 2x2 minors containing (1,1)
    Fddd = Fdd;
    for j = 2:5
        Fddd = [Fddd;det(Jdd([1,j],:))];
    end;

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);
    
    % deflate -- add 2x2 minors containing (1,1)
    Fdddd = Fddd;
    for j = 2:9
        Fdddd = [Fdddd;det(Jddd([1,j],:))];
    end;

    % jacobian is full rank -- find distinct nonzero polynomials
    FB = unique(Fdddd);
    FB = FB(subs(FB,X,rand(size(X))) ~= 0);
    VB = X;
    
return;

function [FC,VC] = deflate53_3_C()
    % deflate using method C -- "kerneling" deflation

    syms x y;
    X = [x y];

    F = [14*x + 33*y - 3*sqrt(5)*(x^2 + 4*x*y + 4*y^2 + 2) + sqrt(7) + x^3 + 6*x^2*y + 12*x*y^2 + 8*y^3;
         41*x - 18*y - sqrt(5) + 8*x^3 - 12*x^2*y + 6*x*y^2 - y^3 + 3*sqrt(7)*(4*x*y - 4*x^2 - y^2 - 2)];

    % first iteration -- Jacobian has 1 dim'l null space at Pt = [5^(1/2)/5 + (2*7^(1/2))/5, (2*5^(1/2))/5 - 7^(1/2)/5];
    J = jacobian(F,X);
    A = J(1,1);
    B = J(1,2);
    C = J(2,1);
    D = J(2,2);
    N = det(A)*(D - C*inv(A)*B);

    % deflate
    Fd = [F;reshape(N,[numel(N),1])];

    % second iteration -- Jacobian still has 1 dim'l null space
    Jd = jacobian(Fd,X);
    A = Jd(1,1);
    B = Jd(1,2);
    C = Jd([2:end],1);
    D = Jd([2:end],2);
    N = det(A)*(D - C*inv(A)*B);
    
    % deflate
    Fdd = [Fd;reshape(N,[numel(N),1])];

    % third iteration -- Jacobian still has 1 dim'l null space
    Jdd = jacobian(Fdd,X);    
	A = Jdd(1,1);
    B = Jdd(1,2);
    C = Jdd([2:end],1);
    D = Jdd([2:end],2);
    N = det(A)*(D - C*inv(A)*B);
    
    % deflate
    Fddd = [Fdd;reshape(N,[numel(N),1])];

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);    
	A = Jddd(1,1);
    B = Jddd(1,2);
    C = Jddd([2:end],1);
    D = Jddd([2:end],2);
    N = det(A)*(D - C*inv(A)*B);    
    
    % deflate
    Fdddd = [Fddd;reshape(N,[numel(N),1])];
    
    % jacobian is full rank -- find distinct nonzero polynomials
    FC = unique(Fdddd);
    FC = FC(subs(FC,X,rand(size(X))) ~= 0);
    VC = X;
    
return;

function [FD,VD] = deflate53_3_D()
    % deflate using method D -- |i|=1 deflation

    syms x y;
    X = [x y];

    F = [14*x + 33*y - 3*sqrt(5)*(x^2 + 4*x*y + 4*y^2 + 2) + sqrt(7) + x^3 + 6*x^2*y + 12*x*y^2 + 8*y^3;
         41*x - 18*y - sqrt(5) + 8*x^3 - 12*x^2*y + 6*x*y^2 - y^3 + 3*sqrt(7)*(4*x*y - 4*x^2 - y^2 - 2)];

    % first iteration -- Jacobian has 1 dim'l null space at Pt = [5^(1/2)/5 + (2*7^(1/2))/5, (2*5^(1/2))/5 - 7^(1/2)/5];
    J = jacobian(F,X);
    A = J(1,1);
    B = J(1,2);
    N = det(A)*[-inv(A)*B;eye(1,1)];

    % deflate
    Fd = [F;J(2,:)*N];

    % second iteration -- Jacobian still has 1 dim'l null space
    Jd = jacobian(Fd,X);

    % deflate
    Fdd = [Fd;Jd(2:end,:)*N];

    % third iteration -- Jacobian still has 1 dim'l null space
    Jdd = jacobian(Fdd,X);    
    
    % deflate
    Fddd = [Fdd;Jdd(2:end,:)*N];

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);    
    
    % deflate
    Fdddd = [Fddd;Jddd(2:end,:)*N];
    
    % jacobian is full rank -- find distinct nonzero polynomials
    FD = unique(Fdddd);
    FD = FD(subs(FD,X,rand(size(X))) ~= 0);
    VD = X;
    
return;
