// %%%%%%%%% PROJECT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

/* Update Needed
1. How to get an output flashing buffer?
2. How to terminate the program during run after an error is thrown?
3. How to get a stringTokenizer to read data from a text file?
4. A large set of random input data to test the consistence of the program.
*/
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// ^^^^^^^
//        ^^^^^^^^
//        ^^^^^^^^^^^
//              ^^^^^^^^^^^^^^^^^
//                          ^^^^^^^^^^^^^^^^^^^^
//

//@@@@@@@@@@@@@@@@@@@@@@@@@ Project A @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
/* Input data:
    a   : a >= -1
    e   : e_i in Z*_2, i = 1,2,...,H
    l   : l_i >= 1, i = 1,2,...,H
    m   : m_j >= 0 are integers, j = 1,2,...,M
    n   : n_k >= 0 are integers, k = 1,2,...,N
*/
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

/* Functions:
    F := function(r)
        F_1 := function(a,c)
        F_2 := function(a,c)
        SortL := function(l,m,n)
        G_1 := function(c,e,l,m,n)
        G_2 := function(c,e,l,m,n)
        I_1 := function(alpha,c)
        A := function(beta,c,e,l,m,n)
        A_r := function(c,e,l,m,n,r)
        S := function(c,a,e,l,m,n)
        W := function(a,e,l,m,n)

*/
// Compute f(r)
F := function(r)
        if r ge 3 then return 1/4;
        elif r ge 1 and r le 2 then return 1/2;
        else return "function f(r) fails for invalid parameter value.";
        end if;
end function;



// Compute f_1(a,c)
F_1 := function(a,c)
        if c ge 1 and c le (a+1) then return 3;
        elif c eq (a+2) then return 1;
        else return "function f_1(a,c) fails for invalid parameter values.";
        end if;
end function;



// Compute f_2(a,c)
F_2 := function(a,c)
        if c gt 0 and c le (a+1) then return 1;
        elif c eq (a+2) then return -1;
        else return "function f_2(a,c) fails for invalid parameter values.";
        end if;
end function;


