# -*- coding: utf-8 -*-
"""
The GeometricXL Algorithm.

Gröbner basis algorithms and related methods like XL are algebraic in
nature. In particular, their complexity is not invariant under a
linear change of coordinates. As an example consider Cyclic-6::

    sage: P.<a,b,c,d,e,f,h> = PolynomialRing(GF(32003))
    sage: I = sage.rings.ideal.Cyclic(P,6).homogenize(h)
    sage: J = Ideal(I.groebner_basis())

The generators of ``J`` form a Gröbner basis and we can use this
property to find a common root for these generators. Now, consider the
same equations but permute the variables in the ring::

    sage: P.<a,b,c,d,e,f,h> = PolynomialRing(GF(32003),order='lex')
    sage: I = sage.rings.ideal.Cyclic(P,6).homogenize(h)
    sage: J = Ideal(I.groebner_basis())
    sage: R = PolynomialRing(GF(32003),P.ngens(),list(reversed(P.variable_names())),order='lex')
    sage: H = Ideal([R(f) for f in J.gens()])

The generators of ``H`` do not form a Gröbner basis in ``R`` which is
``P`` with its variables reversed. If we are only trying to solve a
system of equations choosing the right permutation of variables might
make a significant impact on the performance of our Gröbner basis algorithm::

    sage: t = cputime()
    sage: gb = H.groebner_basis('libsingular:std')
    sage: gb[-1].degree()
    19
    sage: cputime(t) # output random-ish
    25.36...

While in this example it is easy to see which variable permutation is
the cheapest one, this is not necessarily the case in general. The
GeometricXL algorithm [MP07]_ is invariant under any linear change of
coordinates and has the following property:

    Let ``D`` be the degree reached by the algorithm XL to solve a
    given system of equations under the optimal linear change of
    coordinates. Then GeometricXL will also solve this system of
    equations for the degree ``D``, without applying this optimal
    linear change of coordinates first. 

The above behaviour holds under two assumptions:

 * the characteristic of the base field ``K`` is bigger than ``D``
 * the system of equations has one over "very few" solution.

To demonstrate this behaviour, we use a synthetic benchmark which is a
Gröbner basis under a linear change of coordinates::

    sage: e,h = random_example(n=6)

``e`` is the original easy system while ``h`` is the "rotated"
system::

    sage: e.basis_is_groebner()
    True

    sage: max([f.total_degree() for f in e.gens()])
    2

    sage: h.basis_is_groebner()
    False

    sage: max([f.total_degree() for f in h.gens()])
    2

GeometricXL recovers linear factors and thus candidates for common
roots at ``D=2``::

    sage: hH = h.homogenize()
    sage: f = GeometricXL(hH, D=2); f.factor(False)
    0.0...s -- 1. D: 2
    ...
    (-2684) 
    * (-1056*x5 - 2964*x4 - 177*x3 + 6206*x2 + 376*x1 + 6257*x0 + h) 
    * (2957*x5 - 792*x4 - 4323*x3 - 14408*x2 - 2750*x1 - 8823*x0 + h)

While any Gröbner basis algorithm would have to reach at least degree 64::

    sage: gb = h.groebner_basis('libsingular:slimgb')
    sage: gb[-1].degree()
    64

AUTHORS:

- Martin Albrecht - initial, ad-hoc implementation

.. note::

  This implementation is very ad-hoc and not robust by any stretch.

REFERENCES:

.. [MP07] S. Murphy and M.B. Paterson; *A Geometric View of
  Cryptographic Equation Solving*; Journal of Mathematical Cryptology,
  Vol. 2; pages 63-107; 2008. A version is available as Departmental
  Technical Report RHUL-MA-2007-4 at
  http://www.ma.rhul.ac.uk/static/techrep/2007/RHUL-MA-2007-4.pdf
"""

from sage.all import *

