WIP #1056 add https://github.com/Mangara/diophantine solvers

symja_android_library/matheclipse-external/src/main/java/io/github/mangara/diophantine/LinearSolver.java

@@ -0,0 +1,177 @@
+/*
+ * Copyright 2022 Sander Verdonschot <sander.verdonschot at gmail.com>.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
+ * in compliance with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software distributed under the License
+ * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
+ * or implied. See the License for the specific language governing permissions and limitations under
+ * the License.
+ */
+package io.github.mangara.diophantine;
+
+import io.github.mangara.diophantine.iterators.IntegerIterator;
+import io.github.mangara.diophantine.iterators.MappingIterator;
+import io.github.mangara.diophantine.iterators.XYIterator;
+import java.math.BigInteger;
+import java.util.Collections;
+import java.util.Iterator;
+import java.util.function.Function;
+
+/**
+ * Solves linear Diophantine equations in two variables.
+ * <p>
+ * The method is based on https://www.alpertron.com.ar/METHODS.HTM#Linear
+ */
+public class LinearSolver {
+
+  /**
+   * Solves the linear Diophantine equation d x + e y + f = 0.
+   *
+   * @param d
+   * @param e
+   * @param f
+   * @return an iterator over all integer solutions (x, y)
+   */
+  public static Iterator<XYPair> solve(long d, long e, long f) {
+    return solve(BigInteger.valueOf(d), BigInteger.valueOf(e), BigInteger.valueOf(f));
+  }
+
+  /**
+   * Solves the linear Diophantine equation d x + e y + f = 0.
+   *
+   * @param d
+   * @param e
+   * @param f
+   * @return an iterator over all integer solutions (x, y)
+   */
+  public static Iterator<XYPair> solve(BigInteger d, BigInteger e, BigInteger f) {
+    if (d.signum() == 0 && e.signum() == 0) {
+      return solveTrivial(f);
+    } else if (d.signum() == 0) {
+      return solveSingle(e, f, true);
+    } else if (e.signum() == 0) {
+      return solveSingle(d, f, false);
+    } else {
+      return solveGeneral(d, e, f);
+    }
+  }
+
+  private static Iterator<XYPair> solveTrivial(BigInteger f) {
+    if (f.signum() != 0) {
+      return Collections.emptyIterator();
+    }
+
+    return new XYIterator();
+  }
+
+  /**
+   * Solves the linear Diophantine equations g x + f = 0 or g y + f = 0.
+   * <p>
+   * If arbitraryX is true, the equation solved is g y + f = 0 and x ranges over all integers.
+   * Otherwise, it is g x + f = 0, with y ranging over all integers.
+   *
+   * @param g
+   * @param f
+   * @param arbitraryX
+   * @return an iterator over all integer solutions (x, y)
+   */
+  private static Iterator<XYPair> solveSingle(BigInteger g, BigInteger f, boolean arbitraryX) {
+    if (f.mod(g.abs()).signum() != 0) {
+      return Collections.emptyIterator();
+    }
+
+    BigInteger val = f.negate().divide(g);
+
+    Function<BigInteger, XYPair> map;
+    if (arbitraryX) {
+      map = k -> new XYPair(k, val);
+    } else {
+      map = k -> new XYPair(val, k);
+    }
+
+    return new MappingIterator<>(new IntegerIterator(), map);
+  }
+
+  // Pre: d != 0 && e != 0
+  private static Iterator<XYPair> solveGeneral(BigInteger d, BigInteger e, BigInteger f) {
+    Eq reduced = reduce(d, e, f);
+
+    if (reduced == null) {
+      return Collections.emptyIterator();
+    } else {
+      return solveReduced(reduced);
+    }
+  }
+
+  private static Eq reduce(BigInteger d, BigInteger e, BigInteger f) {
+    BigInteger gcd = d.gcd(e).abs();
+
+    if (f.mod(gcd).signum() != 0) {
+      // No solutions, as d x + e y will always be a multiple of gcd for integer x and y
+      return null;
+    }
+
+    return new Eq(d.divide(gcd), e.divide(gcd), f.divide(gcd));
+  }
+
+  // Pre: d != 0 && e != 0 && d and e are co-prime
+  private static Iterator<XYPair> solveReduced(Eq eq) {
+    XYPair solution = findAnySolution(eq);
+
+    final BigInteger dx = eq.e;
+    final BigInteger dy = eq.d.negate();
+
+    return new MappingIterator<>(new IntegerIterator(),
+        (k) -> new XYPair(solution.x.add(dx.multiply(k)), solution.y.add(dy.multiply(k))));
+  }
+
+  private static XYPair findAnySolution(Eq eq) {
+    // Run the extended Euclidean algorithm (
+    // https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm )
+    BigInteger prevR = eq.d;
+    BigInteger r = eq.e;
+    BigInteger prevS = BigInteger.ONE;
+    BigInteger s = BigInteger.ZERO;
+    BigInteger prevT = BigInteger.ZERO;
+    BigInteger t = BigInteger.ONE;
+
+    while (r.signum() != 0) {
+      BigInteger quotient = prevR.divide(r);
+      BigInteger tempR = r;
+      BigInteger tempS = s;
+      BigInteger tempT = t;
+
+      r = prevR.subtract(quotient.multiply(r));
+      s = prevS.subtract(quotient.multiply(s));
+      t = prevT.subtract(quotient.multiply(t));
+
+      prevR = tempR;
+      prevS = tempS;
+      prevT = tempT;
+    }
+
+    // Results are in prevR (which is the gcd = +/- 1), prevS, and prevT
+    // Thus, d * prevS + e * prevT = +/- 1
+    // We want d * x + e * y = -f, so we need to multiply by f or -f, depending on the sign of prevR
+    BigInteger factor = eq.f.negate().multiply(prevR);
+    BigInteger x = factor.multiply(prevS); // -/+ f * prevS
+    BigInteger y = factor.multiply(prevT); // -/+ f * prevT
+
+    return new XYPair(x, y);
+  }
+
+  private static class Eq {
+
+    BigInteger d, e, f;
+
+    Eq(BigInteger d, BigInteger e, BigInteger f) {
+      this.d = d;
+      this.e = e;
+      this.f = f;
+    }
+  }
+}

