WIP #1056 call diophantine solvers from Function FindInstance

- set ToggleFeature.SOLVE_DIOPHANTINE = true to test the experimental implementation
axkr · Aug 24, 2024 · 0a9e4c1 · 0a9e4c1
1 parent 724c550
commit 0a9e4c1
Show file tree

Hide file tree

Showing 5 changed files with 187 additions and 5 deletions.
diff --git a/...roid_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/ToggleFeature.java b/...roid_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/ToggleFeature.java
@@ -58,5 +58,11 @@ public class ToggleFeature {
    */
   public static boolean SERIES_DENOMINATOR = false;
 
+  /**
+   * If <code>true</code>, enable solvers from package <code>io.github.mangara.diophantine</code> to
+   * find some solutions in {@link S#FindInstance}
+   */
+  public static boolean SOLVE_DIOPHANTINE = true;
+
   public static boolean SHOW_STEPS = true;
 }
diff --git a/...oid_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/NumberTheory.java b/...oid_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/NumberTheory.java
@@ -16,6 +16,7 @@
 import java.math.BigInteger;
 import java.util.ArrayList;
 import java.util.HashMap;
+import java.util.Iterator;
 import java.util.List;
 import java.util.Map;
 import java.util.SortedMap;
@@ -79,6 +80,7 @@
 import org.matheclipse.core.interfaces.ISymbol;
 import org.matheclipse.core.numbertheory.GaussianInteger;
 import org.matheclipse.core.numbertheory.Primality;
+import org.matheclipse.core.polynomials.QuarticSolver;
 import org.matheclipse.core.sympy.series.Sequences;
 import org.matheclipse.core.visit.VisitorExpr;
 import com.google.common.math.BigIntegerMath;
@@ -88,8 +90,14 @@
 import edu.jas.arith.ModInteger;
 import edu.jas.arith.ModIntegerRing;
 import edu.jas.poly.GenPolynomial;
+import edu.jas.poly.Monomial;
 import edu.jas.ufd.FactorAbstract;
 import edu.jas.ufd.FactorFactory;
+import io.github.mangara.diophantine.QuadraticSolver;
+import io.github.mangara.diophantine.Utils;
+import io.github.mangara.diophantine.XYPair;
+import io.github.mangara.diophantine.quadratic.ParabolicSolver;
+import io.github.mangara.diophantine.quadratic.PellsSolver;
 
 public final class NumberTheory {
 
@@ -6616,4 +6624,123 @@ private static IAST quadraticIrrationalPlus(IAST plusAST, IASTMutable resultList
     }
     return F.NIL;
   }
