/*************************************************************
 Implementation of practical Key-Recovery against SFLASH
 first described in the PhD thesis of Gilles Macario-Rat

 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.
*************************************************************/


// parameters for SFLASH v3
q := 128;
K := GF(q);
n_public := 3;
n := 67;
theta := 33;
assert GCD(1+q^theta,q^n - 1) eq 1;

printf "q = %o, n=%o\n", q, n;

L := ext<K | n>;
R<lambda,mu> := PolynomialRing(K,2);
R2<X> := PolynomialRing(R);
RL<lambdaL,muL> := PolynomialRing(L,2);
VV := VectorSpace(K,n^2);

V,toV := VectorSpace(L, K);
toL := toV^-1;
E := Basis(V);
SetVerbose("Groebner", 0);

TraceBase := Matrix(n, [Trace(toL(e_i)*toL(e_j),K) : e_i,e_j in E])^-1;

// Generates a C^* public key
function PKGen(S,T,n_pub)
    quad := [ ZeroMatrix(K,n,n) : i in [1..n_pub]];
    for i := 1 to n do
	for j := 1 to n do
	    plop := toV (toL(E[i]*S)*toL(E[j]*S)^(q^theta) )*T;
	    for k := 1 to n_pub do
		quad[k][i,j] := plop[k];
	    end for;
	end for;
    end for;
    return quad;
end function;
    
printf "Key Generation\n";
time PK := PKGen(Random(GL(n,K)),Random(GL(n,K)), n_public);
// test the public key
//x := Random(V);
//assert toV( toL(x*S)^(1 + q^theta))*T eq Vector([ (x*M,x) : M in PK ]);

// ---------------------- Attack Starts Here -----------------------

delta := ((q div 2) - 1)*&+[q^i : i in [0..n-1]] + &+ [ q^(2*i*theta) : i in [(n+1) div 2..n-1]];
d := (n-1) div 2; 

function WorkOutPencil(Q1, Q2)
    phi_P := lambda*(Q1 + Transpose(Q1)) + mu*(Q2 + Transpose(Q2));   
    // Option 1, via the characteristic polynomial [slower]
    //    time Chi := CharacteristicPolynomial(phi_P);  
    //    p_i := [ Sqrt(MonomialCoefficient(R2!Chi, X^(n-2*i))) : i in [0..d]];
    //    Adj_square_root := &+ [ p_i[i+1] * phi_P^(d-i): i in [0..d]]; 
    //    k_i := Vector([ &+ [Adj_square_root[i,j] : j in [1..n]]: i in [1..n]]); 
    // or
    // Option 2, using an implicit Gröbner basis computation [much faster]
    k_i := Vector(R, ElementToSequence(Kernel(phi_P).1));
    P_S := &+ [ RL!k_i[i] * toL(E[i]) : i in [1..n]];  
    P_N := (k_i*(lambda*Q1 + mu*Q2),k_i);
    return P_S, [ toV(x[1]) : x in Roots(UnivariatePolynomial(Evaluate(RL!P_N,muL,1)))]; 
end function;

// --------- main loop --------------

"Working out pencils";
time P_S12, roots_12 := WorkOutPencil(PK[1], PK[2]);
time P_S13, roots_13 := WorkOutPencil(PK[1], PK[3]);

my_T1 := toV(1);
repeat my_T2 := toV(toL(Random(roots_12))*toL(my_T1*TraceBase))*TraceBase^-1; until Rank(Matrix([my_T1,my_T2])) eq 2;

reconstructed_P_S12 := &* [ lambdaL*toL(my_T1*TraceBase)^(q^(2*i*theta))+muL*toL(my_T2*TraceBase)^(q^(2*i*theta)) : i in [(n+1) div 2..n-1]];
Sx := [ toV(x) : x in Coefficients(P_S12)];
Sy := [ toV(x) : x in Coefficients(reconstructed_P_S12)];

printf "Now trying the possible T_3\n";    
time for choice in roots_13 do
    my_T3 := toV(toL(choice)*toL(my_T1*TraceBase))*TraceBase^-1;
  time if 1 eq 1 then
  reconstructed_PS13 := &* [ lambdaL*toL(my_T1*TraceBase)^(q^(2*i*theta))+muL*toL(my_T3*TraceBase)^(q^(2*i*theta)) : i in [(n+1) div 2..n-1]];
  Sx2 := [ toV(x) : x in Coefficients(P_S13)];
  Sy2 := [ toV(x) : x in Coefficients(reconstructed_PS13)];

    full_Sx := Sx cat Sx2;
    full_Sy := Sy cat Sy2;
  
  full_Sx := Sx cat Sx2;// cat Sx3;
  full_Sy := Sy cat Sy2;// cat Sy3;
  
    // find a suitable basis of K^n over which S is known
  if (Rank(Matrix(full_Sx)) lt n) then
    printf "< : %o  ", Rank(Matrix(full_Sx));
  else
    if my_T3 in sub<V | [my_T1, my_T2]> then
      printf "@";
    else
  repeat
      potential_basis_idx := {};
	repeat Include(~potential_basis_idx, Random([1..#full_Sx])); until #potential_basis_idx eq n;
      potential_basis := [ full_Sx[i] : i in potential_basis_idx ];
  until IsInvertible(Matrix(potential_basis)); 
  
    // reconstruct S in the canonical basis
  possible_S := Matrix(potential_basis)^-1*Matrix([ full_Sy[i] : i in potential_basis_idx ]);
  
    // check if S is correct
    A_quad := PKGen(possible_S, IdentityMatrix(K,n),n);
  A_space := sub<VV | [ ElementToSequence(x) : x in A_quad]>;
    if &and [ Vector(ElementToSequence(M)) in A_space : M in PK ] then
	printf "SOLUTION FOUND !!!\n";
	break;
end if;
end if;
    end if;
end for;