symja_android_library/matheclipse-external/src/main/java/io/github/mangara/diophantine/QuadraticSolver.java

@@ -0,0 +1,92 @@
+/*
+ * Copyright 2022 Sander Verdonschot <sander.verdonschot at gmail.com>.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
+ * in compliance with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software distributed under the License
+ * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
+ * or implied. See the License for the specific language governing permissions and limitations under
+ * the License.
+ */
+package io.github.mangara.diophantine;
+
+import java.math.BigInteger;
+import java.util.Collections;
+import java.util.Iterator;
+import io.github.mangara.diophantine.quadratic.EllipticalSolver;
+import io.github.mangara.diophantine.quadratic.HyperbolicSolver;
+import io.github.mangara.diophantine.quadratic.ParabolicSolver;
+import io.github.mangara.diophantine.quadratic.SquareDiscriminantSolver;
+
+/**
+ * Solves quadratic Diophantine equations in two variables.
+ * <p>
+ * The method is based on K. R. Matthews, "Solving the Diophantine equation \(ax^2 + bxy + cy^2 + dx
+ * + ey + f = 0\)" http://www.numbertheory.org/PDFS/general_quadratic_solution.pdf
+ */
+public class QuadraticSolver {
+
+  /**
+   * Solves the quadratic Diophantine equation a x^2 + b xy + c y^2 + d x + e y + f = 0.
+   *
+   * @param a
+   * @param b
+   * @param c
+   * @param d
+   * @param e
+   * @param f
+   * @return an iterator over all integer solutions (x, y)
+   */
+  public static Iterator<XYPair> solve(long a, long b, long c, long d, long e, long f) {
+    return solve(BigInteger.valueOf(a), BigInteger.valueOf(b), BigInteger.valueOf(c),
+        BigInteger.valueOf(d), BigInteger.valueOf(e), BigInteger.valueOf(f));
+  }
+
+  /**
+   * Solves the quadratic Diophantine equation a x^2 + b xy + c y^2 + d x + e y + f = 0.
+   *
+   * @param a
+   * @param b
+   * @param c
+   * @param d
+   * @param e
+   * @param f
+   * @return an iterator over all integer solutions (x, y)
+   */
+  public static Iterator<XYPair> solve(BigInteger a, BigInteger b, BigInteger c, BigInteger d,
+      BigInteger e, BigInteger f) {
+    if (a.signum() == 0 && b.signum() == 0 && c.signum() == 0) {
+      return LinearSolver.solve(d, e, f);
+    }
+
+    BigInteger D = Utils.discriminant(a, b, c);
+
+    if (D.signum() == 0) {
+      return ParabolicSolver.solve(a, b, c, d, e, f);
+    } else if (Utils.isSquare(D)) {
+      return SquareDiscriminantSolver.solve(a, b, c, d, e, f);
+    } else if (Utils.legendreConstant(a, b, c, d, e, f, D).signum() == 0) {
+      return solveTrivialCase(a, b, c, d, e, f, D);
+    } else if (D.signum() < 0) {
+      return EllipticalSolver.solve(a, b, c, d, e, f);
+    } else {
+      return HyperbolicSolver.solve(a, b, c, d, e, f);
+    }
+  }
+
+  // Pre: D = b^2 - 4ac != 0 && D not a perfect square && Legendre's k = 0
+  private static Iterator<XYPair> solveTrivialCase(BigInteger a, BigInteger b, BigInteger c,
+      BigInteger d, BigInteger e, BigInteger f, BigInteger D) {
+    BigInteger alpha = Utils.legendreAlpha(b, c, d, e);
+    BigInteger beta = Utils.legendreBeta(a, b, d, e);
+
+    if (alpha.mod(D).signum() == 0 && beta.mod(D).signum() == 0) {
+      return Collections.singletonList(new XYPair(alpha.divide(D), beta.divide(D))).iterator();
+    } else {
+      return Collections.emptyIterator();
+    }
+  }
+}