Primitive double to RationalNumber about ojalgo HOT 28 CLOSED

This post on stackoverflow.com might provide a starting point:

https://stackoverflow.com/questions/13222664/convert-floating-point-number-into-a-rational-number-in-java

from ojalgo.

But that starting point highlights the problem of big exponents; won't it be better if we kept numerator as a pair of exponent and mantissa? It's still a proper rational number, and all the operations are well defined for it. Admittedly, we get to the land of having two kinds on rational numbers, and growing complexity, but it doesn't sound as a particularly high price to pay.

from ojalgo.

And the big-exponents problem isn't even that tricky to reproduce; if we try 1e7 as a parameter, we get the value of 5368709120000000 * (1<<23) / 4503599627370496 and the numerator is already overflown.

from ojalgo.

Dividing by gcd fixes the problem with 1e7 specifically, but leaves a question of testing:

• what are the "important" values?
• do we want 0.1 to be (3602879701896397 / 36028797018963968) or (1 / 10)?

from ojalgo.

I'd say the resulting rational number should be as exact as possible.

For the valueOf method 0.1 should definitely be represented as (1 / 10). Hard to explain why it should be anything else.

For the add/subtract/multiply/divide methods efficient execution is more important.

Do you have an implementations to test?

from ojalgo.

The implementation, which is quite naive, is at https://github.com/optimatika/ojAlgo/compare/master...alf239:issue49?expand=1; it fails some 233 tests (which I guess answers some of my questions), so I wasn't planning to push it just yet.

As for 0.1d, it's 11111110111001100110011001100110011001100110011001100110011010, which is

• sign: 0, positive
• exponent: 1111111011 (subtracting the bias of 1023, we get 0b1111111011 - 1023 = -4)
• fraction: 1001100110011001100110011001100110011001100110011010 standing for binary (sic!) 1.1001100110011001100110011001100110011001100110011010 which is 0b11001100110011001100110011001100110011001100110011010 * 2**(-52) == 7205759403792794 / 4503599627370496

Now 7205759403792794 and 4503599627370496 have a GCD of 2, so while 7205759403792794 / 4503599627370496 == 1.6 in the normal world, our precise world tells us that it's actually 3602879701896397 / 2251799813685248, and if we account for the exponent (2**-4), we get 3602879701896397 / 36028797018963968, which is almost but not quite equal to 0.1:

>>> 3602879701896397.0 / 36028797018963968
0.1

I will take a look at the BigDecimal's trick around 0.1, but the naïve approach leads us to the fraction above.

from ojalgo.

Updated to use approach from https://stackoverflow.com/a/13224190/562388 — it's less readable, but follows the javadoc, which is a great thing.

from ojalgo.

Looks clean.

I don't follow the logic entirely: how many test cases fail now?

I assume to minimise the number of failing tests one should try to mimic the behaviour of the BigDecimal/BigInteger based implementation - the part with the shift.

from ojalgo.

There are 233 failing tests at the moment; I will look closer today or tomorrow.

e.g.

# CompareImplementations.testASINH
junit.framework.AssertionFailedError: Quaternion vs Rational 
Expected :-0.0711857704764945
Actual   :-0.07665118472290924

# CompareImplementations.testMULTIPLY
junit.framework.AssertionFailedError: Quaternion vs Rational 
Expected :-0.8976069672109909
Actual   :0.8976069672109909

The latter gives hope that there's some ridiculous bug hiding, stay tuned.

from ojalgo.

Regarding the number of failing tests, it's hard to tell: the tests are randomised, and the seed isn't fixed, so what fails one time can easily pass next time.

from ojalgo.

Oh, I see, so you jump to BigIntegers and then divide both numerator and denominator by 2 until they fit safely? Interesting, that may help in many cases.

from ojalgo.

Why do you use

        if (retDenom.signum() > 0) {
            retNumer = retNumer.negate();
            retDenom = retDenom.negate();
        }

shouldn't it be

        if (retDenom.signum() < 0) {
            retNumer = retNumer.negate();
            retDenom = retDenom.negate();
        }

instead?

from ojalgo.

The above was indeed a bug responsible for CompareImplementations.testMULTIPLY failures and quite a few others. Fixed now, checking the rest.

from ojalgo.

org.ojalgo.matrix.BasicMatrixTest#testGetCondition apparently exposes another problem: sometimes,

    public Array1D<Double> getSingularValues() {

        if ((mySingularValues == null) && this.isComputed()) {
            mySingularValues = this.makeSingularValues();
        }

        return mySingularValues;
    }

appears to be in situation where mySingularValues == null and not this.isComputed(), so the result is null, and the test explodes with NPE. Wonder why it isn't computed...

from ojalgo.

So it's there as the matrix is not positive definite, thus cannot be decomposed.

from ojalgo.

Also, RationalNumber.POSITIVE_INFINITY.add(RationalNumber.POSITIVE_INFINITY) throws an exception. Similar problem for NaN, too.

from ojalgo.

Is the problem the zero denominator in POSITIVE_INFINITY and NaN? May have to review parts of the code to make sure it handles zero denominators.

