Reworking ComplexInfinity

src/Infinities.jl

@@ -1,20 +1,28 @@
 module Infinities
 
-import Base: angle, isone, iszero, isinf, isfinite, abs, one, oneunit, zero, isless, inv,
-                +, -, *, ^, ==, <, ≤, >, ≥, fld, cld, div, mod, min, max, sign, signbit,
-                string, show, promote_rule, convert, getindex, tryparse,
+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,
+                isinteger, round, floor, ceil, trunc, float,
                 Bool, Integer
 
 export ∞,  ℵ₀,  ℵ₁, RealInfinity, ComplexInfinity, InfiniteCardinal, NotANumber, PositiveInfinity, NegativeInfinity
 # The following is commented out for now to avoid conflicts with Infinity.jl
 # export Infinity
 
 """
-NotANumber()
+    NotANumber()
 
-represents something that is undefined, for example, `0 * ∞`.
+Construct the undefined value, for example the result of `0 * ∞`.
+
+Every float type has a `NaN` of its own. This one belongs to none of them.
 """
-struct NotANumber <: Number end
+struct NotANumber <: Real end
+
+(::Type{T})(::NotANumber) where {T<:AbstractFloat} = T(NaN)
+float(::NotANumber) = NaN
+Base.hash(::NotANumber, h::UInt)::UInt = hash(NaN, h)
 
 
 """
@@ -88,47 +96,90 @@ 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*∞            # cispi(0.5)∞
+    (1+im)*∞        # cispi(0.25)∞
+    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}()
+# A full turn fills the `UInt64` range, so a half turn is 2^63 units.
+const _HALFTURN = UInt64(2)^63 # the negative real axis
+# `_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(signbit(x))) # 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{Bool}) = y.signbit
-signbit(y::ComplexInfinity{<:Integer}) = !(mod(y.signbit,2) == 0)
-signbit(y::ComplexInfinity) = y.signbit
+# `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
-angle(x::ComplexInfinity) = π*x.signbit
+sign(y::ComplexInfinity) = cispi(_halfturns(y))
+angle(x::ComplexInfinity) = _halfturns(x) * π
+abs(::ComplexInfinity) = ∞
+conj(y::ComplexInfinity) = ComplexInfinity(-y.turns)
 
-show(io::IO, x::ComplexInfinity) = print(io, "exp($(x.signbit)*im*π)∞")
+# 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(_halfturns(x))
+    complex(_ray(c), _ray(s))
+end
+
+# 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)
+    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
@@ -138,12 +189,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
 
 

src/algebra.jl

@@ -6,28 +6,36 @@
 @inline infpromote(x::RealInfinity, y::Union{Integer, Rational}) = (x, float(y))
 @inline infpromote(x::Union{Integer, Rational}, y::RealInfinity) = (float(x), y)
 @inline infpromote(x::RealInfinity, ::InfiniteCardinal) = (x, ∞)
+# `Base` promotes every `Real` to `BigFloat`, which would convert the infinity away.
+@inline infpromote(x::BigFloat, y::Union{Infinity,RealInfinity}) = (x, y)
+@inline infpromote(x::Union{Infinity,RealInfinity}, y::BigFloat) = (x, y)
 
 
 # sign
 +(::Infinity) = RealInfinity()
 -(::Infinity) = RealInfinity(true)
 -(y::RealInfinity) = RealInfinity(!signbit(y))
--(y::ComplexInfinity{B}) where B<:Integer = sign(y) == 1 ? ComplexInfinity(one(B)) : ComplexInfinity(zero(B))
+-(y::ComplexInfinity) = ComplexInfinity(y.turns ⊻ _HALFTURN)
 +(x::InfiniteCardinal) = x
 -(::InfiniteCardinal) = -∞
 
 
 # addition
 @inline toinf(x) = RealInfinity(signbit(x))
-@inline toinf(x::Complex) = ComplexInfinity(angle(x))
+@inline toinf(x::Complex) = ComplexInfinity(_directionof(x))
 @inline toinf(x::ComplexInfinity) = x
 
-@inline _infadd(x, y) = angle(x) == angle(y) ? y : throw(ArgumentError("Angles must be the same to add ∞"))
+# The undefined value that matches the operands, as `Base` returns `NaN` or `NaN + NaN*im`.
+@inline _undefined(x, y) =
+    x isa ExtendedComplex || y isa ExtendedComplex ? ComplexNotANumber : NotANumber()
+
+@inline _infadd(x, y) = angle(x) == angle(y) ? y : _undefined(x, y)
 
 @inline __add(x, y::AllInfinities) = isinf(x) ? _infadd(toinf(x), y) : y
 @inline __add(x::Integer, y::InfiniteCardinal) = max(x, y)
 
-@inline _add(x, y) = __add(infpromote(x, y)...)
+# A `NaN` argument makes the result undefined. Types with no `NaN` fold the test away.
+@inline _add(x, y) = isnan(x) ? _undefined(x, y) : __add(infpromote(x, y)...)
 
 +(x::Number, y::AllInfinities) = _add(x, y)
 +(x::AllInfinities, y::Number) = _add(y, x)
@@ -44,16 +52,21 @@
 
 # multiplication
 
-@inline _sb(x) = signbit(x)
-@inline _sb(x::Complex) = angle(x)/π # overloading `signbit` causes type piracy 
+# The count of the direction a value points in. `_turns` instead reads its argument as a
+# number of half turns.
+@inline _directionof(x::Real) = signbit(x) ? _HALFTURN : zero(UInt64)
+@inline _directionof(x::Complex) = _turns(angle(x) / π) # overloading `signbit` causes type piracy
+@inline _directionof(x::ComplexInfinity) = x.turns
 
-@inline __mul(x, y::AllInfinities) = RealInfinity(_sb(x) ⊻ _sb(y))
-@inline __mul(x, y::ComplexInfinity) = ComplexInfinity(_sb(x) + _sb(y))
-@inline __mul(x, y::ComplexInfinity{Bool}) = ComplexInfinity(_sb(x) ⊻ _sb(y))
-@inline __mul(x::Complex, y::ComplexInfinity{Bool}) = ComplexInfinity(_sb(x) + _sb(y))
+@inline __mul(x, y::AllInfinities) = RealInfinity(signbit(x) ⊻ signbit(y))
+@inline __mul(x, y::ComplexInfinity) = ComplexInfinity(_directionof(x) + _directionof(y))
 @inline __mul(x::Integer, y::InfiniteCardinal) = x > 0 ? y : throw(ArgumentError("Cannot multiply $x * $y"))
 
-@inline _mul(x, y) = iszero(x) ? throw(ArgumentError("Cannot multiply $x * $y")) : __mul(infpromote(x, y)...)
+@inline function _mul(x, y)
+    isnan(x) && return _undefined(x, y)
+    iszero(x) && return _undefined(x, y)
+    __mul(infpromote(x, y)...)
+end
 
 *(x::Number, y::AllInfinities) = _mul(x, y)
 *(x::AllInfinities, y::Number) = _mul(y, x)
@@ -66,24 +79,43 @@
 *(::Infinity, ::Infinity) = ∞
 
 
+# division
+# `\` needs nothing of its own, `Base` defining it as `y / x`.
+@inline _div(x, y) = x * inv(y)
+
+/(x::AllInfinities, y::Number) = _div(x, y)
+/(x::Number, y::AllInfinities) = _div(x, y)
+/(x::AllInfinities, y::AllInfinities) = _undefined(x, y)
+
 # mod
 @inline function _mod(x::Real, y::IntegerInfinities)
+    isnan(x) && return NotANumber()
     signbit(x) == signbit(y) || throw(ArgumentError("mod($x,$y) is unbounded"))
     x
 end
 mod(x::Real, y::IntegerInfinities) = _mod(x, y)
 mod(::IntegerInfinities, ::Real) = NotANumber()
 mod(::IntegerInfinities, ::IntegerInfinities) = NotANumber()
 
+# rem, divrem
+# `rem` keeps the sign of the dividend, so unlike `mod` it stays bounded either way.
+rem(x::Real, ::IntegerInfinities) = isnan(x) ? NotANumber() : x
+rem(::IntegerInfinities, ::Real) = NotANumber()
+rem(::IntegerInfinities, ::IntegerInfinities) = NotANumber()
+# `Base` computes the remainder of two `Integer`s as `a - div(a,b)*b`, which an `InfiniteCardinal` cannot evaluate.
+divrem(x::Real, y::IntegerInfinities) = (div(x, y), rem(x, y))
+divrem(x::IntegerInfinities, y::Real) = (div(x, y), rem(x, y))
+divrem(x::IntegerInfinities, y::IntegerInfinities) = (div(x, y), rem(x, y))
+
 # fld, cld, div
-_divinf(T) = zero(T)
-_fldinf(x) = signbit(x) ? -one(x) : zero(x)
-_cldinf(x) = signbit(x) ? zero(x) : one(x)
-div(::T, ::IntegerInfinities) where T <: Real = _divinf(T)
+_divinf(x) = isnan(x) ? NotANumber() : zero(x)
+_fldinf(x) = isnan(x) ? NotANumber() : signbit(x) ? -one(x) : zero(x)
+_cldinf(x) = isnan(x) ? NotANumber() : signbit(x) ? zero(x) : one(x)
+div(x::Real, ::IntegerInfinities) = _divinf(x)
 fld(x::Real, ::IntegerInfinities) = _fldinf(x)
 cld(x::Real, ::IntegerInfinities) = _cldinf(x)
 
-_inffcd(x, y) = signbit(y) ? -x : x
+_inffcd(x, y) = isnan(y) ? NotANumber() : signbit(y) ? -x : x
 for OP in (:fld,:cld,:div)
     @eval begin
         $OP(x::IntegerInfinities, y::Real) = _inffcd(x, y)
@@ -94,8 +126,9 @@ end
 # power
 # Although the base implementation can cover these cases, it can change overtime and yield inconsistent results.
 # ref: https://github.com/JuliaMath/Infinities.jl/actions/runs/19993302836/
-_infpow(::PositiveInfinity, p) = ifelse(iszero(p), one(p), ifelse(p > 0, +∞, +zero(p)))
+_infpow(::PositiveInfinity, p) = isnan(p) ? NotANumber() : ifelse(iszero(p), one(p), ifelse(p > 0, +∞, +zero(p)))
 function _infpow(x::NegativeInfinity, p)
+    isnan(p) && return NotANumber()
     !isinteger(p) && throw(Base.Math.throw_exp_domainerror(x))
     iszero(p) && return one(p)
     isodd(p) && return ifelse(p > 0, -∞, -zero(p))
@@ -107,4 +140,36 @@ end
 # inv
 inv(::Union{Infinity,InfiniteCardinal}) = 0
 inv(x::RealInfinity) = inv(float(x))
-inv(x::ComplexInfinity) = zero(ComplexF64)
+inv(x::ComplexInfinity) = zero(ComplexF64)
+
+
+# NotANumber
+# Anything computed from an undefined value is undefined again, as it is for `NaN`.
+for op in (:+, :-, :*, :/, :^, :div, :fld, :cld, :mod, :rem, :min, :max)
+    for Typ in NotANumberRivals
+        @eval $op(x::NotANumber, ::$Typ) = x
+        @eval $op(::$Typ, y::NotANumber) = y
+    end
+    for Typ in NotANumberComplexRivals
+        @eval $op(::NotANumber, ::$Typ) = ComplexNotANumber
+        @eval $op(::$Typ, ::NotANumber) = ComplexNotANumber
+    end
+    @eval $op(x::NotANumber, ::NotANumber) = x
+end
+for Typ in NotANumberRivals
+    @eval divrem(x::NotANumber, ::$Typ) = (x, x)
+    @eval divrem(::$Typ, y::NotANumber) = (y, y)
+end
+for Typ in NotANumberComplexRivals
+    @eval divrem(::NotANumber, ::$Typ) = (ComplexNotANumber, ComplexNotANumber)
+    @eval divrem(::$Typ, ::NotANumber) = (ComplexNotANumber, ComplexNotANumber)
+end
+divrem(x::NotANumber, ::NotANumber) = (x, x)
+# `Base` has its own `^(::Number, ::Integer)`, which a literal exponent also routes through.
+^(x::NotANumber, ::Integer) = x
+^(::Integer, y::NotANumber) = y
+^(x::NotANumber, ::Rational) = x
+^(::Irrational{:ℯ}, y::NotANumber) = y
+for f in (:+, :-, :abs, :inv, :sign, :conj)
+    @eval $f(x::NotANumber) = x
+end

src/ambiguities.jl

@@ -1,10 +1,9 @@
 for Typ in (Base.TwicePrecision, AbstractChar, Complex)
-    @eval begin
-        RealInfinity(x::$Typ) = throw(MethodError(RealInfinity, x))
-        ComplexInfinity{T}(x::$Typ) where T<:Real = ComplexInfinity(T(x))
-    end
+    @eval RealInfinity(x::$Typ) = throw(MethodError(RealInfinity, x))
 end
-ComplexInfinity{T}(x::ComplexInfinity{T}) where T<:Real = x
+# `Base` builds a `TwicePrecision` of any type by adding its two halves, which for two opposed
+# directions is undefined.
+ComplexInfinity(x::Base.TwicePrecision) = ComplexInfinity(halfturns = Float64(x))
 
 for Typ in (Rational, BigInt, BigFloat)
     for (op, fop) in ((:<, :_lt), (:≤, :_le))
@@ -19,23 +18,35 @@ for Typ in (Rational, BigInt, BigFloat, Complex, AbstractIrrational)
 end
 
 for Typ in (Complex, Rational, Complex{Bool}, Integer)
-    @eval +(x::AllInfinities, y::$Typ) = _add(y, x)
-    @eval +(x::$Typ, y::AllInfinities) = _add(x, y)
-    @eval -(x::AllInfinities, y::$Typ) = _sub(x, y)
-    @eval -(x::$Typ, y::AllInfinities) = _sub(x, y)
-    @eval *(x::AllInfinities, y::$Typ) = _mul(y, x)
-    @eval *(x::$Typ, y::AllInfinities) = _mul(x, y)
+    # `_add` and `_mul` dispatch on the infinity being second; `_sub` and `_div` delegate to them.
+    for (op, fop) in ((:+, :_add), (:*, :_mul))
+        @eval $op(x::AllInfinities, y::$Typ) = $fop(y, x)
+        @eval $op(x::$Typ, y::AllInfinities) = $fop(x, y)
+    end
+    for (op, fop) in ((:-, :_sub), (:/, :_div))
+        @eval $op(x::AllInfinities, y::$Typ) = $fop(x, y)
+        @eval $op(x::$Typ, y::AllInfinities) = $fop(x, y)
+    end
 end
 
 ^(x::RealInfinity, y::Rational) = _infpow(infpromote(x, y)...)
 
 for Typ in (Rational, )
     @eval mod(::IntegerInfinities, ::$Typ) = NotANumber()
     @eval mod(x::$Typ, y::IntegerInfinities) = _mod(x, y)
+    @eval rem(::InfiniteCardinal, ::$Typ) = NotANumber()
+    @eval rem(x::$Typ, ::IntegerInfinities) = x
     for op in (:fld, :cld, :div)
         @eval $op(x::InfiniteCardinal, y::$Typ) = _inffcd(x, y)
     end
-    @eval div(::T, ::IntegerInfinities) where T <: $Typ = _divinf(T)
+    @eval div(x::$Typ, ::IntegerInfinities) = _divinf(x)
     @eval fld(x::$Typ, ::IntegerInfinities) = _fldinf(x)
     @eval cld(x::$Typ, ::IntegerInfinities) = _cldinf(x)
-end
+end
+
+divrem(x::BigInt, y::IntegerInfinities) = (div(x, y), rem(x, y))
+
+# an `InfiniteCardinal` is an `Integer`, for which `Base` has its own `isapprox`
+isapprox(x::InfiniteCardinal, y::Integer; kwargs...) = x == y
+isapprox(x::Integer, y::InfiniteCardinal; kwargs...) = x == y
+isapprox(x::InfiniteCardinal, y::InfiniteCardinal; kwargs...) = x == y