+
+  public static IAST diophantinePolynomial(final IExpr expr, IAST varList,
+      int maximumNumberOfResults) {
+    VariablesSet varSet = new VariablesSet(varList);
+    IASTMutable result = F.NIL;
+    try {
+      // try to generate a common expression polynomial
+      JASConvert<edu.jas.arith.BigInteger> jas = new JASConvert<edu.jas.arith.BigInteger>(
+          varSet.getArrayList(), edu.jas.arith.BigInteger.ZERO);
+      GenPolynomial<edu.jas.arith.BigInteger> ePoly = jas.expr2JAS(expr, false);
+      result = diophantinePolynomial(ePoly, varList, maximumNumberOfResults);
+      result = QuarticSolver.sortASTArguments(result);
+      return result;
+    } catch (JASConversionException e2) {
+      e2.printStackTrace();
+    }
+    return result;
+  }
+
+  private static IASTAppendable diophantinePolynomial(
+      GenPolynomial<edu.jas.arith.BigInteger> polynomial, IAST varList,
+      int maximumNumberOfResults) {
+    long varDegree = polynomial.degree(0);
+
+    if (polynomial.isConstant()) {
+      return F.ListAlloc(1);
+    }
+    // a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
+    BigInteger a = BigInteger.ZERO;
+    BigInteger b = BigInteger.ZERO;
+    BigInteger c = BigInteger.ZERO;
+    BigInteger d = BigInteger.ZERO;
+    BigInteger e = BigInteger.ZERO;
+    BigInteger f = BigInteger.ZERO;
+    try {
+      if (varDegree <= 2) {
+        if (varList.argSize() == 1) {
+          // x is only variable => b=0;c=0;e=0;
+          for (Monomial<edu.jas.arith.BigInteger> monomial : polynomial) {
+            edu.jas.arith.BigInteger coeff = monomial.coefficient();
+            BigInteger zz = coeff.val;
+            long xExp = monomial.exponent().getVal(0);
+            if (xExp == 2) {
+              a = zz;
+            } else if (xExp == 1) {
+              d = zz;
+            } else if (xExp == 2) {
+              f = zz;
+            } else {
+              throw new ArithmeticException(
+                  "diophantinePolynomial::Unexpected exponent value: " + xExp);
+            }
+          }
+        } else if (varList.argSize() == 2) {
+          // x and y are both variables
+          // a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
+          for (Monomial<edu.jas.arith.BigInteger> monomial : polynomial) {
+            edu.jas.arith.BigInteger coeff = monomial.coefficient();
+            BigInteger zz = coeff.val;
+            int xi = monomial.exponent().varIndex(0);
+            int yi = monomial.exponent().varIndex(1);
+            long xExp = monomial.exponent().getVal(xi);
+            long yExp = monomial.exponent().getVal(yi);
+            if (xExp == 2 && yExp == 0) {
+              a = zz;
+            } else if (xExp == 1 && yExp == 1) {
+              b = zz;
+            } else if (xExp == 0 && yExp == 2) {
+              c = zz;
+            } else if (xExp == 1 && yExp == 0) {
+              d = zz;
+            } else if (xExp == 0 && yExp == 1) {
+              e = zz;
+            } else if (xExp == 0 && yExp == 0) {
+              f = zz;
+            } else {
+              throw new ArithmeticException(
+                  "diophantinePolynomial::Unexpected exponent value: " + xExp);
+            }
+          }
+        }
+        Iterator<XYPair> diophantineSolver = null;
+        IASTAppendable result = F.ListAlloc();
+        if (maximumNumberOfResults == 1 && c.signum() < 0 && f.equals(BigInteger.valueOf(-4))
+            && b.signum() == 0 && d.signum() == 0 && e.signum() == 0) {
+          BigInteger cNegate = c.negate();
+          if (!Utils.isSquare(cNegate)) {
+            // use Pell's equation if -c is not a perfect square
+            XYPair xyPair = PellsSolver.leastPositivePellsFourSolution(cNegate);
+            result.append(F.List(F.Rule(varList.arg1(), F.ZZ(xyPair.x)),
+                F.Rule(varList.arg2(), F.ZZ(xyPair.y))));
+            return result;
+          }
+        }
+        if (a.signum() != 0 || b.signum() != 0 || c.signum() != 0) {
+          // D := b^2 - 4ac && D == 0
+          BigInteger D = b.multiply(b).subtract(a.multiply(c).multiply(BigInteger.valueOf(4)));
+          if (D.signum() == 0) {
+            diophantineSolver = ParabolicSolver.solve(a, b, c, d, e, f);
+          }
+        }
+        if (diophantineSolver == null) {
+          diophantineSolver = QuadraticSolver.solve(a, b, c, d, e, f);
+        }
+
+        int n = 0;
+        while (diophantineSolver.hasNext() && n++ < maximumNumberOfResults) {
+          XYPair xyPair = diophantineSolver.next();
+          result.append(F.List(F.Rule(varList.arg1(), F.ZZ(xyPair.x)),
+              F.Rule(varList.arg2(), F.ZZ(xyPair.y))));
+        }
+        return result;
+      }
+    } catch (ArithmeticException aex) {
+      //
+    }
+
+    return F.NIL;
+  }
 }
diff --git a/...droid_library/matheclipse-core/src/main/java/org/matheclipse/core/convert/JASConvert.java b/...droid_library/matheclipse-core/src/main/java/org/matheclipse/core/convert/JASConvert.java
@@ -356,7 +356,14 @@ private GenPolynomial<C> expr2Poly(final IExpr exprPoly, boolean numeric2Rationa
             // "JASConvert:expr2Poly - invalid exponent: " + ast.arg2().toString());
           }
           try {
-            return fPolyFactory.univariate(base.getSymbolName(), exponent);
+            GenPolynomial<C> v = fPolyFactory.univariate(base.getSymbolName(), 1L);
+            return v.power(exponent);
+            // int indexOf = fVariables.indexOf(base);
+            // if (indexOf >= 0) {
+            // ExpVectorLong v = new ExpVectorLong(fVariables.size(), indexOf, exponent);
+            // return fPolyFactory.valueOf(fRingFactory.getONE(), v);
+            // }
+            // return fPolyFactory.univariate(base.getSymbolName(), exponent);
           } catch (IllegalArgumentException iae) {
             // fall through
           }

