//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

// Cost estimates for isogeny computation using Deuring for the people without Algo 5

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Tp := [
    [ 3, 38, 1 ],
    [ 7, 6, 1 ],
    [ 11, 4, 1 ],
    [ 439, 4, 1 ],    
    [ 347, 2, 1 ],
    [ 257, 2, 1 ],
    [ 137, 2, 1 ],
    [ 59, 2, 1 ],
    [ 5, 3, 1 ],
    [ 359, 1, 1 ],
    [ 353, 1, 1 ],
    [ 311, 1, 1 ],
    [ 281, 1, 1 ],
    [ 197, 1, 1 ],
    [ 173, 1, 1 ],
    [ 139, 1, 1 ],
    [ 89, 1, 1 ],
    [ 47, 1, 1 ],
    [ 31, 1, 1 ],
    [ 19, 1, 1 ],
    [ 13, 1, 1 ],
    [ 577, 1, 1 ],
    [ 733, 1, 1 ],
    [ 853, 1, 1 ],
    [ 983, 1, 1 ],
    [ 1223, 1, 1 ]
    ];


// Mult Cost per extension for k = 1, ... 
Mk := [1, 3, 6, 9, 13, 18, 22, 26, 30, 35, 45, 52, 58, 66, 74, 81, 90, 98, 107, 116, 132, 135, 144, 156, 169, 174, 180, 198,\
208, 210, 232, 243, 270, 270, 286, 294, 306, 321, 348, 348, 360, 396];

// Sqr Cost per extension for k = 1, ... 
Sk := [1, 2, 5, 6, 12, 12, 22, 12, 25, 24, 24, 24, 58, 44, 60, 24, 50, 50, 55, 48, 110, 48, 80, 48, 144, 116, 125, 88, 110, 120,\
120, 48, 120, 100, 150, 100, 160, 170, 110, 96, 180, 220];


// up to 30 extensions. Bound for powers of two to be considered for the ell_i per extension k
// For example for k = 1 we allow ell_i's of up to 16 bits. For k = 2 ell_i's of 13 bits and so on. 
BoundT2 := [16, 13, 10, 13, 10, 10, 10, 12, 11, 10, 10, 10, 10, 10, 10, 16, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 12];

// **************************

// Auxiliary Functions

// **************************


Switch_KPS_Column := function(kps, i, j)
N := #kps;
Temp := [];

for k:= 1 to N do
    Temp[k] := kps[k];
    t:= kps[k][i];
    Temp[k][i] := kps[k][j];
    Temp[k][j] := t;
end for;

return Temp;
end function;




Sort_T := function(kps) // Sort for multiplicity
   // First we find how many multiciplities we have 
    N := #kps; //Total number of \ell_i's to be processed
    Points_multi := []; Ext_A := [];
    ext := 0; Nm := 0;

    while(Nm lt N) do
        for i:= 1 to N do
            if(kps[i][3] notin Ext_A) then
                ext := kps[i][3];
                Append(~Ext_A, ext);
                break;
            end if;
        end for;
        Points_Aux := [];
        for i := 1 to N do
            if(kps[i][3] eq ext) then
                Append(~Points_Aux, kps[i]);
            end if;
        end for; 
        Points_Aux := Sort(Switch_KPS_Column(Points_Aux, 1, 2));
        Points_Aux :=  Switch_KPS_Column(Points_Aux, 1, 2);
        Insert(~Points_multi, Nm+1, Nm+1, Points_Aux);
        Nm := #Points_multi;
    end while;    
    return Points_multi;
end function;


