diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/ToggleFeature.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/ToggleFeature.java
index def8e3fa24..38e2259b89 100644
--- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/ToggleFeature.java
+++ b/symja_android_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 true, enable solvers from package io.github.mangara.diophantine to
+ * find some solutions in {@link S#FindInstance}
+ */
+ public static boolean SOLVE_DIOPHANTINE = true;
+
public static boolean SHOW_STEPS = true;
}
diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/NumberTheory.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/NumberTheory.java
index 45896f6c55..a6a55c01b5 100644
--- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/NumberTheory.java
+++ b/symja_android_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 jas = new JASConvert(
+ varSet.getArrayList(), edu.jas.arith.BigInteger.ZERO);
+ GenPolynomial 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 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 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 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 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/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/convert/JASConvert.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/convert/JASConvert.java
index 78d46b7ac7..8cdd503ced 100644
--- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/convert/JASConvert.java
+++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/convert/JASConvert.java
@@ -356,7 +356,14 @@ private GenPolynomial expr2Poly(final IExpr exprPoly, boolean numeric2Rationa
// "JASConvert:expr2Poly - invalid exponent: " + ast.arg2().toString());
}
try {
- return fPolyFactory.univariate(base.getSymbolName(), exponent);
+ GenPolynomial 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/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Solve.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Solve.java
index 71f15da7a7..da0a8c9b67 100644
--- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Solve.java
+++ b/symja_android_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/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/FindInstanceTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/FindInstanceTest.java
index f6ec16a896..95f713eba8 100644
--- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/FindInstanceTest.java
+++ b/symja_android_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