def random_minors(A, k, count):
    """
    Return a list of ``count`` elements containing ``k``-minors of
    ``A``.

    Let ``A`` be an ``m x n`` matrix and k an integer with ``0 < k <=
    m``, and ``k <= n``. A ``k x k`` minor of ``A`` is the determinant
    of a ``k x k`` matrix obtained from ``A`` by deleting ``m - k``
    rows and ``n - k`` columns.

    INPUT:

    - ``k`` - integer

    - ``count`` - the number of elements returned

    EXAMPLE::

        sage: A = Matrix(ZZ,2,3,[1,2,3,4,5,6]); A
        [1 2 3]
        [4 5 6]
        sage: random_minors(A, 2, 3)
        set([-6, -3])
    """
    all_rows = range(A.nrows())
    all_cols = range(A.ncols())
    m = set()

    total = (binomial(A.nrows(),k) * binomial(A.ncols(),k))
    ratio = ZZ(count)/total

    if count > total/2:
        for rows in combinations_iterator(all_rows,k):
            for cols in combinations_iterator(all_cols,k):
                if random() <= ratio:
                    m.add(A.matrix_from_rows_and_columns(rows,cols).determinant())
        return m
    else:
        C_r, C_c = Combinations(all_rows, k), Combinations(all_cols, k)
        l_r, l_c = C_r.cardinality(), C_c.cardinality()

        while len(m) != count:
            r,c = randint(0,l_r-1), randint(0,l_c-1)
            rows, cols = C_r.unrank(r), C_c.unrank(c)   
            B = A.matrix_from_rows_and_columns(rows,cols)
            B.set_immutable()
            m.add(B)

        return [B.determinant() for B in m]

def random_minor(A, k):
    while True:
        rows, cols = set([randint(0,A.nrows()-1) for _ in range(k)]), set([randint(0,A.ncols()-1) for _ in range(k)])
        if len(rows) != k or len(cols) != k:
            continue
        rows, cols =  sorted(rows), sorted(cols)
        B = A.matrix_from_rows_and_columns(rows,cols)
        return B.determinant()

def min_rank_system(A, r, m):
    F = []
    M = set()
    old = (0,0)
    i = 0

    for i in range(m):
        minor = random_minor(A, k=r+1)
        if minor == 0: 
            continue
        m = minor.monomials()
        M = M.union(m)
        F.append(minor)

        if len(M) > len(F):
            continue
        F = list(Ideal(F).interreduced_basis())
        if F[-1].nvariables() <= 2:
            break
        if old != (len(F), len(M)):
            print "    |F|: %4d |M|: %4d"%(len(F), len(M))
            old = len(F), len(M)
        
    print "    |F|: %4d |M|: %4d"%(len(F), len(M))
    return mq.MPolynomialSystem(F)

def C_f(f,D):
    """
    Return the partial derivative matrix `C_f^D` for the polynomial
    ``f`` and the degree ``D``.

    `C_f^D` is the coefficient matrix of a set of polynomials
    generated by deriving `f` w.r.t. to every monomial of degree `D`.

    The monomials/rows of the coefficient matrix are ordered w.r.t. to
    the term ordering of the parent ring in descending order.

    INPUT:

    - ``f`` - polynomial

    - ``D`` - degree

    OUTPUT:

        A,v such that A is a coefficient matrix and v a monomial vector.

    EXAMPLE::

        sage: P.<x0,x1,x2> = PolynomialRing(GF(37),3,order='lex')
        sage: f = (-12) * (-15*x0 + 9*x1 + x2) * (-9*x0^2 + 10*x0*x1 - 12*x0*x2 - 14*x1^2 - 18*x1*x2 + x2^2)
        sage: A,v = C_f(f,1); A
        [24 31  3 26  8 28]
        [34 15  8 22  6 34]
        [20  8 19  3 31  1]

        sage: A,v = C_f(f,2); A
        [11 31  3]
        [31 15  8]
        [ 3  8 19]
        [15  7  6]
        [ 8  6 31]
        [19 31  2]
    """
    P = f.parent()
    gens = P.gens()
    monomials = sum([mul(map(pow, gens, exps)) for exps in am(D, len(gens))],P(0))
    partial_diffs = []
 
    for e in monomials.exponents():
        fbar = f
        for i in range(len(e)):
            for j in range(e[i]):
                fbar = fbar.derivative(gens[i])
        partial_diffs.append(fbar)

    A,v = mq.MPolynomialSystem(P,partial_diffs).coefficient_matrix()
    return A,v