Sort_kps := function(kps)
    N := #kps;
    kps_odd := []; kps_even := [];
    for i:= 1 to N do
        if (kps[i][3] in [1, 2, 4, 8]) then
            Append(~kps_even, kps[i]);
        else    
            Append(~kps_odd, kps[i]);
        end if;
    end for;
    //kps := Insert(Sort(kps_odd), #kps_odd+1, #kps_odd+1, Sort(kps_even));
    kps := Insert(kps_odd, #kps_odd+1, #kps_odd+1, kps_even);
    return Sort_T(kps);
end function;



Fill_kps_csidh := function(kps, S, k)

    N := #S;
    aux := [];
    for i := 1 to N do
        m := S[i][2];    
        Append(~aux, [S[i][1], m, k]);
    end for; 
    Insert(~kps, #kps +1, #kps +1, Sort(aux));
    return kps;
end function;



// Function that collects T-Torsion
CollectT := function(p, Bound, limit, verbose)
AddT :=[];
kps := [];
Tbits := 0;
kfin := 100; // acting as infinite
for k := 1 to limit do
    S := TrialDivision(p^(2*k) -1, 2^BoundT2[k]);
    //S := TrialDivision(p^(2*k) -1, 2^16);
    Nini := #AddT;
    DiffT := [];
    for i := 2 to #S do
        ell := S[i][1]^S[i][2];
        if(S[i][1] gt 2^BoundT2[k]) then
        //if(S[i][1] gt 2^12) then
            continue;
        end if;
        N := #AddT;
        flag := 1;
        for j := 1 to N do
            if ((AddT[j] mod S[i][1] eq 0)) then
                flag := 0;
                break;
            end if;
        end for;
        if (flag eq 1) then
            Append(~AddT, ell);
            Append(~DiffT, ell);
            //ArithCost +:= CosT(ell, k);
        end if;
    end for;
    Tbits := Log(2, &*AddT);
    if (Nini lt #AddT) then
                    if(verbose eq true) then
                        printf "k =%o:  %o\n", k, Factorization(&*DiffT);
                        printf "Tbits = %o\n", Tbits;
                    end if;
        // Here we generate kps with information about the extensions and the \ell_is
        kps := Fill_kps_csidh(kps, Factorization(&*DiffT), k);
    end if;
    if (Tbits ge 320) then
        kfin := k;
        break;
    end if;
end for;
T := Log(2, &*AddT);

return Sort_kps(kps), Factorization(&*AddT), kfin, Log(2, &*AddT);
end function;



//***************************************************************** //
// Functions from the Deuring for the people stuff adaptation
//***************************************************************** //


// cost of computing minimal polynomial using Shoup
cost_Shoup := function(k)
	Mp2k := Mk[k];
	cS := (4 * Ceiling(Sqrt(2*k)) + 1) * Mp2k + 4*k^2 - k + 1;
	return cS;
end function;

// Cost of computing kernel polynomial and codomain curve
cost_xisog_dftp := function(l, k)
    m := (l-1) div (2*k);
    primroot_l := 1;
    while not IsPrimitive(primroot_l,l) do
        primroot_l := primroot_l+1;
    end while;
    size_primitivel := Floor(Log(2, primroot_l));
    c_xMul := (m-1) * size_primitivel * (8 * Mk[k] + 2*Sk[k]);
    if primroot_l eq 2 then
        c_xMul := (m-1) * 4 * Mk[k];
    end if;
    c_Shoup := m * cost_Shoup(k);
    expKrtsb := Log(2, 3);
    c_prod := Ceiling((m*k)^expKrtsb - m*k^expKrtsb);
    c_xIsog := 4*Floor(Log(2, l)) + 8;
    cost := c_xMul + c_Shoup + c_prod + c_xIsog;

    return Ceiling(cost);
end function;


// Cost of evaluating points by pushing x-coordinates
cost_xeval_dftp := function(l, k)
	lp := (l+1) div 2;
	kr := Floor(Sqrt(lp));
    c_xEval := (4 * (kr-1) + 2) * Mk[k] + (l-1) * k;

    return Ceiling(c_xEval);
end function;

// Cost of computing an \ell isogeny
cost_xisog_velu := function(l, k)
    d := (l -1) div 2;
    kps := 4*(d-1)*Mk[k] + 2*(d-1)*Sk[k];
    xisog := 1*l*Mk[k] + (2*Log(2, l) +4)*Sk[k];
    velu := kps + xisog;

    return Ceiling(velu);
end function;


// Cost of evaluating an \ell isogeny
cost_xeval_velu := function(l, k)
    d := (l -1) div 2;
    xeval := 4*d*Mk[k] + 2*Sk[k];

    return Ceiling(xeval);
end function;


// Cost of computing an \ell isogeny
cost_xisog_vsqt := function(l, k)
    // TODO: Check Jorge & Odalis paper for # of Mults and Sqrs
    b:= Ceiling(Sqrt(l-1) / 2);
    kps := (4*b^Log(2, 3) + 2*Log(2, b) +16*b -2) * Mk[k];
    xisog := (18*b^Log(2, 3) + 6*Log(2, b) -5*b +12) * Mk[k];
    vsqt := kps + xisog;

    return Ceiling(vsqt);
end function;


// Cost of evaluating an \ell isogeny
cost_xeval_vsqt := function(l, k)
    // TODO:  Check Jorge & Odalis paper for # of Mults and Sqrs
    b:= Ceiling(Sqrt(l-1) / 2);
    xeval := 15*b^Log(2, 3) + 5*b +2;
    xeval := Mk[k]*xeval;

    return Ceiling(xeval);
end function;


// Total cost with CSIDH strategies
// LMK is a sequence composed of objects [l, m, k]
// l: isogeny degree, m: multiplicity of the isogeny, k: extension degree
cost_total_csidh := function(LMK, K, verbose)
    Cost_xisog := 0; Cost_xeval := 0; Cost_xmul := 0;
    N_xisog := 0; N_xeval := 0; N_xmul := 0;

    LMK_temp := LMK;
    for lmk in LMK do
        Cost_xisog_dftp := 0; Cost_xisog_velu := 0; Cost_xisog_vsqt := 0;
        Cost_xeval_dftp := 0; Cost_xeval_velu := 0; Cost_xeval_vsqt := 0;

        l,m,k := Explode(lmk);

        if verbose then
            printf "processing point %o\n", lmk; 
        end if;

        // Cost of computing scalar multiplications a la CSIDH
        next_k := 0;
        for i := 2 to #LMK_temp do
            lp,mp,kp := Explode(LMK_temp[i]);
            if (k/kp in [2,4,8]) and (kp in K) then
                break;
            elif k/kp in [1,2,4,8] then
                Cost_xmul +:= mp*Ceiling(Log(2,lp))*(12 * Mk[k] + 6*Sk[k]);
                N_xmul +:= mp;
                next_k := Max(next_k, kp);
            end if;
        end for;

        // Cost of finding m image curves
        Cost_xisog_dftp +:= m*cost_xisog_dftp(l, k);
        Cost_xisog_velu +:= m*cost_xisog_velu(l, k);
        Cost_xisog_vsqt +:= m*cost_xisog_vsqt(l, k);
        N_xisog +:= m;
        if m gt 1 then
            Cost_xmul +:= m*Log(2, m)*Ceiling(Log(2,l))*(12 * Mk[k] + 6*Sk[k]);
            N_xmul +:= Ceiling(m*Log(2, m));

            Cost_xeval_dftp +:= m*Log(2, m)*cost_xeval_dftp(l, k);
            Cost_xeval_velu +:= m*Log(2, m)*cost_xeval_velu(l, k);
            Cost_xeval_vsqt +:= m*Log(2, m)*cost_xeval_vsqt(l, k);
            N_xeval  +:= Ceiling(m*Log(2, m));
        end if;


        Cost_xmul_dftp := 0; Cost_xmul_velu := 0; Cost_xmul_vsqt := 0;
        N_xmul_dftp := 0; N_xmul_velu := 0; N_xmul_vsqt := 0;
        N_xeval_dftp := 0; N_xeval_velu := 0; N_xeval_vsqt := 0;
        ind := 0;

        if next_k notin {0,k} then
            Cost_xmul_sub := Ceiling(Log(2,l))*(12 * Mk[k] + 6*Sk[k]);
            Cost_xeval_sub_dftp := cost_xeval_dftp(l, k);
            Cost_xeval_sub_velu := cost_xeval_velu(l, k);
            Cost_xeval_sub_vsqt := cost_xeval_vsqt(l, k);

            if Cost_xmul_sub le Cost_xeval_sub_dftp then
                Cost_xmul_dftp := Cost_xmul_sub;
                N_xmul_dftp := 1;
            else
                Cost_xeval_dftp +:= Cost_xeval_sub_dftp;
                N_xeval_dftp := 1;
            end if;

            if Cost_xmul_sub le Cost_xeval_sub_velu then
                Cost_xmul_velu := Cost_xmul_sub;
                N_xmul_velu := 1;                
            else
                Cost_xeval_velu +:= Cost_xeval_sub_velu;
                N_xeval_velu := 1;
            end if;

            if Cost_xmul_sub le Cost_xeval_sub_vsqt then
                Cost_xmul_vsqt := Cost_xmul_sub;
                N_xmul_vsqt := 1;
            else
                Cost_xeval_vsqt +:= Cost_xeval_sub_vsqt;
                N_xeval_vsqt := 1;
            end if;

            ind := Index(K, k);
            K[ind] := next_k;
        elif next_k eq 0 then
            Exclude(~K, k); // No need to push the full torsion point in the same ext for the final \ell_i
        end if;

        for kp in K do // pushing all required points
            if (Index(K, kp) ne ind) or (next_k in {0,k}) then
                Cost_xeval_dftp +:= m*cost_xeval_dftp(l, kp);   
                Cost_xeval_velu +:= m*cost_xeval_velu(l, k);
                Cost_xeval_vsqt +:= m*cost_xeval_vsqt(l, k);   
                N_xeval +:= m;
            end if;
        end for;

        Cost_dftp := Cost_xisog_dftp + Cost_xeval_dftp + Cost_xmul_dftp;
        Cost_velu := Cost_xisog_velu + Cost_xeval_velu + Cost_xmul_velu;
        Cost_vsqt := Cost_xisog_vsqt + Cost_xeval_vsqt + Cost_xmul_vsqt;

        min_cost := Min([Cost_dftp, Cost_velu, Cost_vsqt]);

        if (k notin [1,2,4,8]) or (Cost_dftp eq min_cost) then
            Cost_xisog +:= Cost_xisog_dftp;
            Cost_xeval +:= Cost_xeval_dftp;
            Cost_xmul +:= Cost_xmul_dftp;
            N_xeval +:= N_xeval_dftp;
            N_xmul +:= N_xmul_dftp;
        elif Cost_velu eq min_cost then
            Cost_xisog +:= Cost_xisog_velu;
            Cost_xeval +:= Cost_xeval_velu;
            Cost_xmul +:= Cost_xmul_velu;
            N_xeval +:= N_xeval_velu;
            N_xmul +:= N_xmul_velu;
        else
            Cost_xisog +:= Cost_xisog_vsqt;
            Cost_xeval +:= Cost_xeval_vsqt;
            Cost_xmul +:= Cost_xmul_vsqt;
            N_xeval +:= N_xeval_vsqt;
            N_xmul +:= N_xmul_vsqt;
        end if;

        Remove(~LMK_temp, 1); // Removing the \ell_i processed from LMK
    end for;
    
    //printf "\n# of x_isog = %o, # of x_eval = %o, # of xmul = %o", N_xisog, N_xeval, N_xmul;
    //printf "\ncost x_isog = %o, cost x_eval = %o, cost xmul = %o\n", Ceiling(Cost_xisog), Ceiling(Cost_xeval), Ceiling(Cost_xmul);
    return Ceiling(Cost_xeval + Cost_xisog + Cost_xmul);
end function;


pow2 := [1, 2, 4, 8];

sort_LMK := function(LMK, K)
    N := #LMK;
    K_odd := []; K_even := [];
    LMK_odd := []; LMK_even := [];
    for lmk in LMK do
        k := lmk[3];
        if k notin (pow2 cat K_odd) then
            Append(~K_odd, k);
            Append(~LMK_odd, [lmk]);
        elif k in K_odd then
            ind := Index(K_odd, k);
            Append(~LMK_odd[ind], lmk);
        elif k notin K_even then
            Append(~K_even, k);
            Append(~LMK_even, [lmk]);
        else
            ind := Index(K_even, k);
            Append(~LMK_even[ind], lmk);
        end if;
    end for;

    for ind := 1 to #LMK_odd do
        Sort(~LMK_odd[ind]);
        Reverse(~LMK_odd[ind]);

        LMK_k := LMK_odd[ind];
        Cost_big_last := cost_total_csidh(LMK_k, [K_odd[ind]], false);

        Append(~LMK_k, LMK_k[1]);
        Remove(~LMK_k, 1);
        Cost_big_first := cost_total_csidh(LMK_k, [K_odd[ind]], false);

        if Cost_big_first lt Cost_big_last then
            LMK_odd[ind] := LMK_k;
        end if;
    end for;


    for ind := 1 to #K_even do
        Sort(~LMK_even[ind]);
        Reverse(~LMK_even[ind]);

        LMK_k := LMK_even[ind];
        Cost_big_last := cost_total_csidh(LMK_k, [K_even[ind]], false);

        Append(~LMK_k, LMK_k[1]);
        Remove(~LMK_k, 1);
        Cost_big_first := cost_total_csidh(LMK_k, [K_even[ind]], false);

        if Cost_big_first lt Cost_big_last then
            LMK_even[ind] := LMK_k;
        end if;
    end for;

    Sort(~LMK_even, func<LMK_k_i, LMK_k_j | LMK_k_i[1][3] - LMK_k_j[1][3]>);
    Reverse(~LMK_even);
    LMK_even := &cat LMK_even;

    Cost_best := 2^25;

    Permuts_K_odd :=  Permutations(Seqset(K_odd));

    for Permut_k in Permuts_K_odd do
        ind3 := Index(Permut_k, 3);
        ind6 := Index(Permut_k, 6);
        if ind3*ind6 ne 0 and ind3 le ind6 then
            continue;
        end if;

        LMK_permut := [];
        for k in Permut_k do
            ind := Index(K_odd, k);
            Append(~LMK_permut, LMK_odd[ind]);
        end for;


        LMK_permut := &cat LMK_permut;
        LMK_permut := LMK_permut cat LMK_even;
        Cost_permut := cost_total_csidh(LMK_permut, K, false);

        if Cost_permut lt Cost_best then
            Cost_best := Cost_permut;
            LMK_best := LMK_permut;
        end if;
    end for;

    return Cost_best, LMK_best;
end function;



// **************************

// Main

// **************************
f:= 80;
Bound := 2^14;
expi := 3;
nprimes := 8;
delta := 4;
limit := 10;
SMax := 0;
MinSbits := 320;


// T torsion of p_1973, the official prime for SQISign level 1.
// The performance of this prime is our baseline
printf "**************Scoring p_1973 from SQISign Level 1************\n";
T := 3^38*5^3*7^6*11^4*13*19*31* 47 * 59^2 *89 * 137^2 * 139 * 173 * 197 * 257^2 *281 * 311 * 347^2* 353 * 359 * 439^4 *577 * 733 * 853 * 983 * 1223;
kps_T := []; kps_T := Sort_T(Fill_kps_csidh(kps_T, Factorization(T), 1));

//Cost_p1973, LMK_p1973 :=sort_LMK(kps_T, [1]);
//printf "\nLMK_p10 = %o", LMK_p1973;
Cost_p1973 := cost_total_csidh(Tp, [1], false);
printf "\ncost SQISign prime level 1 p_1973 = %o", Cost_p1973;


// Scoring p7 from ApresSQI
printf "**************Scoring p7 from Apres************\n";
p7 := 2^145*3^9*59^3*311^3*317^3*503^3-1;
kps_A, S, kf, bits := CollectT(p7, Bound, limit, true);
Cost_p7apres, LMK_p7apres :=sort_LMK(kps_A, [1,2,3,4,5,6,7,8]);
printf "\nell_i array for ApresSQI prime p7 = %o", LMK_p7apres;
printf "\ncost ApresSQI prime p7 = %o", Cost_p7apres;


// Scoring a new p7 prime
//p7new :=  2^155 * 3^3 *13^3 * 17^3 * 59^3 *173^3 * 317 * 1567*1609 -1; 
p7new := 2^149 * 7^3 * 17^3 * 23^3 * 61^3*157^3 * 823 * 1051 * 1297 -1;
kps_A, S, kf, bits := CollectT(p7new, Bound, limit, true);
Cost_p7new, LMK_p7new :=sort_LMK(kps_A, [1,3, 4,6,7]);
printf "\nell_i array for new prime p7 = %o", LMK_p7new;
printf "\ncost new prime p7 = %o", Cost_p7new;



// Scoring a new p10
//p := 0x2CB080F3AD9271CA97956750695AA8898F3E024FFFFFFFFFFFFFFFFFFFFFFFFF;
p10 := 2^100 * 5^16 *13^8 *127^8 * 433 * 2069 * 2113 -1; 
kps_A, S, kf, bits := CollectT(p10, Bound, limit, true);
Cost_p10, LMK_p10 :=sort_LMK(kps_A, [1,2,3,5]);
printf "\nell_i array for new prime p10 = %o", LMK_p10;
printf "\ncost new prime p10 = %o", Cost_p10;



p_3917 := 24646990764040422654033369043417789695753066731129803480569060745362536398847;
p := p_3917;
kps_A, S, kf, bits := CollectT(p, Bound, limit, true);
Cost_p10, LMK_p10 :=sort_LMK(kps_A, [1,3,4,6]);
e := Ceiling((15/4)*Log(2, p)) + 25; Ceiling(Log(2, p));
f := Factorization(p+1)[1][2];
printf "\nell_i array for new prime p3917 = %o", LMK_p10;
printf "\ncost new prime p3917 = %o", Ceiling(e/f)*Cost_p10;


flag := false;

p1973 := 0x34e29e286b95d98c33a6a86587407437252c9e49355147ffffffffffffffffff;
p := p1973;
kps_A, S, kf, bits := CollectT(p, Bound, limit, flag);
Cost_p10, LMK_p10 :=sort_LMK(kps_A, [1]);
e := Ceiling((15/4)*Log(2, p)) + 25; Ceiling(Log(2, p));
f := Factorization(p+1)[1][2];
printf "\nell_i array for new prime p3917 = %o", LMK_p10;
printf "\ncost new prime p3917 = %o", Ceiling(e/f)*Cost_p10;


p1223 := 0xea6a4dda9518e5c5d50ccdfbd97e4c49efe85e0e09039c7ffffffffffffffff;
p := p1223;
kps_A, S, kf, bits := CollectT(p, Bound, limit, flag);
Cost_p10, LMK_p10 :=sort_LMK(kps_A, [1]);
e := Ceiling((15/4)*Log(2, p)) + 25; Ceiling(Log(2, p));
f := Factorization(p+1)[1][2];
printf "\nell_i array for new prime p1223 = %o", LMK_p10;
printf "\ncost new prime p3917 = %o", Ceiling(e/f)*Cost_p10;


p8011 := 0x31ebc32c245c72c40115748f25c4ba516cb58aaae247ffffffffffffffffffff;
p := p8011;
kps_A, S, kf, bits := CollectT(p, Bound, limit, flag);
Cost_p10, LMK_p10 :=sort_LMK(kps_A, [1]);
e := Ceiling((15/4)*Log(2, p)) + 25; Ceiling(Log(2, p));
f := Factorization(p+1)[1][2];
printf "\nell_i array for new prime p8011 = %o", LMK_p10;
printf "\ncost new prime p3917 = %o", Ceiling(e/f)*Cost_p10;