Reduce overflow via common denominators

Project.toml

@@ -1,7 +1,7 @@
 name = "Ratios"
 uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
 authors = ["Tim Holy <[email protected]>"]
-version = "0.4.1"
+version = "0.4.2"
 
 [deps]
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"

README.md

@@ -2,6 +2,76 @@
 
 [![Build Status](https://travis-ci.org/timholy/Ratios.jl.svg?branch=master)](https://travis-ci.org/timholy/Ratios.jl)
 
-This package provides types similar to Julia's `Rational` type, which make some sacrifices but have better computational performance.
+This package provides types similar to Julia's `Rational` type, which make some sacrifices but have better computational performance at the risk of greater risk of overflow.
 
 Currently the only type provided is `SimpleRatio(num, den)` for two integers `num` and `den`.
+
+Demo:
+
+```julia
+julia> x, y, z = SimpleRatio(1, 8), SimpleRatio(1, 4), SimpleRatio(2, 8)
+(SimpleRatio{Int}(1, 8), SimpleRatio{Int}(1, 4), SimpleRatio{Int}(2, 8))
+
+julia> x+y
+SimpleRatio{Int}(12, 32)
+
+julia> x+z
+SimpleRatio{Int}(3, 8)
+```
+
+`y` and `z` both represent the rational number `1//4`, but when performing arithmetic with `x`
+`z` is preferred because it has the same denominator and is less likely to overflow.
+
+To detect overflow, [SaferIntegers.jl](https://github.com/JeffreySarnoff/SaferIntegers.jl) is recommended:
+
+```julia
+julia> using Ratios, SaferIntegers
+
+julia> x, y = SimpleRatio{SafeInt8}(1, 20), SimpleRatio{SafeInt8}(1, 21)
+(SimpleRatio{SafeInt8}(1, 20), SimpleRatio{SafeInt8}(1, 21))
+
+julia> x + y
+ERROR: OverflowError: 20 * 21 overflowed for type Int8
+Stacktrace:
+[...]
+```
+
+[FastRationals](https://github.com/JeffreySarnoff/FastRationals.jl) is another package with safety and performance characteristics that lies somewhere between `SimpleRatio` and `Rational`:
+
+```julia
+julia> @benchmark x + y setup=((x, y) = (SimpleRatio(rand(-20:20), rand(2:20)), SimpleRatio(rand(-20:20), rand(2:20))))
+BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
+ Range (min … max):  1.727 ns … 28.575 ns  ┊ GC (min … max): 0.00% … 0.00%
+ Time  (median):     1.739 ns              ┊ GC (median):    0.00%
+ Time  (mean ± σ):   1.753 ns ±  0.445 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
+
+            ▂ ▃ ▃ ▃ ▄ ▆ ▇ █ ▆ ▅  ▅ ▇ ▄ ▁
+  ▂▁▂▁▃▁▅▁█▁█▁█▁█▁█▁█▁█▁█▁█▁█▁█▁▁█▁█▁█▁█▁▆▁▃▁▃▁▃▁▃▁▃▁▃▁▃▁▃▁▂ ▄
+  1.73 ns        Histogram: frequency by time        1.76 ns <
+
+ Memory estimate: 0 bytes, allocs estimate: 0.
+
+julia> @benchmark x + y setup=((x, y) = (FastRational(rand(-20:20), rand(2:20)), FastRational(rand(-20:20), rand(2:20))))
+BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
+ Range (min … max):  3.192 ns … 89.167 ns  ┊ GC (min … max): 0.00% … 0.00%
+ Time  (median):     3.215 ns              ┊ GC (median):    0.00%
+ Time  (mean ± σ):   3.307 ns ±  1.820 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
+
+    ▃█▆█▃▂
+  ▄███████▅▄▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▁▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂ ▃
+  3.19 ns        Histogram: frequency by time        3.45 ns <
+
+ Memory estimate: 0 bytes, allocs estimate: 0.
+
+julia> @benchmark x + y setup=((x, y) = (Rational(rand(-20:20), rand(2:20)), Rational(rand(-20:20), rand(2:20))))
+BenchmarkTools.Trial: 10000 samples with 996 evaluations.
+ Range (min … max):  22.385 ns … 81.269 ns  ┊ GC (min … max): 0.00% … 0.00%
+ Time  (median):     32.777 ns              ┊ GC (median):    0.00%
+ Time  (mean ± σ):   33.162 ns ±  4.743 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
+
+                ▁  ▇  ▄▂ ▁▆▃  ▂█▃ ▁ ▇   ▁    ▁
+  ▁▁▁▁▁▄▂▁▆▂▂█▄▃█▅▆█▇▇█████████████▅█▆▂▁█▇▁▂▁█▄▁▂▁▁▆▂▁▁▁▁▃▁▁▁ ▃
+  22.4 ns         Histogram: frequency by time        45.8 ns <
+
+ Memory estimate: 0 bytes, allocs estimate: 0.
+```

src/Ratios.jl

@@ -1,16 +1,51 @@
+"""
+`Ratios` provides `SimpleRatio`, a faster variant of `Rational`.
+Speed comes at the cost of greater vulnerability to overflow.
+
+API summary:
+
+- `r = SimpleRatio(num, den)` is equivalent to `num // den`
+- `common_denominator` standardizes a collection of `SimpleRatio`s to have the same denominator,
+  making some arithmetic operations among them less likely to overflow.
+"""
 module Ratios
 
 import Base: convert, promote_rule, *, /, +, -, ^, ==, decompose
 
 using Requires
 
-export SimpleRatio
+export SimpleRatio, common_denominator
 
 struct SimpleRatio{T<:Integer} <: Real
     num::T
     den::T
 end
 
+"""
+    SimpleRatio(num::Integer, den::Integer)
+
+Construct the equivalent of the rational number `num // den`.
+
+Operations with `SimpleRatio` are faster, but also more vulnerable to integer overflow,
+than with `Rational`. Arithmetic with `SimpleRatio` does not perform any simplification of the resulting ratios.
+The best defense against overflow is to use ratios with the same denominator: in such cases,
+`+` and `-` will skip forming the product of denominators. See [`common_denominator`](@ref).
+
+If overflow is a risk, consider constructing them using SaferIntegers.jl.
+
+# Examples
+
+```julia
+julia> x, y, z = SimpleRatio(1, 8), SimpleRatio(1, 4), SimpleRatio(2, 8)
+(SimpleRatio{$Int}(1, 8), SimpleRatio{$Int}(1, 4), SimpleRatio{$Int}(2, 8))
+
+julia> x+y
+SimpleRatio{$Int}(12, 32)
+
+julia> x+z
+SimpleRatio{$Int}(3, 8)
+```
+"""
 SimpleRatio(num::Integer, den::Integer) = SimpleRatio(promote(num, den)...)
 
 convert(::Type{BigFloat}, r::SimpleRatio{S}) where {S} = BigFloat(r.num)/r.den
@@ -32,8 +67,10 @@ Rational{T}(r::SimpleRatio{S}) where {T<:Integer, S<:Integer} = convert(T, r.num
 /(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den, y.num)
 +(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den + y.num, y.den)
 -(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den - y.num, y.den)
-+(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)
--(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)
++(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num + y.num, x.den) :
+                                                     SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)
+-(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num - y.num, x.den) :
+                                                     SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)
 ^(x::SimpleRatio, y::Integer) = SimpleRatio(x.num^y, x.den^y)
 
 -(x::SimpleRatio{T}) where {T<:Signed} = SimpleRatio(-x.num, x.den)
@@ -61,6 +98,27 @@ end
 
 decompose(x::SimpleRatio) = x.num, 0, x.den
 
+"""
+    common_denominator(x::SimpleRatio, ys::SimpleRatio...)
+
+Return the equivalent of `(x, ys...)` but using a common denominator.
+This can be useful to avoid overflow.
+
+!!! info
+    This function is quite slow.  In performance-sensitive contexts where the
+    ratios are constructed with literal constants, it is better to ensure a common
+    denominator at the time of original construction.
+"""
+function common_denominator(x::SimpleRatio, ys::SimpleRatio...)
+    all(y -> y.den == x.den, ys) && return (x, ys...)
+    cd = gcd(x.den, map(y -> y.den, ys)...)        # common divisor
+    # product of the denominators
+    pd = Base.Checked.checked_mul(cd, mapreduce(z -> z.den ÷ cd, Base.Checked.checked_mul, (x, ys...)))
+    return map((x, ys...)) do z
+        SimpleRatio(z.num * (pd ÷ z.den), pd)
+    end
+end
+
 function __init__()
     @require FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" begin
         using .FixedPointNumbers: FixedPoint, Fixed, Normed, rawone

test/runtests.jl

@@ -26,6 +26,13 @@ using FixedPointNumbers
     @test r == 0.5
     @test 0.5 == r
 
+    r3 = SimpleRatio(7,3)
+    @test r2+r3 === SimpleRatio(9,3)
+    @test r2-r3 === SimpleRatio(-5,3)
+
+    r8 = SimpleRatio{Int8}(1, 20)
+    @test r8 + r8 == SimpleRatio(1, 10)
+
     @test_throws OverflowError -SimpleRatio(0x02,0x03)
 
     @test r + SimpleRatio(0x02,0x03) == SimpleRatio(7,6)
@@ -42,4 +49,13 @@ using FixedPointNumbers
     @test SimpleRatio(5,3) * -0.03Q0f7 == SimpleRatio{Int}(rationalize((5.0*(-0.03Q0f7))/3))
     r = @inferred(SimpleRatio(0.75Q0f7))
     @test r == 3//4 && r isa SimpleRatio{Int16}
+
+    # common_denominator
+    @test common_denominator(SimpleRatio(2,7), SimpleRatio(3,11), SimpleRatio(-1,5)) ===
+        (SimpleRatio(2*11*5,385), SimpleRatio(3*7*5,385), SimpleRatio(-1*7*11,385))
+    @test common_denominator(SimpleRatio(5,12), SimpleRatio(4,15), SimpleRatio(-1,9)) ===
+        (SimpleRatio(75,180), SimpleRatio(48,180), SimpleRatio(-20,180))
+    @test_throws OverflowError common_denominator(SimpleRatio{Int8}(1, 20), SimpleRatio{Int8}(2, 21))
+    @test common_denominator(SimpleRatio{Int8}(1, 20), SimpleRatio{Int8}(3, 20)) ===
+        (SimpleRatio{Int8}(1, 20), SimpleRatio{Int8}(3, 20))
 end