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

[fabioromano1]

On Thu, 1 Aug 2024 19:10:51 GMT, Raffaello Giulietti <rgiulietti at openjdk.org> wrote:

>> fabioromano1 has updated the pull request incrementally with one additional commit since the last revision:
>> 
>>   Last small changes
>
> I guess the overhead is negligible when compared to the arithmetic operation (shifts, divisions, etc.).
> Also, the maximal stack depth for the recursion itself is quite limited and under control.
> If at all, I would rather postpone that effort to a followup PR, and only if there are noticeable advantages without compromising readability.

@rgiulietti In case, the code should look like this:

    /**
     * Assumes {@code intLen > 2 && intLen % 2 == 0
     * && Integer.numberOfLeadingZeros(value[offset]) <= 1}
     */
    private MutableBigInteger[] sqrtRemKaratsuba(boolean needRemainder) {
        MutableBigInteger sqrt, rem;
        // Base case
        {
            long x = ((value[offset] & LONG_MASK) << 32) | (value[offset + 1] & LONG_MASK);
            long s = unsignedLongSqrt(x);

            // Allocate sufficient space to hold the final square root
            sqrt = new MutableBigInteger(new int[intLen >> 1]);

            // Place the partial square root
            sqrt.intLen = 1;
            sqrt.value[0] = (int) s;

            rem = new MutableBigInteger(x - s * s);
        }

        // Iterative step
        int len;
        for (len = 4; len < intLen; len <<= 1) {
            final int blockLen = len >> 2;
            MutableBigInteger dividend = rem;
            dividend.shiftAddDisjoint(getBlockForSqrt(1, len, blockLen), blockLen);

            // Compute dividend / (2*sqrt)
            MutableBigInteger q = new MutableBigInteger();
            MutableBigInteger u = dividend.divide(sqrt, q);
            if (q.isOdd())
                u.add(sqrt);
            q.rightShift(1);

            sqrt.shiftAdd(q, blockLen);
            // Corresponds to ub + a_0 in the paper
            u.shiftAddDisjoint(getBlockForSqrt(0, len, blockLen), blockLen);
            BigInteger qBig = q.toBigInteger(); // Cast to BigInteger to use fast multiplication
            MutableBigInteger qSqr = new MutableBigInteger(qBig.multiply(qBig).mag);

            rem = u;
            if (rem.subtract(qSqr) < 0) {
                MutableBigInteger twiceSqrt = new MutableBigInteger(sqrt);
                twiceSqrt.leftShift(1);

                // Since subtract() performs an absolute difference, to get the correct algebraic sum
                // we must first add the sum of absolute values of addends concordant with the sign of rem
                // and then subtract the sum of absolute values of addends that are discordant
                rem.add(ONE);
                rem.subtract(twiceSqrt);
                sqrt.subtract(ONE);
            }
        }

        // Final step
        final int blockLen = len >> 2;
        MutableBigInteger dividend = rem;
        dividend.shiftAddDisjoint(getBlockForSqrt(1, len, blockLen), blockLen);

       MutableBigInteger q = new MutableBigInteger();
        MutableBigInteger u = dividend.divide(sqrt, q);
        if (q.isOdd())
            u.add(sqrt);
        q.rightShift(1);

        sqrt.shiftAdd(q, blockLen);
        u.shiftAddDisjoint(getBlockForSqrt(0, len, blockLen), blockLen);
        BigInteger qBig = q.toBigInteger();
        MutableBigInteger qSqr = new MutableBigInteger(qBig.multiply(qBig).mag);

        if (needRemainder) {
            rem = u;
            if (rem.subtract(qSqr) < 0) {
                MutableBigInteger twiceSqrt = new MutableBigInteger(sqrt);
                twiceSqrt.leftShift(1);

                rem.add(ONE);
                rem.subtract(twiceSqrt);
                sqrt.subtract(ONE);
            }
        } else {
            rem = null;
            if (u.compare(qSqr) < 0)
                sqrt.subtract(ONE);
        }

        return new MutableBigInteger[] { sqrt, rem };
    }

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

PR Comment: https://git.openjdk.org/jdk/pull/19710#issuecomment-2263853409