// @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ PROJECT D @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
// Setup a polynomial ring

W<X,RP,RP_,P, P_> := PolynomialRing(Rationals(),5);
K := quo<W|RP*RP_-1,RP^2-P,RP_^2-P_, P_*P-1 >;

/////////////////////////////////////////////////////////////////////////////////////////
// Function List
//  L_k := function(k,l)
//  l_k := function(k,l)
//  d_k := function(k,l)
//  v_k := function(k,e,l,p)
//  delta_k := function(k,l,sign)
//  f_1 := function(alpha,pow,l,p)
//  f_2 := function(alpha,pow,l,p)
//  g_k := function(k,l,a1,p,a)
//  I_1_1 := function(e,l,p,a)
//  I_1_i_aeqb := function(i,e,l,a1,a2,p,a)
//  I_1_i_altb := function(i,e,l,a1,a2,p,a,b)
//  I_2_1 := function(e,l,p,a)
//  I_2_i_aeqb := function(i,e,l,a1,a2,p,a)
//  I_2_i_altb := function(i,e,l,a1,a2,p,a,b)
//  R := function(n,e,l,a1,a2,p,a,b)
//  beta_l := function(e,l,a1,a2,p,a,b,l_val);

//  MAIN FUNCTION:
//      alpha_X := function(e,l,a1,a2,p,a,b,l_val);
// ////////////// //////////////////////// /////////////////////// ///////////////////////////

// Power of (sqrt p)
POP :=function(k)
     if k ge 0 then return RP^k;
     else   return RP_^(-1*k);
     end if;
end function;

// L(k,1) = {1 <= i <= m : l_i - k < 0 is odd}
L_k := function(k,l)
    count := 1;
    L := [];
    for i := 1 to #(l) do
        dif := l[i] - k;
        if dif lt 0 and IsOdd(dif) then
            L[count] := i;
            count +:= 1;
        end if;
    end for;
    return L;
end function;

/////////////////////////////////////////////////////////////////////////////////////////

// l(k,1) = #L(k,1)
l_k := func<k,l|#(L_k(k,l))>;

/////////////////////////////////////////////////////////////////////////////////////////

// d(k) = k + 1/2 SUM_{l_i < k} (l_i - k)
d_k := function(k,l)
    val_d := 0;
    for i := 1 to #(l) do
        if l[i] lt k then val_d +:= l[i]-k;
        end if;
    end for;
    // As the value of some polynomial power, d(k) should be return as an "integer".
    // For this program's sake, d(k) is a "1/2 - base" integer.
    // Ex: if d(k) = 3/2, then its return value is 3; if d(k) = 4, then its return
    //  value is 8.

    return (Integers() ! (2*k + val_d));
end function;

/////////////////////////////////////////////////////////////////////////////////////////

// v(k) = (-1|p)^[l(k,1)/2]*PROD_{i in L(k,1)} (e_i|p)
v_k := function(k,e,l,p)
    kron := KroneckerSymbol(-1,p);
    l_k_val := l_k(k,l);
    min_int := 1;
    if IsEven(l_k_val) then min_int := l_k_val/2;
    else min_int := (l_k_val-1)/2;
    end if;

    prd := 1;
    L := L_k(k,l);
    for idx := 1 to #(L) do
        i := L[idx];
        prd *:= KroneckerSymbol(e[i],p);
    end for;
    return kron^(Integers() ! min_int)*prd;
end function;

/////////////////////////////////////////////////////////////////////////////////////////

// delta(k)^(+-)
delta_k := function(k,l,sign)
    l_k_val := l_k(k,l);
    return (1+sign*((-1)^(l_k_val)))/2;
end function;

/////////////////////////////////////////////////////////////////////////////////////////

// f_1(alpha*p^a)
f_1 := function(alpha,pow,l,p)
    l_val := l_k(pow+1,l);
    if IsEven(l_val) then return -1/p;
    else
        kron := KroneckerSymbol(alpha,p);
        return kron*RP_;
    end if;
end function;

/////////////////////////////////////////////////////////////////////////////////////////

// f_2(alpha*p^a)
f_2 := function(alpha,pow,l,p)
    l_val := l_k(pow+1,l);
    if IsOdd(l_val) then return -1/p;
    else
        kron := KroneckerSymbol(-1*alpha,p);
        return kron*RP_;
    end if;
end function;

/////////////////////////////////////////////////////////////////////////////////////////
// g(k)  ???? a1?
g_k := function(k,l,a1,p,a)
    if IsEven(a-k) then return delta_k(k,l,1);
    else return KroneckerSymbol(a1,p)*delta_k(k,l,-1);
    end if;
end function;

/////////////////////////////////////////////////////////////////////////////////////////

// I_(1,1) when a = b or a < b
I_1_1 := function(e,l,p,a)
    sum := 0;
    for k := 1 to a do
        sum +:= v_k(k,e,l,p)*delta_k(k,l,1)*POP(2*k+d_k(k,l))*X^k;
    end for;
    return (1-p^(-2))*sum;
