Improve performance and stability for Polynomial root-finding

Improve performance of polynomial root-finding by solving linear polynomials directly,
instead of handing off to the Laguerre solver.
For nonlinear polynomials, improve the stability of root-finding for polynomials with very large or very small coefficients
by first conditioning the problem, normalizing the constant coefficient to 1.

contrib/src/main/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/Polynomial.java

@@ -256,19 +256,38 @@ public Expiring<Polynomial> max(Polynomial other) {
    * Finds all roots of this function in the future
    */
   private Stream<Duration> findFutureRoots() {
+    // TODO: In some sense, isn't having an infinite coefficient the same as a vertical line,
+    //   hence the same as having a root at x = 0?
+    //   Unless the value itself is non-finite, that is...
     // If this polynomial can never have a root, fail immediately
     if (this.isNonFinite() || this.isConstant()) {
       return Stream.empty();
     }
 
-    // Defining epsilon keeps the Laguerre solver fast and stable for poorly-behaved polynomials.
-    final double epsilon = 2 * Arrays.stream(coefficients).map(Math::ulp).max().orElseThrow();
+    if (coefficients[0] == 0.0) {
+      return Stream.of(ZERO);
+    }
+
+    // If the polynomial is linear, solve it analytically for performance
+    if (this.degree() <= 1) {
+      double t = -getCoefficient(0) / getCoefficient(1);
+      if (t >= -ABSOLUTE_ACCURACY_FOR_DURATIONS / 2 && t <= MAX_SECONDS_FOR_DURATION) {
+        return Stream.of(Duration.roundNearest(t, SECOND));
+      } else {
+        return Stream.empty();
+      }
+    }
+
+    // Condition the problem by dividing through by the first coefficient:
+    double[] conditionedCoefficients = Arrays.stream(coefficients).map(c -> c / coefficients[0]).toArray();
+    // Defining epsilon keeps the Laguerre solver faster and more stable for poorly-behaved polynomials.
+    final double epsilon = 2 * Arrays.stream(conditionedCoefficients).map(Math::ulp).max().orElseThrow();
     final Complex[] solutions = new LaguerreSolver(0, ABSOLUTE_ACCURACY_FOR_DURATIONS, epsilon)
-            .solveAllComplex(coefficients, 0);
+            .solveAllComplex(conditionedCoefficients, 0);
     return Arrays.stream(solutions)
                  .filter(solution -> Math.abs(solution.getImaginary()) < epsilon)
                  .map(Complex::getReal)
-                 .filter(t -> t >= 0 && t <= MAX_SECONDS_FOR_DURATION)
+                 .filter(t -> t >= -ABSOLUTE_ACCURACY_FOR_DURATIONS / 2 && t <= MAX_SECONDS_FOR_DURATION)
                  .sorted()
                  .map(t -> Duration.roundNearest(t, SECOND));
   }

contrib/src/test/java/gov/nasa/jpl/aerie/contrib/streamline/modeling/polynomial/ComparisonsTest.java

@@ -299,6 +299,80 @@ void comparing_converging_nonlinear_terms_with_fine_precision() {
     check_extrema(false, true);
   }
 
+  // Unrepresentable convergence:
+  // These tests reflect polynomials that in theory converge, but do so in timespans
+  // that are too large to represent. Thus, they should be treated as non-converging.
+
+  @Test
+  void comparing_linear_terms_with_convergence_unrepresentable_by_double() {
+    setup(() -> {
+      set(p, polynomial(Double.MAX_VALUE));
+      set(q, polynomial(0, 0.1));
+    });
+
+    check_comparison(p_lt_q, false, false);
+    check_comparison(p_lte_q, false, false);
+    check_comparison(p_gt_q, true, false);
+    check_comparison(p_gte_q, true, false);
+    check_extrema(true, false);
+  }
+
+  @Test
+  void comparing_linear_terms_with_convergence_unrepresentable_by_duration() {
+    setup(() -> {
+      set(p, polynomial(Duration.MAX_VALUE.ratioOver(SECOND)));
+      set(q, polynomial(0, 0.1));
+    });
+
+    check_comparison(p_lt_q, false, false);
+    check_comparison(p_lte_q, false, false);
+    check_comparison(p_gt_q, true, false);
+    check_comparison(p_gte_q, true, false);
+    check_extrema(true, false);
+  }
+
+  @Test
+  void comparing_nonlinear_terms_with_convergence_unrepresentable_by_double() {
+    setup(() -> {
+      set(p, polynomial(Double.MAX_VALUE));
+      set(q, polynomial(0, 0, 0.1));
+    });
+
+    check_comparison(p_lt_q, false, false);
+    check_comparison(p_lte_q, false, false);
+    check_comparison(p_gt_q, true, false);
+    check_comparison(p_gte_q, true, false);
+    check_extrema(true, false);
+  }
+
+  @Test
+  void comparing_nonlinear_terms_with_convergence_unrepresentable_by_duration() {
+    setup(() -> {
+      set(p, polynomial(Duration.MAX_VALUE.ratioOver(SECOND) * Duration.MAX_VALUE.ratioOver(SECOND)));
+      set(q, polynomial(0, 0, 0.1));
+    });
+
+    check_comparison(p_lt_q, false, false);
+    check_comparison(p_lte_q, false, false);
+    check_comparison(p_gt_q, true, false);
+    check_comparison(p_gte_q, true, false);
+    check_extrema(true, false);
+  }
+
+  @Test
+  void comparing_pathological_nonlinear_terms_with_convergence_unrepresentable_by_duration() {
+    setup(() -> {
+      set(p, polynomial(Duration.MAX_VALUE.ratioOver(SECOND) * Duration.MAX_VALUE.ratioOver(SECOND)));
+      set(q, polynomial(0, Duration.MIN_VALUE.ratioOver(SECOND), 1.0 + Math.ulp(1.0)));
+    });
+
+    check_comparison(p_lt_q, false, false);
+    check_comparison(p_lte_q, false, false);
+    check_comparison(p_gt_q, true, false);
+    check_comparison(p_gte_q, true, false);
+    check_extrema(true, false);
+  }
+
   private void check_comparison(Resource<Discrete<Boolean>> result, boolean expectedValue, boolean expectCrossover) {
     reset();
     var resultDynamics = result.getDynamics().getOrThrow();