JDK 14 RFR of JDK-8233452: java.math.BigDecimal.sqrt() with RoundingMode.FLOOR results in incorrect result

[Joe Darcy]

Hello,

Please review the changes to fix

     JDK-8233452: java.math.BigDecimal.sqrt() with RoundingMode.FLOOR 
results in incorrect result
     http://cr.openjdk.java.net/~darcy/8233452.0/

Some background on the problem and fix.

The core of the BigDecimal.sqrt method is a Newton-Raphson loop. Given a 
sufficiently good initial value and given various other side conditions 
which hold for the square root function, the iteration converges 
quadratically, in this case to the square root of the receiver 
BigDecimal. The iteration starts using an initial value computed using 
Math.sqrt, which is sufficiently good for convergence.

Colloquially, quadratic convergence is described as "double the number 
of digits right" at each step. More formally, the error bound is reduced 
such that the error bound of the (k+1)st iteration is a function of the 
square of the error bound of the k-th iteration where the error bound is 
<< 1.

There are three basic rounding policies on BigDecimal operations:

1) Return an exact result (precision = 0 or RoundingMode.UNNECESSARY).
2) Return the nearest result, choosing a policy for ties (HALF_DOWN, 
HALF_EVEN, HALF_UP).
3) Directed rounding  (UP, DOWN, CEILING, FLOOR)

Since sqrt is only defined for non-negative values, UP and CEILING are 
equivalent for that function as are DOWN and FLOOR.

The directed roundings can have an error as much as 1 ulp (unit in the 
last place) as opposed to an error of only 1/2 ulp under the to-nearest 
rounding modes.  For example, say the true answer is between 1 and 2 
when rounding to integers. Under a to-nearest policy the largest 
difference between the true numerical answer and the returned results is 
at most 0.5; the true answer of 1.5 would round to either 1 or 2 
depending on how ties were handled. In contrast, under a directed 
rounding the largest difference can be nearly 1. If the true result is 
1.9999999.... and rounding down, the returned results will be 1 and the 
difference will be 0.9999999....

The Newton iteration reduces the error at each step and, as currently 
written, it can settle on a value like 2, which is closer to the actual 
answer, but *incorrect* according to the requested rounding mode.

There are a few ways to fix this one. One is to run the Newton loop to a 
higher-precision where these finer distinctions can be teased out. If 
the result is to have a p digits of accuracty, if the loop is run to 
have 2p+2 digits of precision, the correct final answer will be 
computed. (Proof of this property for floating-point arithmetic is 
discussed in various numerical references.) This approach has the 
disadvantage of running under higher precision at some performance cost. 
A reference implementation using this technique was added to the 
regression test.

Another approach is to detect when the result is too large/too small and 
then subtract/add an ulp as a fix-up. This is the approach taken for 
BigDecimal.sqrt. It should be a relatively inexpensive detection and 
fixup using a multiply, a compare, and possibility an add/subtract  as 
opposed to another divide operation under higher precision from going 
around the loop again. If enabled, the existing assertions in BigDecimal 
catch the erroneous answer computed by the existing code since it is 
detected as not bracketing the true result. The assertions pass when the 
need code is run against the previously problematic cases.

An approach not explored for this fix would be to arrange for the 
iteration to start from the "right" side of the number line depending on 
the rounding mode and then monotonically increase/decrease to approach 
the square root. For example, if rounding down, to start from a value 
greater than then root and then have each iteration decrease toward the 
final result. I believe this is technically possible, but would require 
some additional analysis to setup.

Thanks,

-Joe