end function;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// case statement I_(1,i) when a = b
I_1_i_aeqb := function(i,e,l,a1,a2,p,a)
    case i:
        // I_(1,1) when a = b
        when 1:
            return I_1_1(e,l,p,a);
        // I_(1,2) when a = b
        when 2:
            return delta_k(1+a,l,1)*v_k(1+a,e,l,p)*(KroneckerSymbol(-1*a1*a2,p)-p^-1)*POP(2*a+d_k(a+1,l))*X^(a+1);
        // I_(1,3) = I_(1,4) = 0 when a = b
        when 3:
            return 0;
        when 4:
            return 0;
        else error "Error. For I_(1,i), i must be in {1..4}.";
    end case;
end function;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// I_(2,1) when a = b or a < b
I_2_1 := function(e,l,p,a)
    sum := 0;
    for k2 in {1..(a-1)} do
        for k1 in {(k2+1)..a} do
            sum +:= v_k(k1,e,l,p)*v_k(k2,e,l,p)*delta_k(k1,l,1)*delta_k(k2,l,1)
                *POP(2*k1+d_k(k1,l)+d_k(k2,l))*X^(k1+k2);
        end for;
    end for;
    return (1-p^(-2))*sum;
end function;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// case statement I_(2,i) when a = b
I_2_i_aeqb := function(i,e,l,a1,a2,p,a)
    case i:
        // I_(2,1) when a = b
        when 1:
            return I_2_1(e,l,p,a);
        // I_(2,2) when a = b
        when 2:
            sum := 0;
            for k := 1 to a do
                sum +:= delta_k(k,l,1)*v_k(k,e,l,p)*POP(2*a+d_k(1+a,l)+d_k(k,l))*X^(a+1+k);
            end for;
             kron :=KroneckerSymbol(-1*a1*a2,p);
            return v_k(1+a,e,l,p)*delta_k(a+1,l,1)*(kron-p^(-1))*sum;
        // I_(2,3) = I_(2,4) = I_(2,5) = 0
        when 3:
            return 0;
        when 4:
            return 0;
        when 5:
            return 0;
        // I_(2_6) when a = b
        when 6:
            kron :=KroneckerSymbol(-1*a1*a2,p);
            return -1*(kron*delta_k(1+a,l,1)+delta_k(1+a,l,-1))*POP(2*(a-1)+2*d_k(1+a,l))*X^(2*a+2);
        // I_(2,7) = 0
        when 7:
            return 0;
        // I_(2,8) when a = b
        when 8:
            sum := 0;
             kron1 :=KroneckerSymbol(-1,p);
            for k := 1 to a do
                sum +:= (delta_k(k,l,1)+p^(-1)*delta_k(k,l,-1))*POP(2*k+2*d_k(k,l))*X^(2*k);
            end for;
            return sum;
        else error "Error. For I_(2,i), i must be in {1..8}";
    end case;
end function;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// case statement I_(1,i) when a < b
I_1_i_altb := function(i,e,l,a1,a2,p,a,b)
    case i:
        // I_(1,1) when a < b
        when 1:
            return I_1_1(e,l,p,a);
        // I_(1,2) when a < b
        when 2:
            kron := KroneckerSymbol(a1,p);
            return v_k(a+1,e,l,p)*(kron*RP_*delta_k(a+1,l,-1)-p^-2*delta_k(a+1,l,1))*POP(2*(a+1)+d_k(a+1,l))*X^(a+1);
        // I_(1,3) when a < b
        when 3:
            sum := 0;
            for k := (a+2) to b do
                sum +:= v_k(k,e,l,p)*g_k(k,l,a1,p,a)*POP(a+k+d_k(k,l))*X^k;
            end for;
            return (1-p^-1)*sum;
        // I_(1,4) when a < b
        when 4:
            cont := 0;
            if ((a-b) mod 2) ne 0 then cont := f_1(a2,b,l,p);
            elif ((a-b) mod 2) eq 0 then cont := KroneckerSymbol(a1,p)*f_2(a2,b,l,p);
            end if;
            return v_k(b+1,e,l,p)*POP(a+b+1+d_k(b+1,l))*cont*X^(b+1);
        else error "Error. For I_(1,i) when a < b, i must be in {1..4}.";
    end case;
