diff --git a/src/conversion.jl b/src/conversion.jl index a1a8bddf..8a32ac63 100644 --- a/src/conversion.jl +++ b/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 """ diff --git a/src/units.jl b/src/units.jl index f6a3f89a..89df1bf4 100644 --- a/src/units.jl +++ b/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) diff --git a/src/utils.jl b/src/utils.jl index e57df1b7..2dd03169 100644 --- a/src/utils.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index 23a42fd4..f85029bd 100644 --- a/test/runtests.jl +++ b/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