# Ticket #11036: trac_11036_improve_solve_mod_perf.patch

File trac_11036_improve_solve_mod_perf.patch, 9.4 KB (added by dsm, 8 years ago)
• ## sage/symbolic/relation.py

```# HG changeset patch
# User D. S. McNeil <dsm054@gmail.com>
# Date 1301110343 -28800
# Node ID 733adae06130875381160304a9845d4ff6a4df47
# Parent  361a4ad7d52c69b64ae2e658ffd0820af0d87e93
Trac 11036: improve solve_mod performance

diff -r 361a4ad7d52c -r 733adae06130 sage/symbolic/relation.py```
 a sage: solve_mod([x^2 == 1, 4*x  == 11], 15) [(14,)] Note that the solution elements belong to the relevant modular ring:: sage: sols = solve_mod([x^2 == 1, 4*x  == 11], 15) sage: sols[0][0] 14 sage: parent(sols[0][0]) Ring of integers modulo 15 Fermat's equation modulo 3 with exponent 5:: sage: var('x,y,z') sage: solve_mod([x^5 + y^5 == z^5], 3) [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)] We can solve with respect to a bigger modulus if it consists only of small prime factors:: sage: [d] = solve_mod([5*x + y == 3, 2*x - 3*y == 9], 3*5*7*11*19*23*29, solution_dict = True) sage: d[y] 8610183 We solve an simple equation modulo 2:: For cases where there are relatively few solutions and the prime factors are small, this can be efficient even if the modulus itself is large:: sage: sorted(solve_mod([x^2 == 41], 10^20)) [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,), (45461397519473547571,), (54538602480526452429,), (61445932736758703821,), (88554067263241296179,), (95461397519473547571,)] We solve a simple equation modulo 2:: sage: x,y = var('x,y') sage: solve_mod([x == y], 2) [(0, 0), (1, 1)] .. warning:: The current implementation splits the modulus into prime powers, then naively enumerates all possible solutions and finally combines the solution using the Chinese Remainder Theorem. The interface is good, but the algorithm is horrible if the modulus has some larger prime factors! Sage *does* have the ability to do something much faster in certain cases at least by using Groebner basis, linear algebra techniques, etc. But for a lot of toy problems this function as is might be useful. At least it establishes an interface. The current implementation splits the modulus into prime powers, then naively enumerates all possible solutions and finally combines the solution using the Chinese Remainder Theorem.  The interface is good, but the algorithm is horrible if the modulus has some larger prime factors! Sage *does* have the ability to do something much faster in certain cases at least by using Groebner basis, linear algebra techniques, etc. But for a lot of toy problems this function as is might be useful. At least it establishes an interface. TESTS: Make sure that we short-circuit in at least some cases:: sage: solve_mod([2*x==1], 2*next_prime(10^50)) [] Try multi-equation cases:: sage: x, y, z = var("x y z") 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) [(0, 0), (4, 4), (0, 3), (4, 7)] 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] sage: solve_mod(eqs, 11) [(8, 5, 6)] Confirm that modulus 1 now behaves as it should:: sage: x, y = var("x y") sage: solve_mod([x==1], 1) [(0,)] 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) [(0, 0)] """ from sage.rings.all import Integer, Integers, PolynomialRing, factor, crt_basis from sage.misc.all import cartesian_product_iterator if modulus < 1: raise ValueError, "the modulus must be a positive integer" vars = list(set(sum([list(e.variables()) for e in eqns], []))) vars.sort(cmp = lambda x,y: cmp(repr(x), repr(y))) n = len(vars) vars.sort(key=repr) factors = [p**i for p,i in factor(modulus)] crt_basis = vector(Integers(modulus), crt_basis(factors)) solutions = [solve_mod_enumerate(eqns, p) for p in factors] if modulus == 1: # degenerate case ans = [tuple(Integers(modulus)(0) for v in vars)] else: factors = [p**i for p,i in factor(modulus)] crt_basis = vector(Integers(modulus), crt_basis(factors)) solutions = [] for p in factors: s = solve_mod_enumerate(eqns, p) if not s: return [] solutions.append(s) ans = [] for solution in cartesian_product_iterator(solutions): solution_mat = matrix(Integers(modulus), solution) ans.append(tuple(c.dot_product(crt_basis) for c in solution_mat.columns())) ans = [] for solution in cartesian_product_iterator(solutions): solution_mat = matrix(Integers(modulus), solution) ans.append(tuple(c.dot_product(crt_basis) for c in solution_mat.columns())) if solution_dict == True: sol_dict = [dict(zip(vars, solution)) for solution in ans] sage: solve_mod([x^5 + y^5 == z^5], 3) [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)] We solve an simple equation modulo 2:: We solve a simple equation modulo 2:: sage: x,y = var('x,y') sage: solve_mod([x == y], 2) .. warning:: Currently this naively enumerates all possible solutions.  The interface is good, but the algorithm is horrible if the modulus is at all large! Sage *does* have the ability to do something much faster in certain cases at least by using the Chinese Remainder Theorem, Groebner basis, linear algebra techniques, etc. But for a lot of toy problems this function as is might be useful. At the very least, it establishes an interface. """ Currently this constructs possible solutions by building up from the smallest prime factor of the modulus.  The interface is good, but the algorithm is horrible if the modulus isn't the product of many small primes! Sage *does* have the ability to do something much faster in certain cases at least by using the Chinese Remainder Theorem, Groebner basis, linear algebra techniques, etc. But for a lot of toy problems this function as is might be useful. At the very least, it establishes an interface. TESTS: Confirm we can reproduce the first few terms of OEIS A187719:: sage: from sage.symbolic.relation import solve_mod_enumerate sage: [sorted(solve_mod_enumerate(x^2==41, 10^i))[0][0] for i in [1..13]] [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429, 13241296179, 19473547571, 2263241296179] Confirm that modulus 1 now behaves as it should:: sage: x, y = var("x y") sage: solve_mod_enumerate([2*x^2 + x*y, -x*y + 2*y^2 + x - 2*y, -2*x^2 + 2*x*y - y^2 - x - y], 1) [(0, 0)] """ from sage.rings.all import Integer, Integers, PolynomialRing from sage.rings.all import factor from sage.symbolic.expression import is_Expression from sage.misc.all import cartesian_product_iterator from sage.modules.all import vector if not isinstance(eqns, (list, tuple)): eqns = [eqns] eqns = [eq if is_Expression(eq) else (eq.lhs()-eq.rhs()) for eq in eqns] modulus = Integer(modulus) if modulus < 1: raise ValueError, "the modulus must be a positive integer" vars = list(set(sum([list(e.variables()) for e in eqns], []))) vars.sort(cmp = lambda x,y: cmp(repr(x), repr(y))) n = len(vars) R = Integers(modulus) S = PolynomialRing(R, len(vars), vars) eqns_mod = [S(eq) if is_Expression(eq) else S(eq.lhs() - eq.rhs()) for eq in eqns] vars.sort(key=repr) if modulus == 1: # degenerate case ans = [tuple(Integers(modulus)(0) for v in vars)] return ans factors = factor(modulus) mrunning = 1 ans = [] for t in cartesian_product_iterator([R]*len(vars)): is_soln = True for e in eqns_mod: if e(t) != 0: is_soln = False break if is_soln: ans.append(t) for fi, (p, m) in enumerate(factors): for mi in xrange(m): mrunning *= p R = Integers(mrunning) S = PolynomialRing(R, len(vars), vars) eqns_mod = [S(eq) for eq in eqns] if fi == 0 and mi == 0: possibles = cartesian_product_iterator([xrange(len(R)) for _ in xrange(len(vars))]) else: shifts = cartesian_product_iterator([xrange(p) for _ in xrange(len(vars))]) pairs = cartesian_product_iterator([shifts, ans]) possibles = (tuple(vector(t)+vector(shift)*(mrunning//p)) for shift, t in pairs) ans = list(t for t in possibles if all(e(t) == 0 for e in eqns_mod)) if not ans: return ans return ans def solve_ineq_univar(ineq):