RFR: 8334755: Asymptotically faster implementation of square root algorithm [v47]

[Raffaello Giulietti]

On Mon, 29 Jul 2024 13:18:16 GMT, fabioromano1 <duke at openjdk.org> wrote:

>> I have implemented the Zimmermann's square root algorithm, available in works [here](https://inria.hal.science/inria-00072854/en/) and [here](https://www.researchgate.net/publication/220532560_A_proof_of_GMP_square_root).
>> 
>> The algorithm is proved to be asymptotically faster than the Newton's Method, even for small numbers. To get an idea of how much the Newton's Method is slow,  consult my article [here](https://arxiv.org/abs/2406.07751), in which I compare Newton's Method with a version of classical square root algorithm that I implemented. After implementing Zimmermann's algorithm, it turns out that it is faster than my algorithm even for small numbers.
>
> fabioromano1 has updated the pull request incrementally with one additional commit since the last revision:
> 
>   If the input is a square, then s0 == 0, so testing for non-zero remainder is redundant

Just some quick reviews.
Expect a more definitive review in the next days.

src/java.base/share/classes/java/math/MutableBigInteger.java line 550:

> 548:      */
> 549:     void safeRightShift(int n) {
> 550:         if (n >= bitLength()) {

The commit message for this reads `More accurate condition for MBI.safeRightShift()`.
If the old version works, please switch back. But if this is a genuine bug, then it needs a separate bug issue and PR.

src/java.base/share/classes/java/math/MutableBigInteger.java line 1891:

> 1889:         int shift = Integer.numberOfLeadingZeros(x.value[x.offset]) & ~1; // shift must be even
> 1890:         if ((x.intLen & 1) != 0)
> 1891:             shift += 32; // x.intLen must be even

Suggestion:

        int shift = (Integer.numberOfLeadingZeros(x.value[x.offset]) & ~1) + ((x.intLen & 1) << 32);

src/java.base/share/classes/java/math/MutableBigInteger.java line 1922:

> 1920:     }
> 1921: 
> 1922:     private static long ulongSqrt(long x) {

Since there is a `unisgnedLongCompare` already, it would be better to name this `unsignedLongSqrt`.

src/java.base/share/classes/java/math/MutableBigInteger.java line 1946:

> 1944:      * @implNote The implementation is based on Zimmermann's works available
> 1945:      * <a href="https://inria.hal.science/inria-00072854/en/">  here</a> and
> 1946:      * <a href="https://www.researchgate.net/publication/220532560_A_proof_of_GMP_square_root">  here</a>

The following variant should be preferred, as it has a readable figure on p. 21, whereas the same figure in the variant linked in the PR seems broken for some reason.
https://inria.hal.science/inria-00072113/document

src/java.base/share/classes/java/math/MutableBigInteger.java line 1948:

> 1946:      * <a href="https://www.researchgate.net/publication/220532560_A_proof_of_GMP_square_root">  here</a>
> 1947:      */
> 1948:     private MutableBigInteger[] sqrtRemZimmermann(int len, boolean needRemainder) {

This should be called `sqrtRemKaratsuba()`. This is the name chosen by Paul Zimmermann himself.

Also, I wonder if it wouldn't be simpler for `len` to represent the `int` length of the square root rather than the `int` length of the argument. It would be more consistent with the Bertot, Magaud, Zimmermann paper and might slightly simplify some later computation. But I'm not sure.

src/java.base/share/classes/java/math/MutableBigInteger.java line 1968:

> 1966:         final int halfLen = len >> 1;
> 1967:         // Recursive invocation
> 1968:         MutableBigInteger[] sr = sqrtRemZimmermann((halfLen & 1) == 0 ? halfLen : halfLen + 1, true);

Suggestion:

        final int halfLen = len >> 1;
        // Recursive invocation
        MutableBigInteger[] sr = sqrtRemZimmermann(halfLen + (halfLen & 1), true);

src/java.base/share/classes/java/math/MutableBigInteger.java line 2018:

> 2016:      * {@code this} number, ending at {@code blockIndex*blockLen} (exclusive).
> 2017:      */
> 2018:     private MutableBigInteger getBlockZimmermann(int blockIndex, int len, int blockLen) {

Should be named `getBlockKaratsuba()` or similar.

test/jdk/java/math/BigInteger/BigIntegerTest.java line 299:

> 297:         /* For every long value n in [0, 2^32) such that x == n * n,
> 298:          * n - 1 <= (long) Math.sqrt(x >= 0 ? x : x + 0x1p64) <= n
> 299:          * must be true.

Suggestion:

         * must be true.
         * This property is used to implement MutableBigInteger.unsignedLongSqrt().

test/micro/org/openjdk/bench/java/math/BigIntegerSquareRoot.java line 85:

> 83:             largeArray[i] = new BigInteger("" + (value / 1000));
> 84:             smallArray[i] = new BigInteger("" + hi);
> 85:         }

I think it would be more representative to have 5 arrays, from extra-small (XS) to extra-large (XL), with elements created by `new BigInteger(r.nextInt(maxLen), r)`, where `maxLen` is 64, 256, 1_024, 4_096, 16_384, resp., and drop the logic with the string parsing.

-------------

PR Review: https://git.openjdk.org/jdk/pull/19710#pullrequestreview-2205484016
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695548820
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695549452
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695549943
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695550713
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695551231
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695558116
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695558368
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695558651
PR Review Comment: https://git.openjdk.org/jdk/pull/19710#discussion_r1695559743