Divide units instead of basefactors in convfact and error on floati…

…ng-point over-/underflow (#648)

src/conversion.jl

@@ -1,38 +1,26 @@
 """
     convfact(s::Units, t::Units)
-Find the conversion factor from unit `t` to unit `s`, e.g., `convfact(m,cm) == 1//100`.
+Find the conversion factor from unit `t` to unit `s`, e.g., `convfact(m, cm) == 1//100`.
 """
 @generated function convfact(s::Units, t::Units)
-    sunits = s.parameters[1]
-    tunits = t.parameters[1]
-
     # Check if conversion is possible in principle
-    sdim = dimension(s())
-    tdim = dimension(t())
-    sdim != tdim && throw(DimensionError(s(),t()))
-
-    # first convert to base SI units.
-    # fact1 is what would need to be multiplied to get to base SI units
-    # fact2 is what would be multiplied to get from the result to base SI units
-
-    inex1, ex1 = basefactor(t())
-    inex2, ex2 = basefactor(s())
-
-    a = inex1 / inex2
-    ex = ex1 // ex2     # do overflow checking?
+    dimension(s()) != dimension(t()) && throw(DimensionError(s(), t()))
 
-    tens1 = mapreduce(tensfactor, +, tunits; init=0)
-    tens2 = mapreduce(tensfactor, +, sunits; init=0)
-
-    pow = tens1-tens2
+    # use absoluteunit because division is invalid for AffineUnits;
+    # convert to FreeUnits first because absolute ContextUnits might still
+    # promote to AffineUnits
+    conv_units = absoluteunit(FreeUnits(t())) / absoluteunit(FreeUnits(s()))
+    inex, ex = basefactor(conv_units)
+    pow = tensfactor(conv_units)
+    inex_orig = inex
 
     fpow = 10.0^pow
-    if fpow > typemax(Int) || 1/(fpow) > typemax(Int)
-        a *= fpow
+    if fpow > typemax(Int) || 1/fpow > typemax(Int)
+        inex *= fpow
     else
         comp = (pow > 0 ? fpow * numerator(ex) : 1/fpow * denominator(ex))
         if comp > typemax(Int)
-            a *= fpow
+            inex *= fpow
         else
             ex *= (10//1)^pow
         end
@@ -41,9 +29,15 @@ Find the conversion factor from unit `t` to unit `s`, e.g., `convfact(m,cm) == 1
     if ex isa Rational && denominator(ex) == 1
         ex = numerator(ex)
     end
-    a ≈ 1.0 ? (inex = 1) : (inex = a)
-    y = inex * ex
-    :($y)
+
+    result = (inex ≈ 1.0 ? 1 : inex) * ex
+    if fp_overflow_underflow(inex_orig, result)
+        throw(ArgumentError(
+            "Floating point overflow/underflow, probably due to large " *
+            "exponents and/or SI prefixes in units"
+        ))
+    end
+    return :($result)
 end
 
 """

src/units.jl

@@ -194,7 +194,7 @@ end
 function basefactor(inex, ex, eq, tens, p)
     # Sometimes (x::Rational)^1 can fail for large rationals because the result
     # is of type x*x so we do a hack here
-    function dpow(x,p)
+    function dpow(x, p)
         if p == 0
             1
         elseif p == 1
@@ -232,19 +232,23 @@ function basefactor(inex, ex, eq, tens, p)
             # Note that sometimes x^1 can cause an overflow error if x is large because
             # of how power_by_squaring is implemented for Rationals, so we use dpow.
             x = dpow(eq*ex*(10//1)^tens, p)
-            return (inex^p, isinteger(x) ? Int(x) : x)
+            result = (inex^p, isinteger(x) ? Int(x) : x)
         else
             x = dpow(ex*(10//1)^tens, p)
-            return ((inex*eq)^p, isinteger(x) ? Int(x) : x)
+            result = ((inex * eq)^p, isinteger(x) ? Int(x) : x)
         end
     else
         if eq_is_exact && can_exact2
-            x = dpow(eq,p)
-            return ((inex * ex * 10.0^tens)^p, isinteger(x) ? Int(x) : x)
+            x = dpow(eq, p)
+            result = ((inex * ex * 10.0^tens)^p, isinteger(x) ? Int(x) : x)
         else
-            return ((inex * ex * 10.0^tens * eq)^p, 1)
+            result = ((inex * ex * 10.0^tens * eq)^p, 1)
         end
     end
+    if fp_overflow_underflow(inex, first(result))
+        throw(ArgumentError("Floating point overflow/underflow, probably due to large exponent ($p)"))
+    end
+    return result
 end
 
 """
@@ -257,18 +261,21 @@ reference units.
 """
 @inline basefactor(x::Unit{U}) where {U} = basefactor(basefactors[U]..., 1, 0, power(x))
 
-function basefactor(x::Units{U}) where {U}
+function basefactor(::Units{U}) where {U}
     fact1 = map(basefactor, U)
-    inex1 = mapreduce(x->getfield(x,1), *, fact1, init=1.0)
-    float_num = mapreduce(x->float(numerator(getfield(x,2))), *, fact1, init=1.0)
-    float_den = mapreduce(x->float(denominator(getfield(x,2))), *, fact1, init=1.0)
-    can_exact = (float_num < typemax(Int))
-    can_exact &= (float_den < typemax(Int))
+    inex1 = mapreduce(first, *, fact1, init=1.0)
+    float_num = mapreduce(x -> float(numerator(last(x))), *, fact1, init=1.0)
+    float_den = mapreduce(x -> float(denominator(last(x))), *, fact1, init=1.0)
+    can_exact = float_num < typemax(Int) && float_den < typemax(Int)
     if can_exact
-        return inex1, mapreduce(x->getfield(x,2), *, fact1, init=1)
+        result = (inex1, mapreduce(last, *, fact1, init=1))
     else
-        return inex1*float_num/float_den, 1
+        result = (inex1 * (float_num / float_den), 1)
+    end
+    if any(fp_overflow_underflow(first(x), first(result)) for x in fact1)
+        throw(ArgumentError("Floating point overflow/underflow, probably due to a large exponent in some of the units"))
     end
+    return result
 end
 
 Base.broadcastable(x::Units) = Ref(x)

src/utils.jl

@@ -253,3 +253,6 @@ struct AffineError <: Exception
 end
 
 Base.showerror(io::IO, e::AffineError) = print(io, "AffineError: $(e.x)")
+
+fp_overflow_underflow(input, result) =
+    isfinite(input) && !isfinite(result) || !iszero(input) && iszero(result)

test/runtests.jl

@@ -42,6 +42,19 @@ using Dates:
 
 const colon = Base.:(:)
 
+macro test_or_throws(extype, ex)
+    return :(
+        try
+            # if the first line throws, go to @test_throws in catch clause
+            # if not: test expression normally
+            result = $ex
+            @test result
+        catch
+            @test_throws $extype $ex
+        end
+    )
+end
+
 @testset "Construction" begin
     @test isa(NoUnits, FreeUnits)
     @test typeof(𝐋) === Unitful.Dimensions{(Unitful.Dimension{:Length}(1),)}
@@ -199,10 +212,25 @@ end
             @test 1u"rps" == 2π/s
             @test 1u"rpm" == 360°/minute
             @test 1u"rpm" == 2π/minute
-
             # Issue 458:
             @test deg2rad(360°) ≈ 2π * rad
             @test rad2deg(2π * rad) ≈ 360°
+            # Issue 647:
+            @test uconvert(u"kb^1000", 1u"kb^1001 * b^-1") === 1000u"kb^1000"
+            @test uconvert(u"kOe^1000", 1u"kOe^1001 * Oe^-1") === 1000u"kOe^1000"
+            # Floating point overflow/underflow in uconvert can happen if the
+            # conversion factor is large, because uconvert does not cancel
+            # common basefactors (or just for really large exponents and/or
+            # SI prefixes). This test makes sure that uconvert does not silently
+            # return NaN, Inf, or 0 in these cases, i.e. either returns a finite
+            # result or throws an error indicating that it cannot handle the
+            # conversion.
+            is_finite_nonzero(x) = isfinite(x) && !iszero(x)
+            @test_or_throws ArgumentError is_finite_nonzero(uconvert(u"kb^12", 1u"b^12"))
+            @test_or_throws ArgumentError is_finite_nonzero(uconvert(u"ab^11", 1u"Tb^11"))
+            @test_or_throws ArgumentError is_finite_nonzero(uconvert(u"Tb^11", 1u"ab^11"))
+            @test_or_throws ArgumentError is_finite_nonzero(uconvert(u"b^11 * eV", 1u"m^22 * J"))
+            @test_or_throws ArgumentError is_finite_nonzero(uconvert(u"m^22 * J", 1u"b^11 * eV"))
         end
     end
 end