def am(D,size):
    """
    Return a list of a exponent tuples of length ``size`` such that
    the degree of the associated monomial is ``D``.

    INPUT:
    
    - ``D`` - degree (must be > 0)
    - ``size`` - length of exponent tuples (must be > 0)

    EXAMPLE:
        sage: am(2,3)
        [(2, 0, 0), (1, 1, 0), (1, 0, 1), (0, 2, 0), (0, 1, 1), (0, 0, 2)]
    """
    res = []
    for d2 in range(D+1):
        d = D-d2
        if size>1:
            for rest in am( d2 , size-1):
                res.append( (d,) + rest )
        else:
            return [ (d,) ]
    return res

class XLDegreeError(Exception):
    pass

def GeometricXL(F,D):
    """
    The GeometricXL algorithm as presented in [MP07]_.

    INPUT:

    - ``F`` - an ``MPolynomialSystem`` or ideal

    - ``D`` - XL degree (> 0)

    EXAMPLE:

    We compute the example from the paper::

        sage: P.<x0,x1,x2> = PolynomialRing(GF(37),3,order='lex')
	sage: x0 > x1 > x2
	True
	sage: f1 = 15*x0^2 + x1^2 + 5*x1*x2
	sage: f2 = 23*x0^2 + x2^2 + 9*x1*x2
	sage: I1 = P * [f1,f2]
	sage: f = GeometricXL(I1, D=2); f
        0.0...s -- 1. D: 2
        ...
        x1^2 + 12*x1*x2 + 9*x2^2

        sage: f.factor(False)
        (9) * (8*x1 + x2) * (18*x1 + x2)

	sage: A = Matrix(P,3,3,[2,26,10,26,4,13,33,21,2]); A
	[  2 -11  10]
	[-11   4  13]
	[ -4 -16   2]
	sage: varmap = (A * Matrix(P,3,1,[x0,x1,x2])).list(); varmap
	[2*x0 - 11*x1 + 10*x2, -11*x0 + 4*x1 + 13*x2, -4*x0 - 16*x1 + 2*x2]
	sage: I2 = P * (f1(*varmap), f2(*varmap))
	sage: f = GeometricXL(I2, D=2); f
        0.000s -- 1. D: 2
        ...
        x0*x1 + 16*x0*x2 - 13*x1^2 + 10*x1*x2 + 10*x2^2

        sage: f.factor(False)
        (10) * (7*x1 + x2) * (9*x0 - 6*x1 + x2)
    """

    try:
	F = mq.MPolynomialSystem(F.ring(),F.gens())
    except AttributeError:
        F = mq.MPolynomialSystem(F)

    D = D if D >= max(f.degree() for f in F) else max(f.degree() for f in F)

    p = Profiler()

    ## 1. Generate the m(binomial(D-2+n)(D-2)) possible polynomials of
    ##    degree D that are formed by multiplying each of the
    ##    polynomials of the original system by some monomial of degree
    ##    D-2.
    p("1. D: %d"%(D,)); print p.print_last()

    ## 2. The degree D is required to be less than the characteristic
    ##    of the finite field F.

    P = F.ring()
    K = P.base_ring()
    gens = P.gens()


    # Note: MatrixF5 would be fine to use here instead of
    # straight-forward XL.
    L = []
    for f in F:
        d = f.total_degree()
        if d > D:
            continue
        
        if d == D:
            L.append(f)
        else:
            M = [mul(map(pow, gens, exps)) for exps in am(D - d,len(gens))]
            if M == []: M = [1]
            L.extend([m*f for m in M])
    assert([f.degree() == D for f in L])

    ## 3. Find a basis S of the linear span of all the polynomials
    ##    generated by the first step.
    p("3. |L|: %d"%(len(L),)); print p.print_last()

    S = mq.MPolynomialSystem(P,L)
    A,v = S.coefficient_matrix()
    A.echelonize()
    S = (A*v).list()

    ## 4. Calculate the matrix C_f^{D-1} of (D - 1)th partial
    ##    derivatives for each polynomial f in S.
    p("4. |S|: %d"%(len(S),)); print p.print_last()

    C = [C_f(f,D-1) for f in S]

    # the coefficent matrices may have different ncols. V is a lookup
    # to map these to a common superspace.

    V = flatten([v.list() for A,v in C] )
    V = sorted( uniq( V ), reverse=True)
    if V[-1] == 0:
	V = V[:-1]
    V = zip( V,range(len(V)) )
    V = dict( V )

    # Todo: Bug in Sage, cannot auto-coerce to this
    Q = PolynomialRing(K, len(C), 'l', order='lex')

    Cbar = []
    for A,v in C:
	if A.is_zero():
	    continue
        Abar = Matrix(Q, A.nrows(), len(V))
        for c in range(A.ncols()):
	    if v[c,0] == 0:
		continue
            cbar = V[v[c,0]]
            for r in range(A.nrows()):
                Abar[r,cbar] = A[r,c]
        Cbar.append(Abar)
    C = Cbar

    min_rank_solutions = []

    for i,Ci in enumerate(C):
	if Ci.change_ring(K).rank() <= 2:
	    w = dict(zip(Q.gens(),[0 for _ in range(len(Q.gens()))]))
	    w[Q.gen(i)] = 1 #(0,0,0,0, ..., 1, ...., 0,0,0,0) 1 at index i.
	    min_rank_solutions.append(w)

    ## 5. Find a linear combination of these partial derivative
    ##    matrices C_f^{D-1} which has rank 2 (or lower) by
    ##    considering the 3-minors or some other method.

    p("5. |min_rank_solutions|: %d"%(len(min_rank_solutions),)); print p.print_last() ; sys.stdout.flush()

    CC = sum( map(mul, zip(Q.gens(),C) ) )
    if True:
        Gbar = min_rank_system(CC, 2, 3*Q.ngens()**3).gens()
        min_rank_solutions += hom_variety(Gbar)
    else:
        # Fall back old code
        G = random_minors(CC, k=3, count=3*Q.ngens()**3) # choose a random subset which we expect to be big enough

        # Linearization
        G = mq.MPolynomialSystem(Q, G)
        A,v = G.coefficient_matrix()
        A.echelonize()
        B = A.matrix_from_rows(range(A.rank()))
        Gbar = (B*v).list()
        min_rank_solutions += hom_variety(Gbar)

    ## 6. Note that this it is not always possible to find such a
    ##    linear combination, and in this case GeometricXL fails for
    ##    degree D.
    p("6. |min_rank_solutions|: %d"%(len(min_rank_solutions),)); print p.print_last(); sys.stdout.flush()

    W = []
    for w in min_rank_solutions:
        try: # homogeneous solution
            # now choose a value != 0 for the free variable
            l0 = w.values()[0].variable(0)

            # and map the depdendent variables accordingly
            wbar = {l0:1}
            for var,val in w.iteritems():
                wbar[var] = val.subs({l0:1})
            w = wbar
        except AttributeError:
            pass
        W.append(w)
    min_rank_solutions = W

    ## 7. Using this linear combination, construct a polynomial in the
    ##    linear span of S that is known to have factors, and then
    ##    factorise this polynomial. This potentially allows the
    ##    elimination of a variable from the original system of
    ##    equations.
    p("7. |min_rank_solutions|: %d"%(len(min_rank_solutions),)); print p.print_last(); sys.stdout.flush()

    w = min_rank_solutions[0]

    f = sum( map(mul, zip( map(lambda x: K(w[x]), Q.gens()), S)) )
    return f

    # 8. The process is repeated on the new smaller system until a
    #    complete solution is found.
    p("8"); print p.print_last()

    pass