// Reorder l_i,m_j,m_j,n_k,n_k into an increasing sequence
// 0 <= l'_1 <= l'_2 <= ... <= l'_n
SortL := function(l,m,n)
        S := [];
        for idx := 1 to #(l) do
            S[idx] := l[idx];
        end for;
        if IsNull(m) eq false then
            for idx1 := 1 to #(m) do
                    S[#(l)+idx1] := m[idx1];
            end for;

            for idx2 := 1 to #(m) do
                    S[#(l)+#(m)+idx2] := m[idx2];
            end for;
        end if;

        if IsNull(n) eq false then
            for idx3 := 1 to #(n) do
                    S[#(l)+2*#(m)+idx3] := n[idx3];
            end for;

            for idx4 := 1 to #(n) do
                    S[#(l)+2*#(m)+#(n)+idx4] := n[idx4];
            end for;
        end if;
        S := Sort(S);

        return S;
end function;


// Compute g_1(c),
G_1 := function(c,e,l,m,n)
        l_p := SortL(l,m,n);
        n := #(l_p);
        g_1 := 0;

        for j := 1 to n do
            g_1 +:= Min(0,l_p[j]-c);
        end for;

        return g_1;
end function;


// Compute g_2(c)
G_2 := function(c,e,l,m,n)
        g1 := G_1(c,e,l,m,n);
        return (-1)^g1;
end function;



// Compute I_1(alpha*2^c)
I_1 := function(alpha,c)
        rval:= 1;
        R<P,P_,X> := PolynomialRing(Rationals(),3);
        K := quo<R|P^4+1,P_^4+1, P_*P-1>;

        if c gt 0 then rval := 1;
        elif c eq 0 then rval := 0;
        elif c lt 0 and IsEven(c) then
            if alpha ge 0 then
                    rval := 2^(Integers() ! (c/2))*P^alpha;
            else rval :=2^(Integers() ! (c/2))*P_^(-1*alpha);
            end if;
        elif c lt 0 and IsOdd(c) then
            if alpha ge 0 then
                    phi := P^(2*alpha);
            else phi := P_^(-2*alpha);
            end if;

            rval := 2^(Integers() ! ((c-1)/2))*(1+phi);
        end if;

        return (R ! (K ! rval));
end function;


// Compute A(beta*2^c,L)
A := function(beta,c,e,l,m,n)
        R<P,P_,X> := PolynomialRing(Rationals(),3);
        K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        prod1 := 1;
        for id := 1 to #(l) do
            prod1 *:= (R ! I_1(e[id]*beta,c+l[id]));
        end for;

        prod2 := 1;
        if IsNull(m) eq false then
            for j := 1 to #(m) do;
                    prod2 *:= 2^(Min(0,c+m[j]));
            end for;
        end if;

        prod3 := 1;
        if IsNull(n) eq false then
            for k := 1 to #(n) do
                    prod3 *:= (-2)^(Min(0,c+n[k]));
            end for;
        end if;

        return (R ! (K ! (prod1*prod2*prod3)));
end function;



// A_sum(c,L) = SUM_(beta in (Z/8)*) A(beta*2^c,L)
A_sum := function(c,e,l,m,n)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
        K := quo<R|P^4+1,P_^4+1, P*P_-1>;
        a_1 := (R ! A(1,c,e,l,m,n));   // A(2^-c,L)
        a_2 := (R ! A(-1,c,e,l,m,n));  // A(-2^-c,L)
        a_3 := (R ! A(3,c,e,l,m,n));       // A(3*2^-c,L)
        a_4 :=  (R ! A(-3,c,e,l,m,n));  // A(-3*2^-c,L)

        return (R ! (K ! (a_1 + a_2 + a_3 + a_4)));
end function;


// S_{<0}(-c, T, L)
S_lt := function(a,c,e,l,m,n)
    return 2^(3*c-5+G_1(c,e,l,m,n))*(3+G_2(c,e,l,m,n))*F_2(a,c);
end function;


// S_r(-c,T,L)
S_r := function(a,c,r,e,l,m,n)
    R<P,P_, X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        sum1 := (R ! (F_1(a,c)*2^(3*c-6-r)*A_sum(r-c,e,l,m,n)));
        sum2 := 1;

        if r eq 0  then
            sum2 := F_2(a,c)*((R ! A(1,-c,e,l,m,n)) + (R ! A(-1,-c,e,l,m,n)));
        elif r eq 1 then
        sum2 := (R ! A(1,-c,e,l,m,n))+( R ! A(-3,-c,e,l,m,n));
        elif r eq 2 then
            sum2 := (R ! A(1,-c,e,l,m,n))+ (R ! A(-1,-c,e,l,m,n));
        elif r ge 3 then
                sum2 := (1/2)*(R ! A_sum(-c,e,l,m,n));
        end if;

        return (R ! (K ! (sum1*sum2)));
end function;


// S_{>=c} (-c,T,L)
S_gec := function(a,c,e,l,m,n)
    R<P,P_, X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        sum0 := 1;
        if c ge 3 then
                sum0 := 2^(2*c-4)*F_1(a,c)*(R ! A_sum(-c,e,l,m,n));
        elif c ge 1 and c lt 3 then
                sum0 := 0;

        for r := c to 2 do
            sum0 +:= (R ! S_r(a,c,r,e,l,m,n));
                end for;
        sum0 := sum0 + 2^(3*c-7)*F_1(a,c)*(R ! A_sum(-c,e,l,m,n));
        end if;

        return  (R ! (K ! sum0));
end function;

// Compute W(X,T,L)
W := function(a,e,l,m,n)
    R<P,P_, X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        s1 := 0;
    for c := 1 to (a+2) do
        s1 +:= (R ! S_lt(a,c,e,l,m,n))*X^(2*c);
    end for;
    s1 := (R ! s1);

        s2 := 0;
        for c := 1 to (a+2) do
            for r := 0 to (c-1) do
            s2 +:=  ( R ! S_r(a,c,r,e,l,m,n))*X^(2*c-r);
        end for;
    end for;
    s2 := (R ! s2);

        s3 := 0;
    for c := 1 to (a+2) do
        s3 +:=(R ! S_gec(a,c,e,l,m,n))*X^c;
    end for;
    s3 := (R ! s3);

    return (R ! (K ! (1+s1+s2+s3)));

 end function;



// S_{<0}(-c, T, L)
P2_S_lt := function(a,c,e,l,m,n)
    return 2^(3*c-5+G_1(c,e,l,m,n))*(3+G_2(c,e,l,m,n))*F_2(a,c);
end function;


// S_r(-c,T,L)
P2_S_r := function(a,c,r,e,l,m,n)
    R<P,P_, X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        sum1 := (R ! (3*2^(3*c-6-r)*A_sum(r-c,e,l,m,n)));
        sum2 := 1;

        if r eq 0  then
            sum2 := ((R ! A(1,-c,e,l,m,n)) + (R ! A(-1,-c,e,l,m,n)));
        elif r eq 1 then
        sum2 := F_2(a,c)*((R ! A(1,-c,e,l,m,n))+( R ! A(-3,-c,e,l,m,n)));
        elif r eq 2 then
            sum2 := F_2(a,c)*( (R ! A(1,-c,e,l,m,n))+ (R ! A(-1,-c,e,l,m,n)));
        elif r ge 3 then
                sum2 := (1/2)* F_2(a,c)*(R ! A_sum(-c,e,l,m,n));
        end if;

        return (R ! (K ! (sum1*sum2)));
end function;


// S_{>=c} (-c,T,L)
P2_S_gec := function(a,c,e,l,m,n)
    R<P,P_, X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        sum0 := 1;
        if c ge 3 then
                sum0 := 3*2^(2*c-4)*F_2(a,c)*(R ! A_sum(-c,e,l,m,n));
        elif c ge 1 and c lt 3 then
                sum0 := 0;

        for r := c to 2 do
            sum0 +:= (R ! P2_S_r(a,c,r,e,l,m,n));
                end for;
        sum0 := sum0 + 3*2^(3*c-7)*F_2(a,c)*(R ! A_sum(-c,e,l,m,n));
        end if;

        return  (R ! (K ! sum0));
end function;

// Compute W(X,T,L)
P2_W := function(a,e,l,m,n)
    R<P,P_, X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

        s1 := 0;
    for c := 1 to (a+2) do
        s1 +:= (R ! P2_S_lt(a,c,e,l,m,n))*X^(2*c);
    end for;
    s1 := (R ! s1);

        s2 := 0;
        for c := 1 to (a+2) do
            for r := 0 to (c-1) do
            s2 +:=  ( R ! P2_S_r(a,c,r,e,l,m,n))*X^(2*c-r);
        end for;
    end for;
    s2 := (R ! s2);

        s3 := 0;
    for c := 1 to (a+2) do
        s3 +:=(R ! P2_S_gec(a,c,e,l,m,n))*X^c;
    end for;
    s3 := (R ! s3);

    return (R ! (K ! (1+s1+s2+s3)));

 end function;


// &&&&&&&&&&&  PROJECT B   &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

// Imput data
// a,c,r,t: integers
// a1,a2:   ai in Z*/8
// b1,b2:   bi in Z*/8


// Reused functions
// I_1 := function(alpha,c)

/* functions

Case1_B1 := function(a,c,r,a1,a2,b1,b2)
Case1_B2 := function(a,c,r,a1,a2,b1,b2)
Case1_B3 := function(a,c,r,a1,a2,b1,b2)
Case2_B1 := function(a,c,t,r,a1,a2,b1,b2)
case2_B2 := function(a,c,t,r,a1,a2,b1,b2)
case2_B3 := function(a,c,t,r,a1,a2,b1,b2)
B := function(a,c,t,r,a1,a2,b1,b2)
*/
// &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

// compute B_1(gamma,T) when t = 0
Case1_B1 := function(a,c,r,a1,a2,b1,b2)
    case1_b1 := 1;
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

    if c le 0 then
        printf "Error thrown as c <= 0. Ignore all the runtime errors.";
        printf "\n";
    elif c gt 0 and c le (a+3) then
        d1 := (-2^(a-c+3)*(a1*b1+2^r*a2*b2) mod 8); // *********d1***********
        prd := 1;
        if d1 lt 0 then
            prd := P_^(-1*d1);
        elif d1 eq 0 then
            prd := 1;
        else prd := P^d1;
        end if;

        case1_b1 := 2^(-3)*(R ! I_1((-1)*a2*b1,a+1-c))*(R ! I_1((-1)*a1*b2,a+3+r-c))*prd;
    elif c gt (a+3) then
        case1_b1 := 0;
    end if;

    return (R ! (K ! case1_b1));
end function;

// compute B_2(gamma,T) when t = 0
Case1_B2 := function(a,c,r,a1,a2,b1,b2)
    case1_b2 := 1;
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

    if c le 0 then
        printf "Error thrown as c <= 0. Ignore all the runtime errors.";
        printf "\n";
    elif c gt 0 and c le (a+3) then
        d2 := (-2^(a-c+3)*(a2*b1+2^r*a1*b2) mod 8); // ********d2***************
        prd := 1;
        if d2 lt 0 then
            prd := P_^(-1*d2);
        elif d2 eq 0 then
            prd := 1;
        else prd := P^d2;
        end if;

        case1_b2 := 2^(-3)*(R ! I_1((-1)*a1*b1,a+1-c))*(R ! I_1((-1)*a2*b2,a+3+r-c))*prd;
    elif c gt (a+3) then
        case1_b2 := 0;
    end if;

    return (R ! (K ! case1_b2));
end function;

// compute B_3(gamma,T) when t = 0
Case1_B3 := function(a,c,r,a1,a2,b1,b2)
    case1_b3 := 1;
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    modu := (a2-a1) mod 4;

    if modu eq 0 then
        if c le 0 then
            printf "Error thrown as c <= 0. Ignore all the runtime errors.";
            printf "\n";
        elif c gt 0 and c le (a+2) then
            d3 := (-2^(a-c+3) mod 8)*a1*(b1+(2^(r+1) mod 8)*b2); // *******d3**************
            prd := 1;
            if d3 eq 0 then
                prd := 1;
            elif d3 gt 0 then
                prd := P^d3;
            else prd := P_^(-1*d3);
            end if;

            case1_b3 := 2^(-3)*prd;
        elif c gt (a+2) then
            case1_b3 := 0;
        end if;
        return case1_b3;
    elif modu ne 0 then
        if c le 0 then
            printf "Error thrown as c <= 0. Ignore all the runtime errors.";
            printf "\n";
        elif c gt (a+1) then
            case1_b3 := 0;
        elif c eq (a+1) then
            case1_b3 :=-1* 2^(-3);
        elif c gt 0 and c le a then
            case1_b3 := 2^(-3);
        end if;
    end if;

    return (R ! (K ! case1_b3));
end function;

// compute B_1(gamma,T) when t > 0
Case2_B1 := function(a,c,t,r,a1,a2,b1,b2)
    case2_b1 := 1;
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;

    if c le 0 then
        printf "Error thrown as c <= 0. Ignore all the runtime errors.";
        printf "\n";
    elif c gt (a+3) then
        case2_b1 := 0;
    elif c gt 0 and c le (a+3) then
        d_1 := -1*(2^(a+3-c) mod 8)*(a1*b1+a2*b2*(2^(t+r) mod 8)); // ******d_1*************
        prd := 1;
        if d_1 eq 0 then
            prd := 1;
        elif d_1 gt 0 then
            prd := P^d_1;
        elif d_1 lt 0 then
            prd := P_^(-1*d_1);
        end if;
        case2_b1 := (1/8)*(R ! I_1(-1*a2*b1,a+t+1-c))*(R ! I_1(-1,a+3-c+r))*prd;
    end if;

    return (R ! (K ! case2_b1));
end function;

// compute B_2(gamma,T) when t > 0
Case2_B2 := function(a,c,t,r,a1,a2,b1,b2)
    case2_b2 := 1;
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    min := Min(a+t+3,a+3+r);

    if c le 0 then
        printf "Error thrown as c <= 0. Ignore all the runtime errors.";
        printf "\n";
    elif c gt 0 and c le min then
           d_2 :=-1*((2^(a+t+3-c) mod 8)*a2*b1+(2^(a+r+3-c) mod 8)*a1*b2); // **d_2*****
        prd := 1;
        if d_2 eq 0 then
            prd := 1;
        elif d_2 gt 0 then
            prd := P^d_2;
        else prd := P_^(-1*d_2);
        end if;
        case2_b2 := (1/8)*(R ! I_1(-1*a1*b1,a+1-c))*(R ! I_1(-1,a+t+3-c+r))*prd;
    elif c gt min then
        case2_b2 := 0;
    end if;

    return (R ! (K ! case2_b2));
end function;

// compute B_3(gamma,T) when t > 0
Case2_B3 := function(a,c,t,r,a1,a2,b1,b2)
    case3_b3 := 1;
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    min := Min(a+t+3,a+r+3);

    alpha := a1+( 2^t mod 8)*a2;


    if c le 0 then
        printf "Error thrown as c <= 0. Ignore all the runtime errors.";
        printf "\n";
    elif c gt 0 and c le min then
        d_3 := (-1*2^(a+t+3-c) mod 8)*a1*a2*b1*alpha-(2^(a+r+3-c) mod 8)*a1*b2-(2^(a+t+r+3-c) mod 8)*a2*b2; // **d_3**
        prd := 1;
        if d_3 eq 0 then
            prd := 1;
        elif d_3 gt 0 then
            prd := P^d_3;
        else prd := P_^(-1*d_3);
        end if;
        case2_b3 := (1/8)*(R ! I_1(-1*(a1*b1+2^t*a2*b1),a+1-c))*prd;
    elif c gt min then
        case2_b3 := 0;
    end if;

    return (R ! (K ! case2_b3));

end function;

// MAIN FUNCTION
// compute B(T,gamma) = B_1(gamma,T)+B_2(gamma,T)+B_3(gamma,T)
B := function(a,c,t,r,a1,a2,b1,b2)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1,P*P_-1>;
    B := 1;

    if t eq 0 then
        B := (R ! Case1_B1(a,c,r,a1,a2,b1,b2))+(R ! Case1_B2(a,c,r,a1,a2,b1,b2))+(R ! Case1_B3(a,c,r,a1,a2,b1,b2));
    elif t gt 0 then
        B := (R ! Case2_B1(a,c,t,r,a1,a2,b1,b2))+(R ! Case2_B2(a,c,t,r,a1,a2,b1,b2))+(R ! Case2_B3(a,c,t,r,a1,a2,b1,b2));
    else
        printf "Error thrown as t < 0. Ignore all the runtime errors.";
        printf "\n";
    end if;

    return (R ! (K ! B));
end function;

// @@@@@@@@@@@@@@@@@@@@@@@@@@@ Project c @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
/*
// Main Function:
    P3_W := function(a,t,e,l,m,n,a1,a2)
// Functions:
    P3_S_lt0 := function(a,c,t,e,l,m,n)
    P3_S_r := function(a,c,t,r,e,l,m,n,a1,a2)
    Phi_4 := function(alpha,beta)
    Beta_t := function(a,c,t,a1,a2,beta,e,l,m,n)
    Help_sum := function(a,c,t,a1,a2,e,l,m,n)
    P3_S_gec := function(a,c,t,a1,a2,e,l,m,n)
*/
// @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

// compute S_{<0}(-c,T,L)
P3_S_lt0 := function(a,c,t,e,l,m,n)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    s := 0;

    if t gt 0 then
        s +:= 2^(3*c-5+G_1(c,e,l,m,n))*(F_1(a,c)+G_2(c,e,l,m,n)*F_2(a,c));
    elif t eq 0 then
        s +:= 2^(3*c-5+G_1(c,e,l,m,n))*(F_1(a,c) +G_2(c,e,l,m,n)*F_2(a,c))*F_2(a,c);

    else
        printf "Error, t has to be greater than or equal to 0.";
        printf "\n";
    end if;

    return (R ! (K ! s));
end function;

// compute S_r(-c,T,L)
P3_S_r := function(a,c,t,r,e,l,m,n,a1,a2)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    p3_s_r := 0;

    if r eq 0 then
        //B := function(a,c,t,r,a1,a2,b1,b2)
        // compute SUM_{gamma in A_{0}^(-c)} [B(T,gamma)A(gamma,L)]
        sum_1 := (R ! B(a,c,t,0,a1,a2,1,-1))*(R ! A(1,-c,e,l,m,n))*(R ! A(-1,-c,e,l,m,n));
        sum_2 := (R ! B(a,c,t,0,a1,a2,1,3))*(R ! A(1,-c,e,l,m,n))*(R ! A(3,-c,e,l,m,n));

        sum_p1 := (R ! B(a,c,t,0,a1,a2,1,1))*(R ! A(1,-c,e,l,m,n))*(R ! A(1,-c,e,l,m,n));
        sum_p2 := (R ! B(a,c,t,0,a1,a2,-1,-1))*(R ! A(-1,-c,e,l,m,n))*(R ! A(-1,-c,e,l,m,n));
        sum_p3 := (R ! B(a,c,t,0,a1,a2,1,-3))*(R ! A(1,-c,e,l,m,n))*(R ! A(-3,-c,e,l,m,n));
        sum_p4 := (R ! B(a,c,t,0,a1,a2,-1,3))*(R ! A(-1,-c,e,l,m,n))*(R ! A(3,-c,e,l,m,n));
        p3_s_r := 2^(3*c-2)*(sum_1+sum_2) + 2^(3*c-3)*(sum_p1+sum_p2+sum_p3+sum_p4);
    elif r eq 1 then
        sum_s1 := 0;
        for beta1 in {1,-3} do
            for beta2 in {1,-1,3,-3} do
                sum_s1 +:= (R ! B(a,c,t,r,a1,a2,beta1,beta2))*(R ! A(beta1,-c,e,l,m,n))*(R ! A(beta2,1-c,e,l,m,n));
            end for;
        end for;
        p3_s_r := 2^(3*c-r-3)*sum_s1;
    elif r eq 2 then
        sum_s2 := 0;
        for beta1 in {1,-1} do
            for beta2 in {1,-1,3,-3} do
                sum_s2 +:= (R ! B(a,c,t,r,a1,a2,beta1,beta2))*(R ! A(beta1,-c,e,l,m,n))*(R ! A(beta2,2-c,e,l,m,n));
            end for;
        end for;
        p3_s_r := 2^(3*c-r-3)*sum_s2;
    elif r ge 3 then
        sum_s3 := 0;
        for beta1 in {1,-1,3,-3} do
            for beta2 in {1,-1,3,-3} do
                sum_s3 +:= (R ! B(a,c,t,r,a1,a2,beta1,beta2))*(R ! A(beta1,-c,e,l,m,n))*(R ! A(beta2,r-c,e,l,m,n));
            end for;
        end for;
        p3_s_r := 2^(3*c-r-4)*sum_s3;
    end if;

    return (R ! (K ! p3_s_r));
end function;

// compute phi(d/4)
Phi_4 := function(alpha,beta)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    d := -1*alpha*beta;
    p_4 := 0;

    if d gt 0 then
        p_4 := P^(2*d);
    elif d eq 0 then
        p_4 := 1;
    else
        p_4 := P_^(-2*d);
    end if;

    return p_4;
end function;


// compute B(T,beta*2^-c) for the cases t= 0 and t >0
Beta_t := function(a,c,t,a1,a2,beta,e,l,m,n)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    beta_t := 0;

    if t eq 0 then
        alpha := a1+a2;
        if c gt (a+3) then
            beta_t := 0;
        elif c eq (a+3) then
            d := -1*alpha*beta;
            if d gt 0 then
                beta_t := P^d;
            elif d eq 0 then
                beta_t := 1;
            else
                beta_t := P_^(-1*d);
            end if;
        elif c eq (a+2) then
            beta_t := (R ! Phi_4(alpha,beta))+(R ! Phi_4(a1,beta))+(R ! Phi_4(a2,beta));
        elif c eq (a+1) then
            beta_t := -1;
        elif c gt 0 and c le a then
            beta_t := 3;
        else
            printf "Error, c must be greater than or equal to zero.";
            printf "\n";
        end if;
    elif t gt 0 then
        if c gt (a+t+3) then
            beta_t := 0;
        elif c gt (a+3) and c le (a+t+3) then

            // compute psi(-a2*beta*2^(a+t-c))
            d1 := (-2^(a+t+3-c) mod 8)*a2*beta;
            ph1 := 1;
            if d1 eq 0 then
                ph1 := 1;
            elif d1 gt 0 then
                ph1 := P^d1;
            elif d1 lt 0 then
                ph1 := P_^(-1*d1);
            end if;
            // sum1 := R ! I_1(-1*a1*beta,a+1-c))*ph1;
             //compute phi(-a1*a2*alpha*beta*2^(a+t-c))
             d2 :=(-2^(a+t+3-c) mod 8)*a1*a2*(a1+2^t*a2)*beta;
             ph2 :=1;
              if d2 eq 0 then
                ph2 := 1;
              elif d2 gt 0 then
                ph2 := P^d2;
              elif d2 lt 0 then
                ph2 := P_^(-1*d2);
              end if;
             //sum2 := (R ! I_1(-1*(a1+2^t*a2)*beta,a+1-c))*ph2;
             //beta_t := sum1 +sum2;

            beta_t :=  (R ! I_1(-1*a1*beta,a+1-c))*ph1
             +(R ! I_1(-1*(a1+2^t*a2)*beta,a+1-c))*ph2;

        elif c gt 0 and c le (a+3) then

         // compute phi(-a1*beta*2^(a-c))
            d0 := -2^(a+3-c)*a1*beta;
            ph0 := 1;
            if d0 eq 0 then
                ph0 := 1;
            elif d0 gt 0 then
                ph0 := P^d0;
            else
                ph0 := P_^(-1*d0);
            end if;
       // compute phi(-a2*beta*2^(a+t-c))
            d1 := -2^(a+t+3-c)*a2*beta;
            ph1 := 1;
            if d1 eq 0 then
                ph1 := 1;
            elif d1 gt 0 then
                ph1 := P^d1;
            else
                ph1 := P_^(-1*d1);
            end if;
     //compute phi(-a1*a2*alpha*beta*2^(a+t-c))
             d2 :=-2^(a+t+3-c)*a1*a2*(a1+2^t*a2)*beta;
             ph2 :=1;
              if d2 eq 0 then
                ph2 := 1;
              elif d2 gt 0 then
                ph2 := P^d2;
              else
                ph2 := P_^(-1*d2);
              end if;

          beta_t := (R ! I_1(-1*a2*beta,a+t+1-c))*ph0+(R ! I_1(-1*a1*beta,a+1-c))*ph1
             +(R ! I_1(-1*(a1+2^t*a2)*beta,a+1-c))*ph2;
      end if;
    elif t lt 0 then
        printf "Error, t must be greater than or equal to zero.";
        printf "\n";
    end if;

    return (R ! (K ! beta_t));
end function;

// compute SUM_{beta in (Z/8)*}[B(T,beta*2^-c)A(beta*2^-c,L)]
Help_sum := function(a,c,t,a1,a2,e,l,m,n)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    hs := 0;

    for beta in {1,-1,3,-3} do
        hs +:= (R ! Beta_t(a,c,t,a1,a2,beta,e,l,m,n) )*(R ! A(beta,-c,e,l,m,n));
    end for;

    return (R ! (K ! hs));
end function;

// compute S_{>=c}(-c,T,L)
P3_S_gec := function(a,c,t,a1,a2,e,l,m,n)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1, P*P_-1>;
    s_g := 0;

    if c ge (t+3) then
        s_g := 2^(2*c-4)*(R ! Help_sum(a,c,t,a1,a2,e,l,m,n));
    elif c gt 0 and c lt (t+3) then
        sum_c := 0;
        for r := c to (t+2) do
            sum_c +:= (R ! P3_S_r(a,c,t,r,e,l,m,n,a1,a2));
        end for;
        s_g := sum_c + 2^(3*c-t-7)*(R ! Help_sum(a,c,t,a1,a2,e,l,m,n));
    elif c le 0 then
        printf "Error, c must be greater than zero.";
        printf "\n";
    end if;

    return (R ! (K ! s_g));
end function;

// main function: compute W(X,T,L)
P3_W := function(a,t,e,l,m,n,a1,a2)
    R<P,P_,X> := PolynomialRing(Rationals(),3);
    K := quo<R|P^4+1,P_^4+1,P*P_-1>;
    w := 0;
    doable := 1;

/*
    if a lt 0 or t lt 0 then
        printf "Error, a, t must be greater than or equal to zero.";
        printf "\n";
        doable := 0;
    end if;



    dd := a1*a2;
    if (dd mod 2) eq 0 then
        printf "Error, a_1 and a_2 have to be odd.";
        printf "\n";
        doable := 0;
    end if;

    for i := 1 to #(e) do
        if (e[i] mod 2) eq 0 then
            printf "Error, e_i has to be odd.";
            printf "\n";
            doable := 0;
            break;
        end if;
    end for;

    if #(e) ne #(l) then
        printf "Error, # of e_i's must equal # of l_i's.";
        printf "\n";
        doable := 0;
    end if;
*/
    if doable eq 1 then
        sum1 := 0;
        for c := 1 to (a+2) do
            sum1 +:= (R ! P3_S_lt0(a,c,t,e,l,m,n))*X^(2*c);
        end for;
        sum2 := 0;
        for c := 1 to (a+t+3) do
            for r := 0 to (c-1) do
                sum2 +:= (R ! P3_S_r(a,c,t,r,e,l,m,n,a1,a2))*X^(2*c-r);
            end for;
        end for;
        sum3 := 0;
        for c := 1 to (a+t+3) do
            sum3 +:= (R ! P3_S_gec(a,c,t,a1,a2,e,l,m,n))*X^c;
        end for;
        w := 1 + sum1 + sum2 + sum3;
    end if;

    return (R ! (K ! w));
end function;


//***************** End of Project c *****************************************************************