end function;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// case statement I_(2,i) when a < b
I_2_i_altb := function(i,e,l,a1,a2,p,a,b)
    case i:
        // I_(2,1) when a < b
        when 1:
            return I_2_1(e,l,p,a);
        // I_(2,2) when a < b
        when 2:
            kron := KroneckerSymbol(a1,p);
            sum := 0;
            for k := 1 to a do
                sum +:= v_k(a+1,e,l,p)*v_k(k,e,l,p)*delta_k(k,l,1)*POP(2*(a+1)+d_k(a+1,l)+d_k(k,l))*X^(a+1+k);
            end for;
            return (kron*RP_*delta_k(a+1,l,-1)-p^-2*delta_k(a+1,l,1))*sum;
        // I_(2,3) when a < b
        when 3:
            sum := 0;
            for k2 in {1..a} do
                for k1 in {(a+2)..b} do
                    sum +:= v_k(k1,e,l,p)*v_k(k2,e,l,p)*delta_k(k2,l,1)*g_k(k1,l,a1,p,a)*POP(a+k1+d_k(k1,l)+d_k(k2,l));
                end for;
            end for;
            return (1-p^-1)*sum;
        // I_(2,4) when a < b
        when 4:
            cont := 0;
            if ((a-b) mod 2) ne 0 then cont := f_1(a2,b,l,p);
            elif ((a-b) mod 2) eq 0 then cont := KroneckerSymbol(a1,p)*f_2(a2,b,l,p);
            end if;
            sum := 0;
            for k := 1 to a do
                sum +:= v_k(b+1,e,l,p)*v_k(k,e,l,p)*delta_k(k,l,1)*POP(a+b+1+d_k(b+1,l)+d_k(k,l))*X^(b+1+k);
            end for;
            return cont*sum;
        // I_(2,5) when a < b
        when 5:
            sum := 0;
            for k := (a+2) to b do
                sum +:= v_k(k,e,l,p)*g_k(k,l,a1,p,a)*POP(a+k+d_k(a+1,l)+d_k(k,l))*X^(a+1+k);
            end for;
            return v_k(a+1,e,l,p)*f_1(a1,a,l,p)*sum;
        // I_(2,6) when a < b
        when 6:
            cont := 0;
            if ((a-b) mod 2) ne 0 then cont := f_1(a2,b,l,p);
            elif ((a-b) mod 2) eq 0 then cont := KroneckerSymbol(a1,p)*f_2(a2,b,l,p);
            end if;
            return v_k(a+1,e,l,p)*v_k(b+1,e,l,p)*f_1(a1,a,l,p)*POP(a+b+1+d_k(a+1,l)+d_k(b+1,l))*cont*X^(a+b+2);
        // I_(2,7) when a < b
        when 7:
            return delta_k(a+1,l,-1)*POP(2*a+2*d_k(a+1,l))*X^(2*a+2);
        // I_(2,8) when a < b
        when 8:
            sum := 0;
            for k := 1 to a do
                sum +:= (delta_k(k,l,1)+p^-1*delta_k(k,l,-1))*POP(2*k+2*d_k(k,l))*X^(2*k);
            end for;
            return sum;
        else error "Error. For I_(2,i) when a < b, i must be in {1..8}.";
    end case;
end function;

////////////////////////////////////////////////////////////////////////////////////////////////

// R_0 or R_1 or R_2
R := function(n,e,l,a1,a2,p,a,b)
    r_val := 0;
    if n eq 0 then r_val := 1;
    elif n eq 1 then
        if a eq b then
            for i := 1 to 4 do
                r_val +:= (W ! I_1_i_aeqb(i,e,l,a1,a2,p,a));
            end for;
        elif a lt b then
            for i := 1 to 4 do
                r_val +:= (W ! I_1_i_altb(i,e,l,a1,a2,p,a,b));
            end for;
        else error "Error. a must be less than or equal to b.";
        end if;
    elif n eq 2 then
        sum := 0;
        val_I26 := 0;
        if a eq b then
            for i := 1 to 8 do
                sum +:= (W ! I_2_i_aeqb(i,e,l,a1,a2,p,a));
            end for;
            val_I26 := I_2_i_aeqb(6,e,l,a1,a2,p,a);
        elif a lt b then
            for i := 1 to 8 do
                sum +:= (W ! I_2_i_altb(i,e,l,a1,a2,p,a,b));
            end for;
            val_I26 := I_2_i_altb(6,e,l,a1,a2,p,a,b);
        else error "Error. a must be less than or equal to b.";
        end if;
        r_val := (1-p^-1)*sum+p^-1*val_I26;
    else error "Error. For R_i(X), i must be in {0,1,2).";
    end if;
    return (W ! (K !  r_val));
end function;

////////////////////////////////////////////////////////////////////////////////////////////////

// beta(X)^l
beta_l := function(e,l,a1,a2,p,a,b,l_val)
    b_val := (1-p^-2)*p^(2*l_val)*X^l_val+(1-p^-1)*p^(3*l_val)*X^(2*l_val)*(1+R(1,e,l,a1,a2,p,a,b));
    return (W ! (K ! b_val));
end function;

/////////////////////////////////////////////////////////////////////////////////////////////////

// alpha(X, T^l,S^l)
alpha_X := function(e,l,a1,a2,p,a,b,l_val)
    val :=  (W ! (K ! (1+p^(2*l_val)*X^(l_val)*R(1,e,l,a1,a2,p,a,b)+p^(3*l_val)*X^(2*l_val)*R(2,e,l,a1,a2,p,a,b)+l_val
        *beta_l(e,l,a1,a2,p,a,b,l_val))));
   // return Evaluate(val,4,p,5,1/p);
    val1 := Evaluate(val, 4, p);
    return Evaluate(val1, 5, 1/p);
end function;

///////////////////////////////////////////////////////////////////////////////////////////////////