diff --git a/..._library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Solve.java b/..._library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Solve.java
@@ -9,9 +9,11 @@
 import org.apache.logging.log4j.Logger;
 import org.hipparchus.linear.FieldMatrix;
 import org.matheclipse.core.basic.Config;
+import org.matheclipse.core.basic.ToggleFeature;
 import org.matheclipse.core.builtin.Algebra;
 import org.matheclipse.core.builtin.BooleanFunctions;
 import org.matheclipse.core.builtin.LinearAlgebra;
+import org.matheclipse.core.builtin.NumberTheory;
 import org.matheclipse.core.builtin.PolynomialFunctions;
 import org.matheclipse.core.builtin.RootsFunctions;
 import org.matheclipse.core.convert.ChocoConvert;
@@ -1289,6 +1291,19 @@ public static IExpr solveIntegers(final IAST ast, IAST equationVariables,
         return F.NIL;
       }
       try {
+        if (ToggleFeature.SOLVE_DIOPHANTINE) {
+          if (equationsAndInequations.argSize() == 1) {
+            IExpr eq1 = equationsAndInequations.arg1();
+            if (eq1.isEqual() && eq1.second().isZero()) {
+              IAST diophantineResult = NumberTheory.diophantinePolynomial(eq1.first(),
+                  equationVariables, maximumNumberOfResults);
+              if (diophantineResult.isPresent()) {
+                return diophantineResult;
+              }
+            }
+          }
+        }
+
         if (equationsAndInequations.isFreeAST(x -> chocoSolver(x))) {
           // choco-solver doesn't handle Power() expressions very well at the moment!
           try {

diff --git a/..._library/matheclipse-core/src/test/java/org/matheclipse/core/system/FindInstanceTest.java b/..._library/matheclipse-core/src/test/java/org/matheclipse/core/system/FindInstanceTest.java
@@ -1,16 +1,43 @@
 package org.matheclipse.core.system;
 
 import org.junit.Test;
+import org.matheclipse.core.basic.ToggleFeature;
 
 /** Tests for FindInstance function */
 public class FindInstanceTest extends ExprEvaluatorTestCase {
 
   @Test
-  public void testDiophantine() {
+  public void testDiophantine001() {
+    if (ToggleFeature.SOLVE_DIOPHANTINE) {
+      // ParabolicSolver
+      check("FindInstance(9*x^2 - 12*x*y + 4*y^2 + 3*x + 2*y == 12,{x,y},Integers, 3)", //
+          "{{x->0,y->-2},{x->1,y->0},{x->2,y->3}}");
+      // QuadraticSolver
+      check("FindInstance(3*x^2 - 8*x*y + 7*y^2 - 4*x + 2*y - 109 == 0,{x,y},Integers, 3)", //
+          "{{x->2,y->-3},{x->2,y->5},{x->14,y->9}}");
+      // Pell 4
+      check("FindInstance(x^2-29986*y^2-4 ==0,{x,y},Integers,1)", //
+          "{{x->135915148103491619905402044543098,y->784889635731418443294120995460}}");
+
+      check("FindInstance(x^2-29986*y^2-4 ==0,{x,y},Integers,3)", //
+          "{{x->-135915148103491619905402044543098,y->784889635731418443294120995460},{x->-\n" //
+              + "2,y->0},{x->135915148103491619905402044543098,y->-784889635731418443294120995460}}");
+    }
+  }
+
+  @Test
+  public void testDiophantine002() {
     // TODO return condition with extra variable C1
-    check("FindInstance(13*x+51*y==0, {x,y}, Integers, 6)", //
-        "{{x->-153,y->39},{x->-102,y->26},{x->-51,y->13},{x->0,y->0},{x->51,y->-13},{x->\n"
-            + "102,y->-26}}");
+    if (ToggleFeature.SOLVE_DIOPHANTINE) {
+      check("FindInstance(13*x+51*y==0, {x,y}, Integers, 6)", //
+          "{{x->-102,y->26},{x->-51,y->13},{x->0,y->0},{x->51,y->-13},{x->102,y->-26},{x->\n" //
+              + "153,y->-39}}");
+    } else {
+      // choco-solver solutions:
+      check("FindInstance(13*x+51*y==0, {x,y}, Integers, 6)", //
+          "{{x->-153,y->39},{x->-102,y->26},{x->-51,y->13},{x->0,y->0},{x->51,y->-13},{x->\n" //
+              + "102,y->-26}}");
+    }
   }
 
   @Test