def hom_variety(T, V=None, v=None):
    if V is None: 
	V = []
    if v is None:
	v = {}

    if T == []:
	return v

    last_candidate = None
    for f in T:
	factors = f.factor(proof=False)
	for factor, _ in factors:
	    if factor.degree() == 1 and factor.nvariables() == 2:
		var = factor.lm()
		factor = factor * factor.coefficient(var)**(-1)
		val = var - factor
		vbar = copy(v)
                if v.get(var,None) == val:
                    continue
		vbar[var] = val                
		Tbar = [f.subs(vbar) for f in T]
		Tbar = [f for f in Tbar if f != 0]
		vbar = hom_variety(Tbar,V,vbar)
		if isinstance(vbar, dict) and vbar not in V:
		    V.append(vbar)
        if len(V) > 0:
            break
    return V

def random_example(K=GF(32003), n=3):
    """
    EXAMPLE::

        sage: e,h = random_example(GF(127), n=3)
        sage: e
        Ideal (-27*x0^2 - 11*x0 - 56, 
               -2*x1^2 + 33*x1*x0 - x1 + 13*x0^2 - 52*x0, 
               -24*x2^2 + 26*x2 + 8*x1*x0 - 50*x0^2 + 34) of Multivariate Polynomial Ring in x2, x1, x0 over Finite Field of size 127
        sage: h
        Ideal (-23*x2^2 - 61*x2*x1 + 13*x2*x0 + 3*x2 + 41*x1^2 + 20*x1*x0 - 54*x1 + 52*x0^2 + 24*x0 - 56, 
               10*x2^2 - 57*x2*x1 + 19*x2*x0 + 43*x2 + 53*x1^2 - 48*x1*x0 + 21*x1 + 34*x0^2 + 8*x0, 
               50*x2^2 - 14*x2*x1 + 26*x2*x0 - 58*x2 + 49*x1^2 - 33*x1*x0 + 11*x1 + 51*x0^2 + 27*x0 + 34) of Multivariate Polynomial Ring in x2, x1, x0 over Finite Field of size 127
    """
    while True:
        l = []
        for i in range(n):
            P = PolynomialRing(K, i+1, list(reversed(["x%d"%j for j in range(i+1)])), order='lex')
            l.append(P.gen(0)**2 + P.random_element())
        P = l[-1].parent()
        l =  [P(f) for f in l]
        if len(Ideal(l).variety()) > 0:
                break

    easy = Ideal(l)

    A = random_matrix(K, n, n)
    assert(A.rank() == n)

    phi = (A * Matrix(P,n,1,P.gens())).list()

    hard = P * [f(*phi) for f in easy.gens()]
    return easy, hard
            


def benchmarketing(K=GF(32003), max_n=8):
    T = []
    singular_crossover = True
    for n in range(2,max_n+1):
        e,h = random_example(K, n)
        hH = h.homogenize()
        t = cputime()
        f = GeometricXL(hH, D=2)
        F = f.factor(False)
        gt = cputime(t)
        assert(F[0][0].degree() == 1)
        assert(F[1][0].degree() == 1)

        if not singular_crossover:
            t = singular.cputime()
            gb = h.groebner_basis('singular')
            st = singular.cputime(t)
            sd = max(f.degree() for f in gb)
        else:
            st, sd = 0.0, 0
        if st > 12.0:
            singular_crossover = True


        t = magma.cputime()
        gb = h.groebner_basis('magma')
        mt = magma.cputime(t)
        md = max(f.degree() for f in gb)

        T.append([n, gt, 2, st, sd, mt, md])
        print T[-1]
        sys.stdout.flush()
    return T