Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MATH-1660 FastMath/AccurateMath.scalb does not handle subnormal results #232

Open
wants to merge 2 commits into
base: master
Choose a base branch
from

Conversation

syenkoc
Copy link

@syenkoc syenkoc commented May 8, 2023

This PR contains a fix for this issue, for both double and float inputs.

I've also included some a unit tests that should provide 100% coverage.

Basic performance testing indicates that this implementation is somewhat faster than Math.scalb (and definitely faster than AccurateMath.scalb). However, I only tried this on OLD Intel hardware. This implementation uses only bit string manipulation, thus will never to any ops on subnormal values, which can be notoriously slow even on modern hardware.

Copy link
Contributor

@aherbert aherbert left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the updated code. The logic all seems correct and the rewrite is more maintainable between double and float by inserting all constants in the top of the method.

It may be useful to have a JMH benchmark between the old and new code for the case of normal-to-normal scaling. However we do not currently have a JMH module for benchmarks in this component so we can leave that for now. Can you insert a comment where this would apply within the method with a TODO, e.g.:

// TODO - speed test of a direct scaling by multiplication when scaleFactor in [-1022, 1023]

If the old method is faster then we could insert the direct multiplication at the top of the code.

I would not use an inline RNG. Just reuse an object that has a robust generation method. Pick one from the Commons RNG component and the UniformRandomProvider implementation can be easily changed in the future.

I prefer to have the tests hit all code paths without requiring randomness. So can you test MATH-1660 explicitly. The rest of the paths should be hit already by the current tests.

Note: You can switch to using JUnit 5 tests and assertions. We just haven't upgraded the old tests yet but new tests can use the new framework. As noted in the comments, it makes the custom assertBitPattern redundant.

The CI build was failing due to an error with the parent pom. This has been corrected in master so please rebase to collect the changes.


// first simple and fast handling when 2^n can be represented using normal numbers
if ((n > -1023) && (n < 1024)) {
return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By dropping this are we removing a faster code execution for this case? Decomposition of the input double is not required. I would think that a typical use of scalb is when you have a number that you think can be scaled to something useful (i.e. a normal to normal scaling). A non-typical use case would be to use scalb to create non-normal numbers. By this logic, this code path may be quite common in practice.

Under your control flow the same case would have: a check for non-finite; a check for non-zero exponent and then the statement to check the result exponent is within the allowed range. For a majority use case all these branches may be very predictable. The speed is then the multiplication plus the composition of the factor (old method), verses the decomposition and then recomposition of the double. It may be interesting to measure the two codes with a JMH benchmark.

@Test
public void doubleNormalValues() {

long m = 1;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do not use an inline RNG. Just use a robust UniformRandomProvider from the Commons RNG project with some fixed seed, e.g.

UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create(0xdeadbeaf);

assertBitPattern(StrictMath.scalb(x, scaleFactor), AccurateMath.scalb(x, scaleFactor));
assertBitPattern(StrictMath.scalb(-x, scaleFactor), AccurateMath.scalb(-x, scaleFactor));

m ^= m << 13;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove

assertBitPattern(StrictMath.scalb(x, scaleFactor), AccurateMath.scalb(x, scaleFactor));
assertBitPattern(StrictMath.scalb(-x, scaleFactor), AccurateMath.scalb(-x, scaleFactor));

// Marsaglia XOR pseudo-random shift.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Replace with a RNG object

assertBitPattern(StrictMath.scalb(x, scaleFactor), AccurateMath.scalb(x, scaleFactor));
assertBitPattern(StrictMath.scalb(-x, scaleFactor), AccurateMath.scalb(-x, scaleFactor));

m ^= m << 13;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Replace with a RNG object


Assert.assertEquals(Float.floatToIntBits(expected), Float.floatToIntBits(actual));
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove extra new line


@Test
public void floatSpecial() {

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove extra new line


long m = 1;
for(int scaleFactor = Double.MIN_EXPONENT * 2; scaleFactor <= (Double.MAX_EXPONENT - Double.MIN_EXPONENT); scaleFactor++) {
for(int count = 0; count < 50000; count++) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a lot of repeats. But you do not include the explicit case in MATH-1660 when rounding has sticky bits 0. I would add some explicit tests to ensure code paths are hit without requiring a lot of randomness. Something like this should hit cases of round-to-even with a trailing guard bit and further 'sticky' bit:

@ParameterizedTest
@MethodSource
void testScalb(double x, int n) {
    Assertions.assertEquals(StrictMath.scalb(x, n), AccurateMath.scalb(x, n));
    Assertions.assertEquals(StrictMath.scalb(-x, n), AccurateMath.scalb(-x, n));
}

static Stream<Arguments> testScalb() {
    Stream.Builder<Arguments> builder = Stream.builder();
    builder.add(Arguments.of(-0x1.8df4a353af3a8p-2, -1024));
    for (int i = 0; i < 8; i++) {
        builder.add(Arguments.of(Double.longBitsToDouble(511L << 52 | i), -512));
    }
    for (int i = 0; i < 8; i++) {
        builder.add(Arguments.of(Double.longBitsToDouble(1023L << 52 | i), -1024));
    }
    return builder.build();
}

The old CM code fails on 1023L << 52 | 2.

Note that you can call assertScalb from other methods too as this can become your default assertion method, e.g. in doubleSpecial do assertScalb(Double.NaN, scaleFactor); etc.

* @param expected Expected value.
* @param actual Actual value.
*/
private static void assertBitPattern(final double expected, final double actual) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not required. Using JUnit 5 this binary equality is performed by Assertions.assertEquals. Internally it just converts the double to a long using doubleToLongBits which collates all NaNs to a single long value.

/**
* Tests for scalb.
*/
public class ScalbTest {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The current AccurateMath tests current reside in commons-math-legacy-core module (due to an artefact of moving around classes for the module structure). For now I would move this to org.apache.commons.math4.legacy.core.jdkmath.AccurateMathScalbTest in that module so all the AccurateMath tests are in the same place.

@syenkoc
Copy link
Author

syenkoc commented May 9, 2023

OK thanks for the review, I should have time this weekend to make the required code changes.

In terms of performance, I tested this using a hacked up version of the old FastMathTestPerformance class. This new implementation definitely faster without the "fast" case code, about 30%.

Here is what I get, comparing StrictMath, Math and the two AccurateMath implementations across 200000000 runs:

Without "fast" case
StrictMath 2.797007 1.0
AccuratMathOld 3.641469 1.3019
Math 2.750160 0.9833
AccurateMathNew 1.4190 0.5073

With "fast" case
StrictMath 2.778882 1.0
AccuratMathOld 3.590121 1.2919
Math 2.709783 0.9751
AccurateMathNew 2.4186 0.8704

These tests were done on a 1.2GHz Intel Broadwell processor from circa 2015. I dumped the JIT-compiled code, but it's a bit hard to make sense of why it is faster.

I am unfamiliar with JMH, but I will try this weekend to get it going.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants