Add a ForwardDiff.jl extension

Project.toml

@@ -3,18 +3,30 @@ uuid = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 authors = ["Sheehan Olver <solver@mac.com>"]
 version = "0.1.13"
 
+[weakdeps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+
+[extensions]
+InfinitiesForwardDiffExt = "ForwardDiff"
+InfinitiesStaticExt = "Static"
+
 [compat]
 Aqua = "0.8"
 Base64 = "1"
+ForwardDiff = "1"
 JET = "0.9, 0.10, 0.11, 0.12"
+Static = "1.4"
 Test = "1"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Base64", "JET", "Test"]
+test = ["Aqua", "Base64", "ForwardDiff", "JET", "Static", "Test"]

README.md

@@ -16,6 +16,59 @@ This Julia package is used to represent infinities, including:
 
 Note that we subtype based on interfaces, rather than strict mathematical definitions. For example,  `ℵ₀ isa Integer` as `Integer` is often used to represent the size of a set or vector. Similarly, `∞ isa Real`.
 
+## Static.jl integration
+
+Loading [Static.jl](https://github.com/SciML/Static.jl) enables an optional package extension.
+The values `∞`, `+∞`, `-∞`, `InfiniteCardinal{k}()`, and `NotANumber()` are already
+fully determined by their types, so they need no separate static representation:
+
+```julia
+using Infinities, Static
+
+static(∞) === ∞                       # true
+is_static(typeof(-∞)) === True()      # true
+known(typeof(ℵ₀)) === ℵ₀               # true
+Static.lt(static(2), ∞) === True()     # true
+static(2) + ℵ₀ === ℵ₀                 # true
+```
+
+Mixed arithmetic and comparisons use the same rules as the corresponding ordinary
+numbers, including undefined results and NaN propagation. Base operations retain their
+ordinary return types: `one(∞)` is `1` and `isinf(∞)` is `true`. Use `static` on a result
+or Static.jl's comparison functions when a static result type is needed.
+
+`ComplexInfinity` stores its direction as a value, not in its type. Its arithmetic
+accepts static operands, but `is_static(ComplexInfinity)` is `False()` and
+`static(im*∞)` is unsupported.
+
+## ForwardDiff.jl integration
+
+Loading [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) enables an optional
+extension for mixed `+`, `-`, `*`, `/`, `mod`, `rem`, `min`, `max`, and comparisons
+of dual numbers with `∞`, `+∞`, `-∞`, and `NotANumber()`. Primal values retain
+Infinities' scalar results and exceptions, without converting infinities to floats.
+Derivative rules use the same scalar arithmetic, so a zero tangent multiplied by an
+infinity becomes `NotANumber()`. Bounded `mod` and `rem` preserve the dividend's
+tangents, while an undefined remainder has undefined tangents.
+
+```julia
+using Infinities, ForwardDiff
+
+ForwardDiff.derivative(input -> input * ∞, 2.0) # +∞
+ForwardDiff.derivative(input -> input + ∞, 2.0) # 1.0
+ForwardDiff.derivative(input -> input / ∞, 2.0) # 0.0
+```
+
+`Dual(∞)` preserves the infinity. Explicitly requesting a floating scalar type, such
+as `Dual{Nothing, Float64}(∞)`, converts it. Mixed symbolic results can use `Real`
+as the dual's scalar parameter, with a corresponding loss of type specialization.
+
+This is not general support for arbitrary symbolic differentiation: ForwardDiff
+operations such as unary negation can reject heterogeneous symbolic partials, and
+powers remain subject to the existing scalar and ForwardDiff limitations. There is
+no floating-point fallback. The extension does not define mixed dual-number operations
+for `InfiniteCardinal` or `ComplexInfinity`.
+
 ## Similar packages
 
 This package is meant to eventually replace [Infinity.jl](https://github.com/cjdoris/Infinity.jl) and the definitions of `∞` in [InfiniteArrays.jl](https://github.com/JuliaArrays/InfiniteArrays.jl). We do not yet support Infinity.jl's notions of `InfExtendedReal` but we hope to add this soon.

ext/InfinitiesForwardDiffExt.jl

@@ -0,0 +1,51 @@
+module InfinitiesForwardDiffExt
+
+using Infinities: Infinity, RealInfinity, NotANumber
+import ForwardDiff: Dual, Partials, value, partials
+
+function symbolic_dual(::Dual{Tag}, primal, tangents::NTuple{N, Any}) where {Tag, N}
+    Dual{Tag, Real, N}(primal, Partials{N, Real}(tangents))
+end
+
+for Typ in (Infinity, RealInfinity, NotANumber)
+    for op in (:^, :min, :max, :(==), :isequal, :<, :<=, :isless)
+        @eval Base.$op(dual::Dual, inf::$Typ) = invoke($op, Tuple{Dual, Real}, dual, inf)
+        @eval Base.$op(inf::$Typ, dual::Dual) = invoke($op, Tuple{Real, Dual}, inf, dual)
+    end
+
+    for (op, forward, reverse) in ((:+, :tangent, :tangent),
+        (:-, :tangent, :(-tangent)), (:*, :(tangent * inf), :(tangent * inf)),
+        (:/, :(tangent / inf), :(-(primal / value(dual)) * tangent)))
+        @eval function Base.$op(dual::Dual, inf::$Typ)
+            primal = $op(value(dual), inf)
+            symbolic_dual(dual, primal, map(tangent -> $forward, Tuple(partials(dual))))
+        end
+        @eval function Base.$op(inf::$Typ, dual::Dual)
+            primal = $op(inf, value(dual))
+            symbolic_dual(dual, primal, map(tangent -> $reverse, Tuple(partials(dual))))
+        end
+    end
+
+    for op in (:mod, :rem)
+        @eval function Base.$op(dual::Dual, inf::$Typ)
+            primal = $op(value(dual), inf)
+            symbolic_dual(dual, primal, map(tangent -> isnan(primal) ? NotANumber() : tangent, Tuple(partials(dual))))
+        end
+        @eval Base.$op(inf::$Typ, dual::Dual) =
+            symbolic_dual(dual, $op(inf, value(dual)), map(_ -> NotANumber(), Tuple(partials(dual))))
+    end
+
+    @eval begin
+        Dual(inf::$Typ) = Dual{Nothing}(inf, ())
+        Dual{Tag}(inf::$Typ) where {Tag} = Dual{Tag}(inf, ())
+        Dual{Tag, Value}(inf::$Typ) where {Tag, Value} = Dual{Tag, Value, 0}(inf)
+        Dual{Tag, Value, N}(inf::$Typ) where {Tag, Value, N} =
+            Dual{Tag, Value, N}(convert(Value, inf), zero(Partials{N, Value}))
+        Dual{Tag}(inf::$Typ, tangents::Partials{N, Value}) where {Tag, N, Value} =
+            Dual{Tag, Real, N}(inf, Partials{N, Real}(Tuple(tangents)))
+        Dual{Tag}(inf::Value, tangents::Partials{N, Value}) where {Tag, N, Value<:$Typ} =
+            Dual{Tag, Real, N}(inf, Partials{N, Real}(Tuple(tangents)))
+    end
+end
+
+end

ext/InfinitiesStaticExt.jl

@@ -0,0 +1,37 @@
+module InfinitiesStaticExt
+
+using Infinities: Infinity, PositiveInfinity, NegativeInfinity, InfiniteCardinal, NotANumber
+using Infinities: AllInfinities, AllRealInfinities, IntegerInfinities, ComplexInfinity, RealInfinity
+using Static: Static, dynamic, StaticNumber, StaticInteger, StaticFloat64, StaticInt, True
+
+for Typ in (Infinity, PositiveInfinity, NegativeInfinity, NotANumber)
+    @eval begin
+        Static.static(x::$Typ) = x
+        Static.is_static(::Type{$Typ}) = True()
+        Static.known(::Type{$Typ}) = $Typ()
+    end
+end
+
+Static.static(x::InfiniteCardinal) = x
+Static.is_static(::Type{InfiniteCardinal{N}}) where {N} = True()
+Static.known(::Type{InfiniteCardinal{N}}) where {N} = InfiniteCardinal{N}()
+
+for (ops, Types, StaticTypes) in (
+    ((:+, :-, :*, :/, :(==)), (AllInfinities,), (StaticNumber,)),
+    ((:isequal,), (NotANumber,), (StaticNumber,)),
+    ((:div, :fld, :cld, :divrem), (IntegerInfinities,), (StaticNumber,)),
+    ((:<, :<=, :>, :>=), (AllInfinities,), (StaticInteger, StaticFloat64)),
+    ((:isless,), (AllRealInfinities, InfiniteCardinal, NotANumber), (StaticInteger, StaticFloat64)),
+    ((:min, :max), (AllInfinities, NotANumber), (StaticInteger,)),
+    ((:min, :max), (IntegerInfinities, ComplexInfinity, NotANumber), (StaticFloat64,)),
+    ((:mod,), (IntegerInfinities, NotANumber), (StaticNumber,)),
+    ((:rem,), (IntegerInfinities, NotANumber), (StaticInteger, StaticFloat64)),
+    ((:*,), (InfiniteCardinal,), (StaticInt{0},)),
+), op in ops, Typ in Types, StaticTyp in StaticTypes
+    @eval Base.$op(x::$Typ, y::$StaticTyp) = $op(x, dynamic(y))
+    @eval Base.$op(x::$StaticTyp, y::$Typ) = $op(dynamic(x), y)
+end
+
+Base.:^(x::RealInfinity, y::StaticNumber) = x^dynamic(y)
+
+end

src/Infinities.jl

@@ -3,7 +3,7 @@ module Infinities
 import Base: angle, isone, iszero, isinf, isfinite, isnan, isreal, abs, one, oneunit, zero, isless, isequal, inv,
                 +, -, *, /, ^, ==, <, ≤, >, ≥, fld, cld, div, mod, rem, divrem, min, max,
                 sign, signbit, isapprox,
-                string, show, promote_rule, convert, getindex, tryparse, conj,
+                string, show, promote_rule, convert, getindex, tryparse, conj, complex,
                 isinteger, round, floor, ceil, trunc, float,
                 Bool, Integer
 
@@ -44,7 +44,7 @@ _convert(::Type{T}, ::Infinity) where {T<:Real} = convert(T, Inf)::T
 (::Type{T})(x::Infinity) where {T<:Real} = _convert(T, x)
 
 sign(y::Infinity) = 1
-angle(x::Infinity) = 0
+angle(x::Infinity) = 0.0
 signbit(::Infinity) = false
 
 one(::Type{Infinity}) = 1
@@ -96,57 +96,99 @@ zero(::Type{RealInfinity}) = 0.0
 # ComplexInfinity
 #######
 
-# angle is π*a where a is (false==0) and (true==1)
-
 """
-ComplexInfinity(signbit)
+    ComplexInfinity(turns::UInt64)
+    ComplexInfinity(; halfturns::Real = 0)
+
+Construct an infinity in the complex plane, pointing in a direction held as a count of
+`2^-64` turns.
+
+The count wraps at a full turn, so the stored `UInt64` and the directions are bijective.
+`0x0` points along the positive real axis. Values increase counterclockwise.
+`0x8000000000000000` points along the negative real axis. Use `reinterpret(UInt64, x)` to
+read out the exact value.
+
+Multiplying by `∞` takes the direction from the other operand, which usually reads better
+than naming an angle:
 
-represents an infinity in the complex plane with the angle
-specified by `π * signbit`. The use of the name `signbit` is
-for consistency with `RealInfinity`.
+    im*∞            # 0 + ∞*im
+    (1+im)*∞        # ∞ + ∞*im
+    exp(im*π/4)*∞   # the same direction again
+
+Those forms and the `halfturns` keyword go through `angle`, so they round. It is exact on
+the axes and at a quarter turn, but `exp(im*π/8)*∞` lands 256 counts past an eighth turn.
+Provide the `UInt64` when you have an off-axis value where accuracy matters.
 """
-struct ComplexInfinity{T<:Real} <: Number
-    signbit::T
+struct ComplexInfinity <: Number
+    turns::UInt64
+    ComplexInfinity(turns::UInt64) = new(turns)
 end
 
-ComplexInfinity{T}() where T = ComplexInfinity(zero(T))
-ComplexInfinity() = ComplexInfinity{Bool}()
-ComplexInfinity{T}(::Infinity) where T<:Real = ComplexInfinity{T}()
+# Half of the `typemax(UInt64) + 1` counts of a full turn, a number that would overflow.
+const _HALFTURN = typemax(UInt64) ÷ 2 + 1 # the negative real axis
+# The signs of the parts along each of the eight rays, the axes and diagonals.
+const _EIGHTH = _HALFTURN ÷ 4
+const _RAYPARTS = ((1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1))
+# `_turns` and `_halfturns` are inverse: half turns in, count out, and back again.
+# `mod` returns 2 itself for a tiny negative angle, since 2 + x rounds back to 2. A full
+# turn is the direction zero. Testing `== 2` rather than `< 2` still lets `NaN` throw.
+@inline _turns(halfturns::Real) = round(UInt64, (h = mod(halfturns, 2); h == 2 ? zero(h) : h) * 0x1p63)
+# Scaling by 2^63 needs 63 bits beyond the numerator, so `Int128` has room for any `Int64`.
+_turns(halfturns::Rational) = round(BigInt, mod(halfturns, 2) * big(2)^63) % UInt64
+_turns(halfturns::Rational{<:Base.BitInteger64}) =
+    round(Int128, mod(halfturns, 2) * Int128(2)^63) % UInt64
+# `Base` puts an angle in `(-π, π]`, so past the half turn the count reads as negative.
+@inline _halfturns(x::ComplexInfinity) =
+    x.turns == _HALFTURN ? 1.0 : reinterpret(Int64, x.turns) / 0x1p63
+
+ComplexInfinity(; halfturns::Real = 0) = ComplexInfinity(_turns(halfturns))
 ComplexInfinity(::Infinity) = ComplexInfinity()
-ComplexInfinity{T}(x::RealInfinity) where T<:Real = ComplexInfinity{T}(signbit(x))
-ComplexInfinity(x::RealInfinity) = ComplexInfinity(signbit(x))
-ComplexInfinity{T}(x::ComplexInfinity) where T<:Real = ComplexInfinity(T(x.signbit)) # ambiguity fix
+ComplexInfinity(x::RealInfinity) = ComplexInfinity(_directionof(x))
+ComplexInfinity(x::ComplexInfinity) = x
+
+signbit(y::ComplexInfinity) = y.turns == _HALFTURN
+isreal(y::ComplexInfinity) = iszero(y.turns) || signbit(y)
 
-signbit(y::ComplexInfinity) = mod(y.signbit, 2) == 1
+# `Base` converts a `Complex` to a `Real` the same way, and throws the same error off the axis.
+RealInfinity(x::ComplexInfinity) = isreal(x) ? RealInfinity(signbit(x)) :
+                                   throw(InexactError(:RealInfinity, RealInfinity, x))
 
-convert(::Type{ComplexInfinity{T}}, ::Infinity) where T = ComplexInfinity{T}()
 convert(::Type{ComplexInfinity}, ::Infinity) = ComplexInfinity()
-convert(::Type{ComplexInfinity{T}}, x::RealInfinity) where T = ComplexInfinity{T}(x)
 convert(::Type{ComplexInfinity}, x::RealInfinity) = ComplexInfinity(x)
 
 
-sign(y::ComplexInfinity{<:Integer}) = mod(y.signbit, 2) == 0 ? 1 : -1
-sign(y::ComplexInfinity) = cispi(y.signbit)
-angle(x::ComplexInfinity) = π*x.signbit
+sign(y::ComplexInfinity) = cispi(_halfturns(y))
+angle(x::ComplexInfinity) = _halfturns(x) * π
 abs(::ComplexInfinity) = ∞
-conj(y::ComplexInfinity{<:Integer}) = y # an integer factor points along the real axis
-conj(y::ComplexInfinity) = ComplexInfinity(mod(-y.signbit, 2))
+conj(y::ComplexInfinity) = ComplexInfinity(-y.turns)
 
 # An exact zero has to stay finite, `Inf * 0` being a `NaN`.
 @inline _ray(c) = iszero(c) ? c : copysign(Inf, c)
 # `Complex` reaches only the eight rays of its two saturating parts, so the direction lands on the nearest of them.
 function float(x::ComplexInfinity)
-    s, c = sincospi(x.signbit)
+    s, c = sincospi(_halfturns(x))
     complex(_ray(c), _ray(s))
 end
 
-show(io::IO, x::ComplexInfinity) = print(io, "exp($(x.signbit)*im*π)∞")
+# The rays print as `Base` prints an infinite `Complex`. Off the rays the readable form
+# names an angle, which recovers most counts but not all, so it is used only where reading
+# it back gives the same direction.
+function show(io::IO, x::ComplexInfinity)
+    k, offset = divrem(x.turns, _EIGHTH)
+    if iszero(offset)
+        r, i = _RAYPARTS[k + 1]
+        return print(io, ("-∞", "0", "∞")[r + 2], (" - ∞*im", " + 0im", " + ∞*im")[i + 2])
+    end
+    h = _halfturns(x)
+    _directionof(cispi(h)) == x.turns ? print(io, "cispi($h)∞") :
+                                    print(io, "ComplexInfinity(", repr(x.turns), ")")
+end
 
-one(::Type{<:ComplexInfinity}) = one(ComplexF64)
-oneunit(::Type{<:ComplexInfinity}) = oneunit(ComplexF64)
+one(::Type{ComplexInfinity}) = one(ComplexF64)
+oneunit(::Type{ComplexInfinity}) = oneunit(ComplexF64)
 oneunit(::ComplexInfinity) = oneunit(ComplexF64)
 zero(::ComplexInfinity) = zero(ComplexF64)
-zero(::Type{<:ComplexInfinity}) = zero(ComplexF64)
+zero(::Type{ComplexInfinity}) = zero(ComplexF64)
 
 
 # `isequal` implies equal hashes, so the infinities have to hash like the float
@@ -156,12 +198,11 @@ Base.hash(::Infinity, h::UInt)::UInt = hash(Inf, h)
 Base.hash(::PositiveInfinity, h::UInt)::UInt = hash(Inf, h)
 Base.hash(::NegativeInfinity, h::UInt)::UInt = hash(-Inf, h)
 
-# Equality of ComplexInfinity is equality of the angle, hence so is the hash.
+# The two real directions have to hash like the real infinities they compare equal to.
 function Base.hash(x::ComplexInfinity, h::UInt)::UInt
-    θ = angle(x)
-    θ == angle(PositiveInfinity()) && return hash(Inf, h)
-    θ == angle(NegativeInfinity()) && return hash(-Inf, h)
-    hash(ComplexInfinity, hash(θ, h))
+    iszero(x.turns) && return hash(Inf, h)
+    x.turns == _HALFTURN && return hash(-Inf, h)
+    hash(ComplexInfinity, hash(x.turns, h))
 end