/*************************************************************
 Implementation of the QFSE (aka IP1S) algorithm described in
 "Practical Cryptanalysis of the Identification Scheme Based
    on the Isomorphism of Polynomial with One Secret Problem",
 Charles Bouillaguet, Jean-Charles Faugère, Pierre-Alain Fouque, Ludovic Perret,
 Public Key Cryptography 2011: 473-493,
 http://dx.doi.org/10.1007/978-3-642-19379-8_29

 This program practically breaks q=2, n=128 (using a lot of RAM)

 Author: Charles Bouillaguet <charles.bouillaguet@ens.fr>
 Please send me an email if you use this file.

 This work is licensed under the Creative Commons Attribution 3.0
 Unported License. To view a copy of this license, visit
 http://creativecommons.org/licenses/by/3.0/ or send a letter to
 Creative Commons, 444 Castro Street, Suite 900, Mountain View,
 California, 94041, USA.
*************************************************************/


// ---------- generates a random instance of the problem ------------
function RandomInstance(K, N, U)
    S := Random(GL(N, K));
    a := [ UpperTriangularMatrix([ Random(K) : i in [1..N*(N+1) div 2]]) : j in [1..U]];
    b := [ S*M*Transpose(S) : M in a];
    for k := 1 to U do for i := 1 to N do for j := i+1 to N do
	b[k][i,j] +:= b[k][j,i];
	b[k][j,i] := 0;
    end for;end for; end for;
    return S,a,b;
end function;


// ---------------- solver algorithm --------------------
function QFSESolver(a,b)
    K := BaseRing(a[1]);
    N := NumberOfRows(a[1]);

    assert #a eq #b;
    U := #a;

    if IsVerbose("User1") then 
      printf "linear algebra (computation of Ker M_phi)\n"; 
    end if;

    // compute the polar forms
    Pol_a := [ M + Transpose(M) : M in a];
    Pol_b := [ M + Transpose(M) : M in b];
    
    // build the matrix M_phi 
    In := IdentityMatrix(K,N);
    M_phi := Transpose(BlockMatrix([ [KroneckerProduct(In, Pol_a[i]), KroneckerProduct(Pol_b[i], -1*In)] : i in [1..U]]));

    // finds its kernel
    V := Kernel(M_phi);
    l := Dimension(V);

    if IsVerbose("User1") then 
      printf "Kernel has dimension %o \n", l;  
    end if;

    S_basis := [ Matrix(K, N, ElementToSequence(V.i)[1..N^2]) : i in [1..l]];
    S_basis_inv := [ Transpose(Matrix(K, N, ElementToSequence(V.i)[1+N^2..2*N^2])) : i in [1..l]];

     if IsVerbose("User1") then	
       printf "Generating the quadratic equations\n"; 
     end if;

    // most of the work should in fact be here
    M_rt := [ [ [S_basis[t] * a[i] * Transpose(S_basis[r]) : t in [1..l]] : r in [1..l]] : i in [1..U]];

    R_bot := PolynomialRing(K,l, "grevlex");
    // now assemble the multivariate polynomials

   Squad_uu := [
	b[i][u,u]
	- &+ [  M_rt[i][r][r][u,u] * R_bot.r * R_bot.r : r in [1..l]]
	- &+ [  (M_rt[i][r][t][u,u] + M_rt[i][t][r][u,u]) * R_bot.r * R_bot.t : r,t in [1..l] | r lt t]
	: u in [1..N], i in [1..U]];

    Squad_uv := &cat [
	[
	b[i][u,v]
	- &+ [  (M_rt[i][r][r][u,v] + M_rt[i][r][r][v,u] ) * R_bot.r * R_bot.r : r in [1..l]]
	- &+ [  (M_rt[i][r][t][u,v] + M_rt[i][r][t][v,u] + M_rt[i][t][r][u,v] + M_rt[i][t][r][v,u] ) * R_bot.r * R_bot.t : r,t in [1..l] | r lt t]
	: u,v in [1..N] | u lt v]
	: i in [1..U]];
    
    eqs := Squad_uu cat Squad_uv;    

    if IsVerbose("User1") then
      printf "Going Gröbner\n"; 
    end if;
        
    // small field, use field equations
    if (#K le 5) then 
	if IsVerbose("User1") then printf "small field ---> adding field equations\n"; end if;
	eqs := eqs cat [R_bot.i^(#K) - R_bot.i : i in [1..Rank(R_bot)]];
    end if;

    if (#K le 3) or (l le 16) then
	J := Ideal(eqs);
	Jlex := ChangeOrder(J, "lex");

	if IsVerbose("User1") then
	    printf "GRevLex GB computation....\n";
	    time Groebner(J);
	    printf "GRevLex --> Lex Order Conversion....\n";
	    time Groebner(Jlex : Al := "Walk");
	else
	    Groebner(J);
	    Groebner(Jlex : Al := "Walk");
	end if;

	Solutions := [ &+[ v[i] * S_basis[i] : i in [1..l]] : v in VarietySequence(J)];

    else
	// too big for direct GB. hybrid solving
	Solutions := [];
	i := 0;
	if IsVerbose("User1") then printf "problem too big ---> hybrid solving (there will be %o iterations)\n", #K; end if;
	for x in K do
	    i +:= 1;
	    if IsVerbose("User1") then printf "iteration %o / %o\n", i, #K; end if;
	    J := Ideal(eqs cat [R_bot.1 - x]);
	    this_it_sols := [ &+[ v[i] * S_basis[i] : i in [1..l]] : v in VarietySequence(J)];
	    if not IsEmpty(this_it_sols) and IsVerbose("User1") then printf "%o solution found\n", #this_it_sols; end if;
	    Solutions := Solutions cat this_it_sols;
	end for;
    end if;


    // Ouf, GB done
    // check solutions
    for Sol in Solutions do
	for i := 1 to U do
	    check := Sol*a[i]*Transpose(Sol) - b[i];
	    assert IsZero(check+Transpose(check));
	end for;
    end for;

    if IsVerbose("User1") then printf "%o solutions found and tested\n", #Solutions; end if;
    return Solutions;
end function;


// --------------------- demo script ----------------------------

// these parameters correspond to the biggest Challenge proposed by Patarin in 1996

N := 16;  // number of variables
U := 2;   // number of equations
q := 2;  
K:=FiniteField(q);


SetVerbose("User1", 1);
SetVerbose("Groebner", 0);

    printf "N = %o\n", N;
// Patarin's biggest Challenge
S,a,b := RandomInstance(K,N,U);
result := QFSESolver(a,b);
assert S in result;