Quadratic deal with a=0 and speed up

[FrogGuaGuaGua]

add an double array, _factorialLnCache,

[jkalias]

Please add some unit tests which exercise these new additions

[jkalias]

Why not write here directly Complex.NaN?

[FrogGuaGuaGua]

Because "Complex.NaN" is not available in net48 and netstandard2.0, it is introduced since net5.0, "Complex.Zero / Complex.Zero" or "Complex(double.NaN, double.NaN)" is used as a substitute for "Complex.NaN" in early versions.

[jkalias]

Ok, understood. Then I would preferably write Complex(double.NaN, double.NaN)

[jkalias]

I find this a bit confusing, I would prefer the standard textbook notation of b^2 - 4ac

[FrogGuaGuaGua]

The computational cost of Complex( , ).SquareRoot() and Complex division(return ( , c/q)) is too high and completely unnecessary. So I use Math.Sqrt() to perform the square root operation of real numbers, which has much lower computational cost. As for why it's not written in the well-known form x1,x2=(-b±sqrt(b * b-4 * a * c))/(2 * a), it is because we need to reduce the number of floating-point multiplication and division operations, especially division(which costs too much time), and temporarily store some intermediate results to avoid repeated calculations.

[jkalias]

It doesn't seem to me that performance here is a big deal, I think the discussion leans more towards premature optimization. Did you at least make any benchmarking tests in various platforms and architectures which backup this claim?

Just to make my earlier comment clear: I wouldn't reject the PR purely on this, I think unit tests are more critical, and even though performance "might" be better currently (a big IF), I think the readability is more important in such a big code base.

[FrogGuaGuaGua]

public static (Complex, Complex) FasterQuadratic(double c, double b, double a)
{
   Complex x1, x2;
   if (a == 0.0)
   {
	   if (b == 0.0)
	   {
		   x1 = new Complex(double.NaN, double.NaN);  // Complex.NaN;
		   x2 = x1;
	   }
	   else
	   {
		   x1 = new Complex(-c / b, 0.0);
		   x2 = x1;
	   }
   }
   else
   {
	   a = 1.0 / a;
	   b = -0.5 * b * a;
	   c = c * a;
	   double delta = b * b - c;
	   if (delta < 0.0)
	   {
		   double sqrtDelta = Math.Sqrt(-delta);
		   x1 = new Complex(b, sqrtDelta);
		   x2 = new Complex(b, -sqrtDelta);
	   }
	   else
	   {
		   double sqrtDelta = Math.Sqrt(delta);
		   x1 = new Complex(b + sqrtDelta, 0.0);
		   x2 = new Complex(b - sqrtDelta, 0.0);
	   }
   }
   return (x1, x2);
}

[TestCase(70000000)]
public void TestQuadraticSpeed(int N)
{
   Complex x1, x2;
   double sum1 = 0.0;
   double a, b, c;
   Stopwatch stopwatch = new Stopwatch();
   for (int i = 11; i < N; i++)
   {
	   a = -0.1 * i + 1.0;
	   b = 0.3464 * i + 5.0;
	   c = -0.3 * i + 7.0;
	   (x1, x2) = FindRoots.Quadratic(c, b, a);
	   sum1 += x1.Real + x2.Real + Math.Abs(x2.Imaginary);
   }
   Console.WriteLine($"sum1={sum1}");
   stopwatch.Stop();
   Console.WriteLine($"FindRoots.Quadratic time: {stopwatch.ElapsedMilliseconds * 0.001:F3}s");

   double sum2 = 0.0;
   stopwatch.Restart();
   for (int i = 11; i < N; i++)
   {
	   a = -0.1 * i + 1.0;
	   b = 0.3464 * i + 5.0;
	   c = -0.3 * i + 7.0;
	   (x1, x2) = FasterQuadratic(c, b, a);
	   sum2 += x1.Real + x2.Real + Math.Abs(x2.Imaginary);
   }
   Console.WriteLine($"sum2={sum2}");
   stopwatch.Stop();
   Console.WriteLine($"FasterQuadratic time: {stopwatch.ElapsedMilliseconds * 0.001:F3}s");
   Assert.AreEqual(sum1, sum2, 1e-12 * sum2);
}

[TestCase(7000000)]
public void TestFasterQuadraticCorrectness(int N)
{
   Complex x1, x2, x3, x4;
   double a, b, c, b_a, c_a;
   for (int i = 11; i < N; i++)
   {
	   a = -0.1 * i + 1.0;
	   b = 0.3464 * i + 5.0;
	   c = -0.3 * i + 7.0;
	   b_a = -b / a;
	   c_a = c / a;
	   (x1, x2) = FasterQuadratic(c, b, a);
	   (x3, x4) = FindRoots.Quadratic(c, b, a);
	   AssertComplexEqual(x1 + x2, b_a, 1e-14);
	   AssertComplexEqual(x1 * x2, c_a, 1e-14);
	   AssertComplexEqual(x3 + x4, x1 + x2, 1e-14);
	   AssertComplexEqual(x3 * x4, x1 * x2, 1e-14);
   }
}

[FrogGuaGuaGua]

FasterQuadratic is 3.8 times faster than FindRoots.Quadratic, and the latter does not handle the case of a=0, which is a bug.

[jkalias]

Isn't this the exact same implementation of FindRoots.Quadratic?

[jkalias]

I don't think these tests bring some value to the quality control of the code base. Performance benchmarks are indeed significant, but the results are not tested. And the correctness test is checking the same results, since the faster quadratic implementation is the same as the production code

[FrogGuaGuaGua]

I gave up committing and don't want to explain such an obvious problem anymore.