Ticket #11036: trac_11036-solve-mod.patch

File trac_11036-solve-mod.patch, 9.1 KB (added by mderickx, 8 years ago)

removed a lot of bugs and refactored the patch

  • sage/symbolic/relation.py

    # HG changeset patch
    # User User D. S. McNeil <dsm054@gmail.com>
    # Date 1312300545 -3600
    # Node ID 6d3f8e75a73d6b333e5f47d00b7ca6e574edd116
    # Parent  b2c0ed2e09a453751e4b2fe092667487f530eac3
    Trac 11036: improve solve_mod performance
    
    diff --git a/sage/symbolic/relation.py b/sage/symbolic/relation.py
    a b  
    769769        sage: d[y]
    770770        8610183
    771771
     772    For cases where there are relatively few solutions and the prime
     773    factors are small, this can be efficient even if the modulus itself
     774    is large::
     775
     776        sage: sorted(solve_mod([x^2 == 41], 10^20))
     777        [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,),
     778        (45461397519473547571,), (54538602480526452429,), (61445932736758703821,),
     779        (88554067263241296179,), (95461397519473547571,)]
     780
    772781    We solve a simple equation modulo 2::
    773782   
    774783        sage: x,y = var('x,y')
     
    777786
    778787    .. warning::
    779788
    780        The current implementation splits the modulus into prime powers,
    781        then naively enumerates all possible solutions and finally combines
    782        the solution using the Chinese Remainder Theorem.
    783        The interface is good, but the algorithm is horrible if the modulus
    784        has some larger prime factors! Sage *does* have the ability to do
    785        something much faster in certain cases at least by using Groebner
    786        basis, linear algebra techniques, etc. But for a lot of toy problems
    787        this function as is might be useful. At least it establishes an interface.
     789       The current implementation splits the modulus into prime
     790       powers, then naively enumerates all possible solutions
     791       (starting modulo primes and then working up through prime
     792       powers), and finally combines the solution using the Chinese
     793       Remainder Theorem.  The interface is good, but the algorithm is
     794       very inefficient if the modulus has some larger prime factors! Sage
     795       *does* have the ability to do something much faster in certain
     796       cases at least by using Groebner basis, linear algebra
     797       techniques, etc. But for a lot of toy problems this function as
     798       is might be useful. At least it establishes an interface.
     799
     800
     801    TESTS:
     802
     803    Make sure that we short-circuit in at least some cases::
     804
     805        sage: solve_mod([2*x==1], 2*next_prime(10^50))
     806        []
     807
     808    Try multi-equation cases::
     809
     810        sage: x, y, z = var("x y z")
     811        sage: solve_mod([2*x^2 + x*y, -x*y+2*y^2+x-2*y, -2*x^2+2*x*y-y^2-x-y], 12)
     812        [(0, 0), (4, 4), (0, 3), (4, 7)]
     813        sage: eqs = [-y^2+z^2, -x^2+y^2-3*z^2-z-1, -y*z-z^2-x-y+2, -x^2-12*z^2-y+z]
     814        sage: solve_mod(eqs, 11)
     815        [(8, 5, 6)]
     816
     817    Confirm that modulus 1 now behaves as it should::
     818   
     819        sage: x, y = var("x y")
     820        sage: solve_mod([x==1], 1)
     821        [(0,)]
     822        sage: solve_mod([2*x^2+x*y, -x*y+2*y^2+x-2*y, -2*x^2+2*x*y-y^2-x-y], 1)
     823        [(0, 0)]
     824   
     825
    788826    """
    789     from sage.rings.all import Integer, Integers, PolynomialRing, factor, crt_basis
     827    from sage.rings.all import Integer, Integers, crt_basis
     828    from sage.symbolic.expression import is_Expression
    790829    from sage.misc.all import cartesian_product_iterator
    791830    from sage.modules.all import vector
    792831    from sage.matrix.all import matrix
    793832
    794833    if not isinstance(eqns, (list, tuple)):
    795834        eqns = [eqns]
     835    eqns = [eq if is_Expression(eq) else (eq.lhs()-eq.rhs()) for eq in eqns]
    796836    modulus = Integer(modulus)
    797837    if modulus < 1:
    798838        raise ValueError, "the modulus must be a positive integer"
    799839    vars = list(set(sum([list(e.variables()) for e in eqns], [])))
    800     vars.sort(cmp = lambda x,y: cmp(repr(x), repr(y)))
    801     n = len(vars)
     840    vars.sort(key=repr)
     841   
     842    if modulus == 1: # degenerate case
     843        ans = [tuple(Integers(1)(0) for v in vars)]
     844        return ans
    802845
    803     factors = [p**i for p,i in factor(modulus)]
    804     crt_basis = vector(Integers(modulus), crt_basis(factors))
    805     solutions = [solve_mod_enumerate(eqns, p) for p in factors]
     846    factors = modulus.factor()
     847    crt_basis = vector(Integers(modulus), crt_basis([p**i for p,i in factors]))
     848    solutions = []
     849   
     850    has_solution = True
     851    for p,i in factors:
     852        solution =_solve_mod_prime_power(eqns, p, i, vars)
     853        if len(solution) > 0:
     854            solutions.append(solution)
     855        else:
     856            has_solution = False
     857            break
     858           
    806859
    807860    ans = []
    808     for solution in cartesian_product_iterator(solutions):
    809         solution_mat = matrix(Integers(modulus), solution)
    810         ans.append(tuple(c.dot_product(crt_basis) for c in solution_mat.columns()))
     861    if has_solution:
     862        for solution in cartesian_product_iterator(solutions):
     863            solution_mat = matrix(Integers(modulus), solution)
     864            ans.append(tuple(c.dot_product(crt_basis) for c in solution_mat.columns()))
    811865
    812866    # if solution_dict == True:
    813867    # Relaxed form suggested by Mike Hansen (#8553):
     
    817871    else:
    818872        return ans
    819873
    820 def solve_mod_enumerate(eqns, modulus):
     874def _solve_mod_prime_power(eqns, p, m, vars):
    821875    r"""
    822     Return all solutions to an equation or list of equations modulo the
    823     given integer modulus. Each equation must involve only polynomials
     876    Internal help function for solve_mod, does little checking since it expects
     877    solve_mod to do that
     878
     879    Return all solutions to an equation or list of equations modulo p^m.
     880    Each equation must involve only polynomials
    824881    in 1 or many variables.
    825882
    826883    The solutions are returned as `n`-tuples, where `n`
    827     is the number of variables appearing anywhere in the given
    828     equations. The variables are in alphabetical order.
    829    
     884    is the number of variables in vars.
     885       
    830886    INPUT:
    831887   
    832888   
    833889    -  ``eqns`` - equation or list of equations
    834890   
    835     -  ``modulus`` - an integer
     891    -  ``p`` - a prime
     892
     893    -  ``i`` - an integer > 0
     894
     895    -  ``vars`` - a list of variables to solve for
    836896   
    837897   
    838898    EXAMPLES::
     
    860920       
    861921    .. warning::
    862922   
    863        Currently this naively enumerates all possible solutions.  The
    864        interface is good, but the algorithm is horrible if the modulus
    865        is at all large! Sage *does* have the ability to do something
    866        much faster in certain cases at least by using the Chinese
    867        Remainder Theorem, Groebner basis, linear algebra techniques,
    868        etc. But for a lot of toy problems this function as is might be
    869        useful. At the very least, it establishes an interface.
     923       Currently this constructs possible solutions by building up
     924       from the smallest prime factor of the modulus.  The interface
     925       is good, but the algorithm is horrible if the modulus isn't the
     926       product of many small primes! Sage *does* have the ability to
     927       do something much faster in certain cases at least by using the
     928       Chinese Remainder Theorem, Groebner basis, linear algebra
     929       techniques, etc. But for a lot of toy problems this function as
     930       is might be useful. At the very least, it establishes an
     931       interface.
     932
     933    TESTS:
     934   
     935    Confirm we can reproduce the first few terms of OEIS A187719::
     936   
     937        sage: from sage.symbolic.relation import _solve_mod_prime_power
     938        sage: [sorted(_solve_mod_prime_power([x^2==41], 10, i, [x]))[0][0] for i in [1..13]]
     939        [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429,
     940        13241296179, 19473547571, 2263241296179]
     941
    870942    """
    871     from sage.rings.all import Integer, Integers, PolynomialRing
    872     from sage.symbolic.expression import is_Expression
     943    from sage.rings.all import Integers, PolynomialRing
     944    from sage.modules.all import vector
    873945    from sage.misc.all import cartesian_product_iterator
    874946   
    875     if not isinstance(eqns, (list, tuple)):
    876         eqns = [eqns]
    877     modulus = Integer(modulus)
    878     if modulus < 1:
    879         raise ValueError, "the modulus must be a positive integer"
    880     vars = list(set(sum([list(e.variables()) for e in eqns], [])))
    881     vars.sort(cmp = lambda x,y: cmp(repr(x), repr(y)))
    882     n = len(vars)
    883     R = Integers(modulus)
    884     S = PolynomialRing(R, len(vars), vars)
    885     eqns_mod = [S(eq) if is_Expression(eq) else
    886                 S(eq.lhs() - eq.rhs()) for eq in eqns]
     947    mrunning = 1
    887948    ans = []
    888     for t in cartesian_product_iterator([R]*len(vars)):
    889         is_soln = True
    890         for e in eqns_mod:
    891             if e(t) != 0:
    892                 is_soln = False
    893                 break
    894         if is_soln:
    895             ans.append(t)
     949    for mi in xrange(m):
     950        mrunning *= p
     951        R = Integers(mrunning)
     952        S = PolynomialRing(R, len(vars), vars)
     953        eqns_mod = [S(eq) for eq in eqns]
     954        if mi == 0:
     955            possibles = cartesian_product_iterator([xrange(len(R)) for _ in xrange(len(vars))])
     956        else:
     957            shifts = cartesian_product_iterator([xrange(p) for _ in xrange(len(vars))])
     958            pairs = cartesian_product_iterator([shifts, ans])
     959            possibles = (tuple(vector(t)+vector(shift)*(mrunning//p)) for shift, t in pairs)
     960        ans = list(t for t in possibles if all(e(*t) == 0 for e in eqns_mod))
     961        if not ans: return ans
    896962
    897963    return ans
    898964