# -----------------------------------
#
#    Algorithm v1.1
#    Nonmaximal Order Constructor
#
#    Author: Jordan Wiebe
#    Date:   10/18/2018
#
#    Copyright Jordan Wiebe 2018
#
# -----------------------------------



# The following algorithm takes as inputs the discriminant of a definite rational quaternion algebra 
# and an admissible level for an order, and constructs an order with the provided level in a 
# quaternion algebra with the provided discriminant.  This algorithm is the implementation of the 
# work in "Constructing non-maximal orders in quaternion algebras", which can be accessed via the 
# author's website at http://math.ou.edu/~jwiebe, or on the arXiv at https://arxiv.org/abs/1810.05249
#
# Any questions or comments may be directed to the author via email using jwiebe [at] ou [dot] edu





def order(disc,lev):
    x = disc/2;o
    
    n1 = 1;
    n2 = 1;
    m1 = 1;
    m2 = 1;
    for a in disc.prime_factors():
        even_val = lev.valuation(a)/2;
        if even_val.is_integer():
            n2*=a;
        else:
            n1*=a;
    
    lev2=lev.prime_to_m_part(disc);
    for b in lev2.prime_factors():
        even_val = lev.valuation(b)/2;
        if even_val.is_integer():
            m2*=b;
        else:
            m1*=b;
    
    if not(x.is_integer()):
        y = lev.valuation(2)/2;
        if y.is_integer():
            #Case1
            prod = n1*n2*m1;
            if prod.mod(4)==3:
                #print "Case 1a";
                a = -n1*n2*m1;
                Q = [1];
                P = [8];
                f=1;
                g=1;
                h=1;
                for p in (n1).prime_factors():
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*kronecker(-1,p)];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^((lev.valuation(p)-1)/2);
                for p in (n2).prime_factors():
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*kronecker(-1,p)];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^(lev.valuation(p)/2-1);
                    h*=p;
                for p in m1.prime_factors():
                    M1 = [x for x in range(1,p) if kronecker(x,p)==kronecker(-1,p)];
                    Q.append(M1[0]);
                    P.append(p);
                    g*=p^((lev.valuation(p)-1)/2);
                for p in m2.prime_factors():
                    if kronecker(a,p)==1:
                        f*=p^(lev.valuation(p)/2);
                    else:
                        g*=p^(lev.valuation(p)/2);
                q = CRT(Q,P);
                while not(q.is_prime()):
                    q+=8*n1*n2*m1;
                b=-1*q;
                
                #Quaternion algebra stuff
                A.<i,j,k>=QuaternionAlgebra(a,b);
                
                if A.discriminant()!=disc:
                    print "Error.  The algebra created has a discriminant which does not match the discriminant provided.  Are you sure you've chosen a valid discriminant?";
                    return;

                x = Mod(a,q).sqrt();
                
                #Hermite normal form stuff
                d,u,v = xgcd(q,2*x);
                y = u*q+v*(-2);
                
                while (y-2*q)>0:
                    y-=2*q;
                
                O=A.quaternion_order([1,g*(1+i)/2,h*f*g*(j+y*k)/(2*q),f*k]);
                if O.discriminant()==lev:
                    return O;
                else:
                    print "Error.  The order created has level which does not match the level provided.  Are you sure you've chosen a viable level for the discriminant provided?";
                    return;
            else:
                #print "Case 1b";
                b = -n1*n2*m1;
                if n2.mod(4)==1:
                    Q = [3];
                else:
                    Q = [1];
                P = [4];
                f=1;
                g=1;
                h=1;
                
                for p in (n1).prime_factors():
                    l = kronecker(-1,p);
                    for q in (n2).prime_factors():
                        l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^((lev.valuation(p)-1)/2);
                for p in (n2).prime_factors():
                    l = kronecker(-1,p);
                    for q in (n1*n2*m1).prime_factors():
                        if q in (n1*m1).prime_factors():
                            l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^(lev.valuation(p)/2-1);
                for p in m1.prime_factors():
                    l = kronecker(-1,p);
                    for q in (n2).prime_factors():
                        l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==l];
                    Q.append(N[0]);
                    P.append(p);
                    g*=p^((lev.valuation(p)-1)/2);
                for p in m2.prime_factors():
                    if p!=2:
                        if kronecker(a,p)==1:
                            f*=p^(lev.valuation(p)/2);
                        else:
                            g*=p^(lev.valuation(p)/2);
                
                q = CRT(Q,P);
                while not(q.is_prime()):
                    q+=4*n1*n2*m1;
                a=-1*q*n2;
                
                if 2 in m2.prime_factors():
                    if a.mod(8)==1:
                        f*=2^(lev.valuation(2)/2);
                    elif a.mod(8)==5:
                        g*=2^(lev.valuation(2)/2);
                
                #Quaternion algebra stuff
                A.<i,j,k>=QuaternionAlgebra(a,b);
                
                if A.discriminant()!=disc:
                    print "Error.  The algebra created has a discriminant which does not match the discriminant provided.  Are you sure you've chosen a valid discriminant?";
                    return;
                
                x = Mod(b,q).sqrt();
                x = Integer(x);
                
                #Hermite normal form stuff
                d,y,z = xgcd(-1*q,2*x);
                z*=2;
                c = 2*z;

                while (c-2*q)>0:
                    c-=2*q;
                
                O=A.quaternion_order([(q+i+z*k)/(2*q),(2*i+c*k)/(2*q),f*g*(h*j+k)/2,f*g*k]);
                
                if O.discriminant()==lev:
                    return O;
                else:
                    print "Error.  The order created has level which does not match the level provided.  Are you sure you've chosen a viable level for the discriminant provided?";
                    return;
        else:
            #Case2
            b = -n1*n2*m1;
            
            # 2 Stuff
            Q=[];
            if n2.mod(8)==1:
                Q.append(7);
            elif n2.mod(8)==3:
                Q.append(5);
            elif n2.mod(8)==5:
                Q.append(3);
            elif n2.mod(8)==7:
                Q.append(1);
            P = [8];
            f=1;
            g=1;
            h=1;
            
            for p in (n1).prime_factors():
                l = kronecker(-1,p);
                for q in (n2).prime_factors():
                    l*=kronecker(q,p);
                N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                Q.append(N[0]);
                P.append(p);
                f*=p^((lev.valuation(p)-1)/2);
            for p in (n2).prime_factors():
                l = kronecker(-1,p);
                for q in (n1*n2*m1).prime_factors():
                    if q in (n1*m1).prime_factors():
                        l*=kronecker(q,p);
                N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                Q.append(N[0]);
                P.append(p);
                f*=p^(lev.valuation(p)/2-1);
            for p in m1.prime_factors():
                if p!=2:
                    l = kronecker(-1,p);
                    for q in (n2).prime_factors():
                        l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==l];
                    Q.append(N[0]);
                    P.append(p);
                    g*=p^((lev.valuation(p)-1)/2);
                else:
                    g*=p^((lev.valuation(p)-1)/2);
            for p in m2.prime_factors():
                if p!=2:
                    if kronecker(a,p)==1:
                        f*=p^(lev.valuation(p)/2);
                    else:
                        g*=p^(lev.valuation(p)/2);
            
            q = CRT(Q,P);
            while not(q.is_prime()):
                q+=4*n1*n2*m1;
            a=-1*q*n2;
            
            if 2 in m2.prime_factors():
                if a.mod(8)==1:
                    f*=2^(lev.valuation(2)/2);
                elif a.mod(8)==5:
                    g*=2^(lev.valuation(2)/2);
            
            #Quaternion algebra stuff
            A.<i,j,k>=QuaternionAlgebra(a,b);
            
            if A.discriminant()!=disc:
                print "Error.  The algebra created has a discriminant which does not match the discriminant provided.  Are you sure you've chosen a valid discriminant?";
                return;
            
            x = Mod(b,q).sqrt();
            x = Integer(x);
            
            #Hermite normal form stuff
            d,y,z = xgcd(-1*q,2*x);
            z*=2;
            c = 2*z;

            while (c-2*q)>0:
                c-=2*q;
            
            O=A.quaternion_order([(q+i+z*k)/(2*q),(2*i+c*k)/(2*q),f*g*(h*j+k)/2,f*g*k]);
            
            if O.discriminant()==lev:
                return O;
            else:
                print "Error.  The order created has level which does not match the level provided.  Are you sure you've chosen a viable level for the discriminant provided?";
                return;
    else:
        #Case3
        
        f=1;
        g=1;
        h=1;
        Q = [];
        P = [];
            
        if lev.valuation(2)==1:
            # 2M
            b = -n1*n2*m1;
            a = -n2;
            P.append(8);
        elif lev.valuation(2)==2:
            # 4M
            prod = Integer(n1*n2*m1/2);
                
            P.append(4);
            
            if prod.mod(4)==1:
                a = -Integer(n2/2);
                b = -prod;
            else:
                # It's possible to obtain an order of the correct level, but a systematic solution is not viable
                # Better to go in a case-by-case basis here.  For a generic order of level 4M, we need a,b==3 mod 4
                print "This cannot be done systematically! Recommend working by hand to choose a,b";
                return;
            
        elif lev.valuation(2)==3:
            # 8M
            b = -n1*n2*m1;
            a = -n2;
            P.append(8);
        else:
            b = -n1*n2*m1;
            a = -n2;
            P.append(8);
        
        #2 Stuff
        if lev.valuation(2)==1:
            # need to make a==5 mod 8, so q*n2==3 mod 8
            if n2.mod(8)==1:
                Q.append(3);
            elif n2.mod(8)==3:
                Q.append(1);
            elif n2.mod(8)==5:
                Q.append(7);
            else:
                Q.append(5);
        elif lev.valuation(2)==2:
            if Integer(n2/2).mod(4)==1:
                Q.append(1);
            else:
                Q.append(3);
        elif lev.valuation(2)==3:
            # need to make a==5 mod 8, so q*n2==3 mod 8
            if n2.mod(8)==1:
                Q.append(3);
            elif n2.mod(8)==3:
                Q.append(1);
            elif n2.mod(8)==5:
                Q.append(7);
            else:
                Q.append(5);
        else:
            if is_even(n2):
                b1=Integer(b/2);
                a1=Integer(a/2);
                if b1.mod(8)==1:
                    # need a'==3,5 mod 8
                    if a1.mod(8)==1:
                        Q.append(3);
                    elif a1.mod(8)==3:
                        Q.append(1);
                    elif a1.mod(8)==5:
                        Q.append(1);
                    else:
                        Q.append(3);
                elif b1.mod(8)==3:
                    # need a'==1,3 mod 8
                    if a1.mod(8)==1:
                        Q.append(1);
                    elif a1.mod(8)==3:
                        Q.append(1);
                    elif a1.mod(8)==5:
                        Q.append(5);
                    else:
                        Q.append(5);
                elif b1.mod(8)==5:
                    # need a'==1,7 mod 8
                    if a1.mod(8)==1:
                        Q.append(1);
                    elif a1.mod(8)==3:
                        Q.append(3);
                    elif a1.mod(8)==5:
                        Q.append(3);
                    else:
                        Q.append(1);
                else:
                    # need a'==5,7 mod 8
                    if a1.mod(8)==1:
                        Q.append(5);
                    elif a1.mod(8)==3:
                        Q.append(5);
                    elif a1.mod(8)==5:
                        Q.append(1);
                    else:
                        Q.append(1);
                    
            else:
                b1=Integer(b/2);
                if b1.mod(4)==1:
                    # need a==3,5 mod 8
                    # choosing a==3 mod 8 for matching order to previous case
                    if a.mod(8)==1:
                        Q.append(3);
                    elif a.mod(8)==3:
                        Q.append(1);
                    elif a.mod(8)==5:
                        Q.append(7);
                    else:
                        Q.append(5);
                else:
                    # need a==5,7 mod 8
                    # choosing a==7 mod 8 for matching order to previous case
                    if a.mod(8)==1:
                        Q.append(7);
                    elif a.mod(8)==3:
                        Q.append(5);
                    elif a.mod(8)==5:
                        Q.append(3);
                    else:
                        Q.append(1);
                    
                
                    
        
        if lev.valuation(2)!=2:
            for p in (n1).prime_factors():
                if p!=2:
                    l = kronecker(-1,p);
                    for q in (n2).prime_factors():
                        l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^((lev.valuation(p)-1)/2);
                else:
                    if lev.valuation(2)==1:
                        # 2M
                        f*=p^((lev.valuation(p)-1)/2);
                    elif lev.valuation(2)==2:
                        # 4M
                        f*=p^((lev.valuation(p)-1)/2);
                    elif lev.valuation(2)==3:
                        # 8M
                        f*=p^((lev.valuation(p)-3)/2);
                    else:
                        f*=p^((lev.valuation(p)-3)/2);
            for p in (n2).prime_factors():
                if p!=2:
                    l = kronecker(-1,p);
                    for q in (n1*n2*m1).prime_factors():
                        if q in (n1*m1).prime_factors():
                            l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^(lev.valuation(p)/2-1);
                else:
                    if lev.valuation(2)==1:
                        # 2M
                        f*=p^(lev.valuation(p)/2-1);
                        h*=p;
                    elif lev.valuation(2)==2:
                        # 4M
                        f*=p^(lev.valuation(p)/2-1);
                        h*=p;
                    elif lev.valuation(2)==3:
                        # 8M
                        f*=p^(lev.valuation(p)/2-1);
                        h*=p;
                    else:
                        f*=p^(lev.valuation(p)/2-2);
            for p in m1.prime_factors():
                l = kronecker(-1,p);
                for q in (n2).prime_factors():
                    l*=kronecker(q,p);
                N = [x for x in range(1,p) if kronecker(x,p)==l];
                Q.append(N[0]);
                P.append(p);
                g*=p^((lev.valuation(p)-1)/2);
            for p in m2.prime_factors():
                if kronecker(a,p)==1:
                    f*=p^(lev.valuation(p)/2);
                else:
                    g*=p^(lev.valuation(p)/2);
                    
        else:
            # 4M
            
            for p in (n1).prime_factors():
                l = kronecker(-1,p);
                for q in (Integer(n2/2)).prime_factors():
                    l*=kronecker(q,p);
                N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                Q.append(N[0]);
                P.append(p);
                f*=p^((lev.valuation(p)-1)/2);
            for p in (n2).prime_factors():
                if p!=2:
                    l = kronecker(-1,p);
                    for q in (n1*n2*m1).prime_factors():
                        if q in (n1*m1).prime_factors():
                            l*=kronecker(q,p);
                    N = [x for x in range(1,p) if kronecker(x,p)==-1*l];
                    Q.append(N[0]);
                    P.append(p);
                    f*=p^(lev.valuation(p)/2-1);
            for p in m1.prime_factors():
                l = kronecker(-1,p);
                for q in (Integer(n2/2)).prime_factors():
                    l*=kronecker(q,p);
                N = [x for x in range(1,p) if kronecker(x,p)==l];
                Q.append(N[0]);
                P.append(p);
                g*=p^((lev.valuation(p)-1)/2);
            for p in m2.prime_factors():
                if kronecker(a,p)==1:
                    f*=p^(lev.valuation(p)/2);
                else:
                    g*=p^(lev.valuation(p)/2);
        
        q = CRT(Q,P);
        
        while not(q.is_prime()):
            q+=8*n1*n2*m1;
        
        a*=q;
        
        #Quaternion algebra stuff
        A.<i,j,k>=QuaternionAlgebra(a,b);
        
        if A.discriminant()!=disc:
            print "Error.  The algebra created has a discriminant which does not match the discriminant provided.  Are you sure you've chosen a valid discriminant?";
            return;
        
        x=Integer(Mod(b,q).sqrt());
        
        #Hermite normal form stuff
        if is_even(n2):
            d,s,t = xgcd(q,h*x);
        else:
            d,s,t = xgcd(q,h*x);
            
        
        if lev.valuation(2)==1:
            d,y,z = xgcd(-1*q,2*x);
            z*=2;
            c = 2*z;
            while (c-2*q)>0:
                c-=2*q;
                
            O=A.quaternion_order([(q+i+z*k)/(2*q),(2*i+c*k)/(2*q),f*g*(h*j+k)/2,f*g*k]);
        elif lev.valuation(2)==2:
            O=A.quaternion_order([1,g*h*(i+t*k)/q,f*g*j,f*k]);
        elif is_even(n2):
            O=A.quaternion_order([1,g*h*(i+t*k)/q,f*g*j,f*k]);
        else:
            O=A.quaternion_order([1,g*h*(i+t*k)/q,f*g*j,f*k]);
        
        if O.discriminant()==lev:
            return O;
        else:
            print "Error.  The order created has level which does not match the level provided.  Are you sure you've chosen a viable level for the discriminant provided?";
            return;
        
        #-----------------------
        
    return;



# Sample Order Calculation
#order(2*3*11,2*3^2*11^5*5);



# Testing function
def test(disc_limit,lev_limit):
    for x in range(2,disc_limit):
        # x is our discriminant
        P = Integer(x).prime_factors();
        if not(is_even(len(P))) and Integer(x).is_squarefree():
            y=1;
            while x*y<lev_limit:
                x = Integer(x);
                y = Integer(y);
                order(x,x*y);
                y+=1;
    return;