Automatically divide by gcd when needed

[MilesCranmer]

I am wondering whether Ratios.jl might be able to automatically divide by the gcd whenever an overflow would otherwise occur. Perhaps this could be done through an extension to SafeIntegers.jl?

With @oscardssmith we've been looking at using SimpleRatio{Int8} for storing dimensions in DynamicQuantities.jl: SymbolicML/DynamicQuantities.jl#4. In other words, the power of each physical dimension would be a SimpleRatio, like (length)^(::SimpleRatio{Int8}).

However, this can easily hit overflows, especially if a user is performing many calculations on a single quantity. So I am wondering if there is a way we automatically divide by the gcd, but only when needed. (Rather than every calculation, with Rational, or never, with the normal SimpleRatio).

[timholy]

We'd have to check the performance implications pretty carefully, but that's something we could consider. One thing to know: the reason for the existence of this package is Interpolations.jl, and I won't make changes that cause it significant harm. However, as with the discussion of FixedPointNumbers, perhaps one could add another type parameter and make the arithmetic configurable.

[i9e1]

In between: Always divide by greatest common power of two.

function reduce_base2(R)
  # den > 0 => trailing_zeros is well defined
  t = trailing_zeros(R.num | R.den)
  
  SimpleRatio(R.num>>t, R.den>>t)
end

This should cost <10 simple instructions. No branching at all.

[timholy]

Right, but interpolation is f*a + (1-f)*b. That's much less than 10 instructions. Now, maybe it would never be triggered, but this is what I mean when I say we'd have to check the performance implications for Interpolations.jl quite carefully.

If we find that there's a conflict in goals, keep in mind that the whole reason I called this SimpleRatio is that we could also add SlightlySlowerButSaferThanSimpleRatio to this package, duplicating the code for SimpleRatio and then adding these extra operations.

[i9e1]

For instruction read RISC-like machine instruction. All operations will get a small but constant extra cost. As an advantage, we should safe conditional evaluation (which can take longer than a few instructions on most modern machines).

Some weeks ago I tried profile my approach against the current version implemented here.
When interpolating random data the performance looks very similiar, which seems a little bit suspicious.

Than I tried to compare the errors with each implementation against exact Rational{BigInt}. For the next few months I dont have the time to do a proper comparison, but if anyone is interested I will keep it on my personal backlog.

[MilesCranmer]

So what I ended up doing in DynamicQuantities.jl was just having a fixed denominator rational number (@oscardssmith suggestion). (Kind of like FixedPointNumbers but a different emphasis)

https://github.com/SymbolicML/DynamicQuantities.jl/blob/12ecee20fa3eb90ed842f2ddfb9f111c230d2c40/src/fixed_rational.jl

The default type used to represent the exponent of a physical dimension is FixedRational{Int32, 2^4 * 3^2 * 5^2 * 7} (an Int32 divided by 25,200). But you can even get away with, e.g., FixedRational{Int8, 6} for many calculations.

It's fast at the things that are important for physical quantities, like equality checks (=> check equality of numerators = fast), and addition/subtraction (=> add the numerators = fast). But it's slower at multiplication and division because it has to divide by the denominator, but those are less common here because it's only used in power laws. It also avoids having to compute the gcd frequently which is good.

It's a bit different in terms of what numbers it can represent. e.g., 1/11 is not possible to represent exactly with the above defaults (a tradeoff in what numbers are representable with a fixed number of bits). However, it's fine for unit calculations.

I'd be happy to move it somewhere else like Ratios.jl if @timholy is interested and thinks it would be of general interest. Having it in a specific package for rational numbers seems like a more maintainable strategy than having an internal type I have to verify and test.

[timholy]

Seems reasonable. My only concern is whether someone might create an unlimited number of denominators? In your application, are there only a few relevant denominators or is it "anything goes"? (Type explosions are a bit concerning...)

[oscardssmith]

For representing the exponent of units I think people will create essentially 1 denominator (the default one). It's very rare to need meters^129/127ths for example.

[timholy]

I think we could add that to this package if you'd like. We could include a warning in the docs about the potential for abuse and encourage people to limit the number of different denominators.

[MilesCranmer]

I should mention that there are no methods available for FixedRationals that have different denominators, so to run into type explosions you would have to do so explicitly.

(The design assumption is that someone would pick a single denominator and use it throughout their code.)