Really bad regime selection on quadratic formula

[msullivan]

I tried the web demo out on the two cases of the quadratic formula, and the negative case (-b - sqrt(b^2 - 4*a*c)) / (2a) produced this output: https://herbie.uwplse.org/demo/7c3b374107c5177485bdae893a16bc32fd361b84.b7c88a7d879b959f14e63c59c9c933cbdc874ea9/graph.html#reproduce

This features the answer -0.5 * (2 * (c/b)) (the result of taylor expanding around -inf) for b < -something*10^-93. This produces super wrong results for totally normal negative b values (for example, at a=1, b=-1, c=-6, it produces -6 instead of the correct -2).

When I tried the positive case, it worked great, splitting into 4 regimes: https://herbie.uwplse.org/demo/500fb1a758f004adc50b440bb74877e556b2c3b7.b7c88a7d879b959f14e63c59c9c933cbdc874ea9/graph.html.

[pavpanchekha]

You know, that's pretty interesting and I'll look into it. (Both quadratics are in our nightlies and generally do equally well, but looks like with this seed Herbie didn't find the right branch to add.)

[kamiyo]

@msullivan Just to let you know, if you set pre-conditions, it works much better. I limited a, b, and c to go from -1e6 to 1e6, and it gave me much better results, but still not ideal.

[bksaiki]

Just checked this. This still happens although rarely. For seeds 0-29 (inclusive), bad regime selection only happens once (seed 12).

[switham]

I (independently) tried the -b +sqrt(...) case and got the terrible answer below, three times-- maybe it just cached the answer it already gave me? By the way, the error graph indicates my version and herbie's version are both pretty accurate on ordinary arguments... !? My test cases have a, b and c uniform-random between -10 and +10 (only using cases where the determinant >= 0 and abs(a) > .01).

Numerical Recipes in C 2nd Ed. p. 184 has this method:

q = -.5 * (b + sgn(b) * sqrt(b*b - 4*a*c))
solutions = [q / a,  c / q]

Here is my input and herbie's output:

double code(double b, double a, double c) {
	return (-b + sqrt((b * b) - ((4.0 * a) * c))) / (2.0 * a);
}
↓
double code(double b, double a, double c) {        /* Mainly, this works iff b < 0.  */
	double tmp;
	if (b <= -5.907209789018034e+89) {
		tmp = -b / a;
	} else if (b <= 4.782539630045316e-37) {
		tmp = ((sqrt(fma(a, (c * -4.0), (b * b))) - b) * 0.5) / a;
	} else {
		tmp = -(c / b);
	}
	return tmp;
}

[pavpanchekha]

Interesting. I know that Herbie can produce the NR answer, and in fact usually does so. In the original issue, it does (mostly) find the NR answer, but decides not to use it. This feels like an issue in regimes. I know that I've fixed some regimes bugs recently, but I don't think either of them would fix the original post, or this one. It needs more triage. @bksaiki's seed survey is definitely helpful, we should do more of those for Herbie.