| 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 | |
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 | |
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 | |
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 | |
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 |
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 |