The condition number is derived from the singular values decomposition. That should work regardless of "positive definite" or not.

That must be a bug related the retDenom.signum(). The sign of the rational number should always be encoded in the numerator.

from ojalgo.

The problem is with POSITIVE_INFINITY specifically. POSITIVE_INFINITY + POSITIVE_INFINITY = 1 / 0 + 1 / 0 = 2 / 0 = throw new ArithmeticException()

NaNs create a problem with doubleValue which is calculated via BigDecimal, and

The BigDecimal class provides no representation for NaN, +∞ or -∞.
(https://stackoverflow.com/a/28066001/562388)

from ojalgo.

Funnily enough, one of the failures is in the CompareImplementations.testMIN

Now the actual case it fails on is,

junit.framework.AssertionFailedError: Big vs Complex, -0.9645598407999487, 0.7898303544540544 
Expected :-0.9645598407999487
Actual   :0.7898303544540544

And it fails because for complex numbers we compare "norms", which are insensitive to the sign.

from ojalgo.

Raised PR #74

from ojalgo.

I'm surprised that tests pass repeatedly on the CI server, but so be it. Locally, I have three tests that fail regularly (though not always, it seems).

from ojalgo.

A better approach would be to use rational fraction aproximation, as in https://uk.mathworks.com/help/matlab/ref/rat.html?s_tid=gn_loc_drop - after all, we're not that bound to powers of 2 and 10.

Alternatively, we can keep the binary exponent around. That would let us keep the whole accuracy of double and then some.

from ojalgo.

It's a painful decision by the way. If we go with the binary exponent, we match doubles perfectly, but then 0.1 will be 7205759403792794 / 4503599627370496 just because of the limits of the mantissa.

On the other hand, if we go with the rational approximation, we'll have that 1 / 10 (haven't checked), but we'll lose any number greater than Long.MAX_VALUE, i.e. 2^63-1.

Do we want to use both and fall back to the exact representation when we fall out of range? Do we want to follow Mathematica's route and go with BigRationals instead?

from ojalgo.

Everything in the org.ojalgo.scalar package should be primitives based. RationalNumber was previously implemented using BigInteger. We're moving away from that.

Having 0.1 represented as 7205759403792794 / 4503599627370496 seems a horrible solution. Don't you agree?

Rational approximation is really interesting (perhaps combined with some way to store magnitude separately). I think it's important to not be too clever here. When users choose to work with RationalNumber it is important that they get something that "is" a rational number.

from ojalgo.

Good point regarding using primitives, fair enough.

As for 0.1, I'm a bit puzzled here. I admit that the expectation is indeed to have 0.1 equal to 1/10 - people tend not to know that 0.1 is actually 0b1.1001100110011001100110011001100110011001100110011010 * 2^{-4}, and are not likely to appreciate us forcing this knowledge on them.

Still, every now and again that abstraction leaks, for example for 0.1 + 0.2:

>>> 0.1 + 0.2
0.30000000000000004
>>> 

Mathematica deals with that problem by clearly separating the exact math (including rationals) and the approximate math; multiplying a rational by double gives one a double, not a rational.

In[1]:= 1 / 10

Out[1]= 1/10

In[2]:= 0.1

Out[2]= 0.1

In[4]:= Rationalize[0.1]

Out[4]= 1/10

In[5]:= 0.1 * (1/10)

Out[5]= 0.01

In[6]:= Rationalize[0.01]

Out[6]= 1/100

In[7]:= 0.1 * 10

Out[7]= 1.

In[8]:= (1/10) * 10


Out[8]= 1

In[9]:= 0.1 / (1/10)

Out[9]= 1.

If we follow this approach, then we should probably make the rational approximation a required call and retire the methods in the original post, as they invite the user in the land of approximation even though we promised it would be a rational.

from ojalgo.

OK, I think maybe we should try do disregard how things "look" and focus on what works best in terms of the computations.

Regarding large numbers: When the denominator gets large (compared to the numerator) then it doesn't seem problematic to approximate this with "zero". But, when it is the numerator that gets large then maybe sometimes returning the largest possible number (denominator 1) is preferable to returning infinity (denominator 0). Once you return infinity or NaN all further operations are defunct.

BTW I want to release v45 as soon as possible.

from ojalgo.

Let's see how the rational approximation approach works. I've drafted it as described, no optiization, nothing, and it seems that so far it... just works.

Once you return infinity or NaN all further operations are defunct.

But then if I don't, the result is incorrect, and heavily biased. For example, instead of 2^100 - 2^102 == -3/4 * 2^102 we would get Long.MAX_VALUE - Long.MAX_VALUE = ... well, either 0 or Long.MAX_VALUE depending on how exactly we do the subtraction.

from ojalgo.

Run complete. Total time: 00:01:21

Benchmark Mode Cnt Score Error Units
RationalCreation.rational thrpt 20 81612.897 ± 1092.180 ops/s
RationalCreation.valueOf thrpt 20 132289752.456 ± 10159005.133 ops/s

from ojalgo.