// PostProcessingTM_GivenConductor.m

/*
INPUT: 
    C:= [c_1, c_2, c_3, c_4], the coefficients of the Thue-Mahler form, F, where
        F:= c_1*x^3 + c_2*x^2*y + c_3*x*y^2 + c_4*y^3, having discriminant DF:= +/- 2^{alpha_0} * 3^{beta_0} * N_1
    
    N:= the conductor of the elliptic curves we seek
    
    p:= the set of rational primes p_1, ..., p_s such that 
        F(u,v) = 2^{alpha_1} * 3^{beta_1} * p_3^{k_3} * p_4^{k_4} *...* p_s^{k_s}
    alpha_1beta_1:= [ [ [alpha_{1,1},A_1], .., [alpha_{1,r}, A_r] ], [ [beta_{1,1},B_1], .., [beta_{1,r}, B_r] ] ] where
        alpha_{1,i}:= the possible value of alpha_1 (ie. the exponent on 2) based on alpha_0 and alpha
        A_1:= 0 or 1, such that 
            if A_1 == 0 then
                alpha_1 is exactly equal to alpha_{1,i}
            if A_1 == 1 then
                alpha_1 is greater or equal to alpha_{1,i}
        beta_{1,i}:= the possible value of beta_1 (ie. the exponent on 3) based on beta_0 and beta
        B_1:= 0 or 1, such that 
            if B_1 == 0 then
                beta_1 is exactly equal to beta_{1,i}
            if B_1 == 1 then
                beta_1 is greater or equal to beta_{1,i}
    
     Soln:= [sol_1, ..., sol_t], where
        sol_1:= [u,v, exp_3, ..., exp_s] such that
            F(u,v) = 2^{alpha_1} * 3^{beta_1} * p_3^{exp_3} * p_4^{exp_4} *...* p_s^{exp_s}
    
    RSolnFile:= text file storing the solutions of the form F which may lead to elliptic curves of conductor N
    
OUTPUT:
     Soln2:= [sol_1, ..., sol_t], where
        sol_1:= [u,v, exp_1, ..., exp_s] such that
            F(u,v) = 2^{alpha_1} * 3^{beta_1} * p_3^{k_3} * p_4^{k_4} *...* p_s^{k_s}
        such that GCD(u,v) = 1 and such that alpha_1, beta_1 are as in Theorem 1 of the reference
                
COMMENTS:
    Given a Thue-Mahler form F, F:= c_1*x^3 + c_2*x^2*y + c_3*x*y^2 + c_4*y^3, having discriminant DF:= +/- 2^{alpha_0} * 3^{beta_0} * N_1, and conductor 
        N = 2^{alpha} * 3^{beta} * N_0, 
    one can compute solutions to 
        F(u,v) = 2^{alpha_1} * 3^{beta_1} * p_3^{k_3} * p_4^{k_4} *...* p_s^{k_s}, where p_3, ..., p_s | N_0
    leading to elliptic curves having good reduction outside p, as per Theorem 1 of the reference. 
    
    This algorithm determines a potential upper bound on alpha_1, beta_1 as per Theorem 1 of the reference based on alpha_1, beta_1, alpha, beta.
    
REFERENCE:
    M. Bennett, A. Gherga, A. Rechnitzer. Computing Elliptic Curves over Q. Submitted, University of British Columbia, 2017.

EXAMPLE:
    > C:=[ 2, 12, 0, 13 ]; 
    > N:= 10^6 -1;
    > A:=1;
    > RSolnFile:= "RelevantSoln.txt";
    > p, alpha_1beta_1, EnterTM, Bounda, Boundb:= PreProcessingTM_GivenConductor(C,N,RSolnFile);
    > p;
    [ 2, 3, 7, 11, 13, 37 ];
    > alpha_1beta_1;
    [
        [
            [ 0, 0 ],
            [ 3, 0 ]
        ],
        [
            [ 0, 0 ],
            [ 1, 0 ]
        ]
    ]
    > Soln:= [
    >     [ 18, 23, 0, 0, 2, 1, 1, 1 ],
    >     [ -18, -23, 0, 0, 2, 1, 1, 1 ],
    >     [ -6, 1, 0, 0, 0, 0, 1, 0 ],
    >     [ -2, -1, 0, 0, 1, 1, 0, 0 ],
    >     [ 6, -1, 0, 0, 0, 0, 1, 0 ],
    >     [ 2, 1, 0, 0, 1, 1, 0, 0 ],
    >     [ 0, 1, 0, 0, 0, 0, 1, 0 ],
    >     [ 0, -1, 0, 0, 0, 0, 1, 0 ],
    >     [ -4, 1, 0, 0, 1, 1, 0, 0 ],
    >     [ 4, -1, 0, 0, 1, 1, 0, 0 ],
    >     [ -1, 9, 0, 0, 1, 0, 0, 2 ],
    >     [ 5, -3, 0, 0, 1, 1, 1, 0 ],
    >     [ 1, -9, 0, 0, 1, 0, 0, 2 ],
    >     [ -5, 3, 0, 0, 1, 1, 1, 0 ],
    >     [ -7, -1, 0, 2, 0, 1, 1, 0 ],
    >     [ 740, 167, 0, 2, 6, 1, 2, 0 ],
    >     [ 31, -5, 0, 3, 0, 1, 0, 0 ],
    >     [ 43, 25, 0, 5, 3, 1, 0, 0 ],
    >     [ -1073, 205, 0, 16, 0, 1, 0, 0 ],
    >     [ 7, 1, 0, 2, 0, 1, 1, 0 ],
    >     [ -740, -167, 0, 2, 6, 1, 2, 0 ],
    >     [ -31, 5, 0, 3, 0, 1, 0, 0 ],
    >     [ -43, -25, 0, 5, 3, 1, 0, 0 ],
    >     [ 1073, -205, 0, 16, 0, 1, 0, 0 ],
    >     [ -26, 1, 0, 3, 1, 1, 1, 0 ],
    >     [ 26, -1, 0, 3, 1, 1, 1, 0 ],
    >     [ 4, 1, 0, 2, 0, 0, 0, 1 ],
    >     [ -16, -19, 0, 2, 0, 3, 1, 0 ],
    >     [ -5, 1, 0, 2, 1, 0, 0, 0 ],
    >     [ -17, -47, 0, 2, 1, 1, 3, 0 ],
    >     [ -43, 17, 0, 2, 1, 2, 0, 1 ],
    >     [ 1, 1, 0, 3, 0, 0, 0, 0 ],
    >     [ 4, -5, 0, 3, 1, 0, 1, 0 ],
    >     [ 1549, -251, 0, 4, 1, 0, 0, 2 ],
    >     [ 8, -1, 0, 5, 0, 0, 0, 0 ],
    >     [ 155, 11, 0, 5, 1, 0, 2, 1 ],
    >     [ 25, 7, 0, 6, 0, 2, 0, 0 ],
    >     [ 136, 145, 0, 6, 1, 1, 0, 2 ],
    >     [ 1076, -175, 0, 9, 0, 0, 1, 1 ],
    >     [ -193, 257, 0, 11, 2, 0, 0, 1 ],
    >     [ -4, -1, 0, 2, 0, 0, 0, 1 ],
    >     [ 16, 19, 0, 2, 0, 3, 1, 0 ],
    >     [ 5, -1, 0, 2, 1, 0, 0, 0 ],
    >     [ 17, 47, 0, 2, 1, 1, 3, 0 ],
    >     [ 43, -17, 0, 2, 1, 2, 0, 1 ],
    >     [ -1, -1, 0, 3, 0, 0, 0, 0 ],
    >     [ -4, 5, 0, 3, 1, 0, 1, 0 ],
    >     [ -1549, 251, 0, 4, 1, 0, 0, 2 ],
    >     [ -8, 1, 0, 5, 0, 0, 0, 0 ],
    >     [ -155, -11, 0, 5, 1, 0, 2, 1 ],
    >     [ -25, -7, 0, 6, 0, 2, 0, 0 ],
    >     [ -136, -145, 0, 6, 1, 1, 0, 2 ],
    >     [ -1076, 175, 0, 9, 0, 0, 1, 1 ],
    >     [ 193, -257, 0, 11, 2, 0, 0, 1 ],
    >     [ 130, -23, 0, 4, 0, 1, 1, 1 ],
    >     [ -130, 23, 0, 4, 0, 1, 1, 1 ],
    >     [ 16778, -2077, 0, 7, 1, 1, 5, 1 ],
    >     [ -16778, 2077, 0, 7, 1, 1, 5, 1 ],
    >     [ -34, 5, 0, 2, 1, 2, 0, 0 ],
    >     [ 34, -5, 0, 2, 1, 2, 0, 0 ],
    >     [ -13, 5, 0, 4, 1, 0, 1, 0 ],
    >     [ -3497, 589, 0, 5, 1, 5, 1, 0 ],
    >     [ 247, 283, 0, 12, 1, 1, 1, 0 ],
    >     [ 13, -5, 0, 4, 1, 0, 1, 0 ],
    >     [ 3497, -589, 0, 5, 1, 5, 1, 0 ],
    >     [ -247, -283, 0, 12, 1, 1, 1, 0 ],
    >     [ -13, 6, 1, 0, 0, 1, 1, 1 ],
    >     [ 13, -6, 1, 0, 0, 1, 1, 1 ],
    >     [ -1, 0, 1, 0, 0, 0, 0, 0 ],
    >     [ 23, -6, 1, 0, 2, 0, 2, 0 ],
    >     [ 1, 0, 1, 0, 0, 0, 0, 0 ],
    >     [ -23, 6, 1, 0, 2, 0, 2, 0 ],
    >     [ 65, -6, 1, 0, 1, 3, 1, 0 ],
    >     [ -65, 6, 1, 0, 1, 3, 1, 0 ],
    >     [ 33, -38, 1, 0, 1, 0, 3, 1 ],
    >     [ -33, 38, 1, 0, 1, 0, 3, 1 ],
    >     [ -5, 4, 1, 4, 0, 1, 0, 0 ],
    >     [ 191, -34, 1, 6, 1, 1, 1, 0 ],
    >     [ 5, -4, 1, 4, 0, 1, 0, 0 ],
    >     [ -191, 34, 1, 6, 1, 1, 1, 0 ],
    >     [ -157, 44, 1, 2, 1, 0, 0, 3 ],
    >     [ 1, 4, 1, 2, 2, 0, 0, 0 ],
    >     [ -7, 2, 1, 3, 0, 1, 0, 0 ],
    >     [ -1, -28, 1, 3, 0, 1, 1, 1 ],
    >     [ 71, -10, 1, 3, 2, 0, 0, 1 ],
    >     [ -133, 20, 1, 4, 0, 0, 3, 0 ],
    >     [ -281, 52, 1, 4, 3, 2, 0, 0 ],
    >     [ 1, 10, 1, 8, 0, 0, 0, 0 ],
    >     [ 157, -44, 1, 2, 1, 0, 0, 3 ],
    >     [ -1, -4, 1, 2, 2, 0, 0, 0 ],
    >     [ 7, -2, 1, 3, 0, 1, 0, 0 ],
    >     [ 1, 28, 1, 3, 0, 1, 1, 1 ],
    >     [ -71, 10, 1, 3, 2, 0, 0, 1 ],
    >     [ 133, -20, 1, 4, 0, 0, 3, 0 ],
    >     [ 281, -52, 1, 4, 3, 2, 0, 0 ],
    >     [ -1, -10, 1, 8, 0, 0, 0, 0 ],
    >     [ -13, 2, 1, 2, 0, 0, 1, 0 ],
    >     [ 65, -34, 1, 3, 4, 0, 1, 0 ],
    >     [ 13, -32, 1, 5, 1, 1, 1, 0 ],
    >     [ 13, -2, 1, 2, 0, 0, 1, 0 ],
    >     [ -65, 34, 1, 3, 4, 0, 1, 0 ],
    >     [ -13, 32, 1, 5, 1, 1, 1, 0 ],
    >     [ 205, -32, 1, 2, 1, 1, 1, 1 ],
    >     [ -205, 32, 1, 2, 1, 1, 1, 1 ],
    >     [ 1, -2, 1, 2, 1, 0, 0, 0 ],
    >     [ 457, -74, 1, 4, 1, 1, 1, 0 ],
    >     [ 17, 8, 1, 5, 1, 0, 1, 0 ],
    >     [ -1, 2, 1, 2, 1, 0, 0, 0 ],
    >     [ -457, 74, 1, 4, 1, 1, 1, 0 ],
    >     [ -17, -8, 1, 5, 1, 0, 1, 0 ],
    >     [ 36, 46, 3, 0, 2, 1, 1, 1 ],
    >     [ -36, -46, 3, 0, 2, 1, 1, 1 ],
    >     [ -12, 2, 3, 0, 0, 0, 1, 0 ],
    >     [ -4, -2, 3, 0, 1, 1, 0, 0 ],
    >     [ 12, -2, 3, 0, 0, 0, 1, 0 ],
    >     [ 4, 2, 3, 0, 1, 1, 0, 0 ],
    >     [ 0, 2, 3, 0, 0, 0, 1, 0 ],
    >     [ 0, -2, 3, 0, 0, 0, 1, 0 ],
    >     [ -8, 2, 3, 0, 1, 1, 0, 0 ],
    >     [ 8, -2, 3, 0, 1, 1, 0, 0 ],
    >     [ -2, 18, 3, 0, 1, 0, 0, 2 ],
    >     [ 10, -6, 3, 0, 1, 1, 1, 0 ],
    >     [ 2, -18, 3, 0, 1, 0, 0, 2 ],
    >     [ -10, 6, 3, 0, 1, 1, 1, 0 ],
    >     [ -14, -2, 3, 2, 0, 1, 1, 0 ],
    >     [ 1480, 334, 3, 2, 6, 1, 2, 0 ],
    >     [ 62, -10, 3, 3, 0, 1, 0, 0 ],
    >     [ 86, 50, 3, 5, 3, 1, 0, 0 ],
    >     [ -2146, 410, 3, 16, 0, 1, 0, 0 ],
    >     [ 14, 2, 3, 2, 0, 1, 1, 0 ],
    >     [ -1480, -334, 3, 2, 6, 1, 2, 0 ],
    >     [ -62, 10, 3, 3, 0, 1, 0, 0 ],
    >     [ -86, -50, 3, 5, 3, 1, 0, 0 ],
    >     [ 2146, -410, 3, 16, 0, 1, 0, 0 ],
    >     [ -52, 2, 3, 3, 1, 1, 1, 0 ],
    >     [ 52, -2, 3, 3, 1, 1, 1, 0 ],
    >     [ 8, 2, 3, 2, 0, 0, 0, 1 ],
    >     [ -32, -38, 3, 2, 0, 3, 1, 0 ],
    >     [ -10, 2, 3, 2, 1, 0, 0, 0 ],
    >     [ -34, -94, 3, 2, 1, 1, 3, 0 ],
    >     [ -86, 34, 3, 2, 1, 2, 0, 1 ],
    >     [ 2, 2, 3, 3, 0, 0, 0, 0 ],
    >     [ 8, -10, 3, 3, 1, 0, 1, 0 ],
    >     [ 3098, -502, 3, 4, 1, 0, 0, 2 ],
    >     [ 16, -2, 3, 5, 0, 0, 0, 0 ],
    >     [ 310, 22, 3, 5, 1, 0, 2, 1 ],
    >     [ 50, 14, 3, 6, 0, 2, 0, 0 ],
    >     [ 272, 290, 3, 6, 1, 1, 0, 2 ],
    >     [ 2152, -350, 3, 9, 0, 0, 1, 1 ],
    >     [ -386, 514, 3, 11, 2, 0, 0, 1 ],
    >     [ -8, -2, 3, 2, 0, 0, 0, 1 ],
    >     [ 32, 38, 3, 2, 0, 3, 1, 0 ],
    >     [ 10, -2, 3, 2, 1, 0, 0, 0 ],
    >     [ 34, 94, 3, 2, 1, 1, 3, 0 ],
    >     [ 86, -34, 3, 2, 1, 2, 0, 1 ],
    >     [ -2, -2, 3, 3, 0, 0, 0, 0 ],
    >     [ -8, 10, 3, 3, 1, 0, 1, 0 ],
    >     [ -3098, 502, 3, 4, 1, 0, 0, 2 ],
    >     [ -16, 2, 3, 5, 0, 0, 0, 0 ],
    >     [ -310, -22, 3, 5, 1, 0, 2, 1 ],
    >     [ -50, -14, 3, 6, 0, 2, 0, 0 ],
    >     [ -272, -290, 3, 6, 1, 1, 0, 2 ],
    >     [ -2152, 350, 3, 9, 0, 0, 1, 1 ],
    >     [ 386, -514, 3, 11, 2, 0, 0, 1 ],
    >     [ 260, -46, 3, 4, 0, 1, 1, 1 ],
    >     [ -260, 46, 3, 4, 0, 1, 1, 1 ],
    >     [ 33556, -4154, 3, 7, 1, 1, 5, 1 ],
    >     [ -33556, 4154, 3, 7, 1, 1, 5, 1 ],
    >     [ -68, 10, 3, 2, 1, 2, 0, 0 ],
    >     [ 68, -10, 3, 2, 1, 2, 0, 0 ],
    >     [ -26, 10, 3, 4, 1, 0, 1, 0 ],
    >     [ -6994, 1178, 3, 5, 1, 5, 1, 0 ],
    >     [ 494, 566, 3, 12, 1, 1, 1, 0 ],
    >     [ 26, -10, 3, 4, 1, 0, 1, 0 ],
    >     [ 6994, -1178, 3, 5, 1, 5, 1, 0 ],
    >     [ -494, -566, 3, 12, 1, 1, 1, 0 ]
    > ];
    > Soln4:= PostProcessingTM(C,N,p,alpha_1beta_1,Soln,RSolnFile);
    > Soln4;
    [
        [ -18, -23, 0, 0, 2, 1, 1, 1 ],
        [ -1, 9, 0, 0, 1, 0, 0, 2 ],
        [ 1, -9, 0, 0, 1, 0, 0, 2 ],
        [ 18, 23, 0, 0, 2, 1, 1, 1 ]
    ]

*/

function PostProcessingTM_GivenConductor(C,N,p,alpha_1beta_1,Soln,RSolnFile)
    
    alpha_1:= alpha_1beta_1[1];
    beta_1:= alpha_1beta_1[2];
    R<x,y>:= PolynomialRing(Integers(),2); 
    F:= C[1]*x^3 + C[2]*x^2*y + C[3]*x*y^2 + C[4]*y^3;  // generates the relevant Thue-Mahler polynomial, F
    DF:= -27*C[1]^2*C[4]^2 + C[2]^2*C[3]^2 + 18*C[1]*C[2]*C[3]*C[4] - 4*C[1]*C[3]^3 - 4*C[2]^3*C[4];    // computes the discriminant of F
    DF_fac:= [f[1] : f in Factorization(DF)];           // computes the prime factors of the discriminant of F
    
    Modp:= [];  // stores the prime factors corresponding to local obstructions; ie. if v_q(N) = 1, for q >= 3, q|( D_F*F(u,v) )
    for q in p do
        if (q ge 3) and (Valuation(N,q) eq 1) then
            Append(~Modp, q);
        end if;
    end for;
    
    Soln2:= [];         // stores only those solutions [u,v] corresponding to the condition that if v_q(N) = 1, for q >= 3, q|( D_F*F(u,v) )
    for s in Soln do
        F_uv:= Integers()! ( Evaluate(F,[s[1],s[2]]));  // computes F(u,v)
        if ( DF*(F_uv) mod &*Modp ) eq 0 then           // verifies whether ( D_F*F(u,v) ) = 0 mod Modp
            Append(~Soln2, s);
        end if;
    end for;

    alpha_10s:= [a[1] : a in alpha_1 | a[2] eq 0];      // collects only those values of alpha_1, a, for which Valuation(F(u,v),2) == a
    alpha_11s:= [a[1] : a in alpha_1 | a[2] eq 1];      // collects only those values of alpha_1, a, for which Valuation(F(u,v),2) >= a
    beta_10s:= [b[1] : b in beta_1 | b[2] eq 0];        // collects only those values of beta_1, b, for which Valuation(F(u,v),3) == b
    beta_11s:= [b[1] : b in beta_1 | b[2] eq 1];        // collects only those values of beta_1, b, for which Valuation(F(u,v),3) >= b
    
    Soln3:= [];         // stores only those solutions [u,v] corresponding to the conditions on alpha_1, beta_1, as in the Theorem (see BeGhRe)
    for s in Soln2 do
        F_uv:= Integers()! ( Evaluate(F,[s[1],s[2]]) );         // computes F(u,v)
        StillProcessing2:= true;        // determines whether or not still processing alpha_1
        StillProcessing3:= true;        // determines whether or not still processing beta_1
        
        i:= 1;
        while (i le #alpha_10s) do      // runs through values of alpha_1, a, corresponding to Valuation(F(u,v),2) == a
            a:= alpha_10s[i];
            
            if ( Valuation(F_uv,2) eq a ) then 
                StillProcessing2:= false;       // sets StillProcessing2 to false; [u,v] is a viable solution
                i:= 2*(#alpha_10s);     // sets i:= 2*(#alpha_10s) to exist while loop on alpha_10s      
            else
                i:= i + 1;
            end if;
        end while;
        
        if (StillProcessing2 eq true) then      // if StillProcessing2 == true, [u,v] is not a viable solution as per alpha_10s
                                                    // [u,v] may be a viable solution as per alpha_11s
            i:= 1;
            while (i le #alpha_11s) do  // runs through values of alpha_1, a, corresponding to Valuation(F(u,v),2) >= a
                a:= alpha_11s[i];
                
                if ( Valuation(F_uv,2) ge a ) then 
                    StillProcessing2:= false;   // sets StillProcessing2 to false; [u,v] is a viable solution
                    i:= 2*(#alpha_11s);         // sets i:= 2*(#alpha_11s) to exist while loop on alpha_11s  
                else
                    i:= i + 1;
                end if;
            end while;
        end if;
            
        if (StillProcessing2 eq false) then     // if StillProcessing2 eq false, [u,v] is a viable solution as per alpha_1
                                                    // verifies if [u,v] is also a viable solution as per beta_1
            j:= 1;
            while (j le #beta_10s) do   // runs through values of beta_1, b, corresponding to Valuation(F(u,v),3) == b
                b:= beta_10s[j];
            
                if ( Valuation(F_uv,3) eq b ) then 
                    StillProcessing3:= false;   // sets StillProcessing3 to false; [u,v] is a viable solution
                    j:= 2*(#beta_10s);  // sets j:= 2*(#beta_10s) to exist while loop on beta_10s 
                else
                    j:= j + 1;
                end if;
            end while;
        
            if (StillProcessing3 eq true) then  // if StillProcessing3 == true, [u,v] is not a viable solution as per beta_10s
                                                    // [u,v] may be a viable solution as per beta_11s 
                j:= 1;
                while (j le #beta_11s) do       // runs through values of beta_1, b, corresponding to Valuation(F(u,v),3) >= b
                    b:= beta_11s[j];
                
                    if ( Valuation(F_uv,3) ge b ) then 
                        StillProcessing3:= false;       // sets StillProcessing3 to false; [u,v] is a viable solution
                        j:= 2*(#beta_11s);      // sets j:= 2*(#beta_11s) to exist while loop on beta_11s
                    else
                        j:= j + 1;
                    end if;
                end while;
            end if;
        end if;
            
        if (StillProcessing2 eq false) and (StillProcessing3 eq false) then     // [u,v] is a viable solution as per alpha_1 and beta_1
            Append(~Soln3, s);
        end if;
    end for;
        
    
    Soln4:= [];         // stores only those solutions [u,v] such that GCD(u,v) = 1
    for s in Soln3 do
        if GCD(s[1],s[2]) eq 1 then
            Append(~Soln4,s);
        end if;
    end for;
    
    Sort(~Soln4);
    
    fprintf RSolnFile, "Relevant solutions corresponding to conductor N = %o, primes p:= %o, and form C:=%o \n ", N, p, C;
    fprintf RSolnFile, "=========================================================================== \n";
    
    fprintf RSolnFile, "Soln:= %o", Soln4;
    
    return Soln4;
end function;

