Re: [julia-users] Complex infinity
Julia's support for complex arithmetic is written purely in Julia, so I don't see how my suggestion to implement a new type is necessarily slower than hard-coding in special values indicating complex infinity. The code branches that one would need would simply be put in different places. Are there LLVM intrinsics and/or hardware support for complex arithmetic? The story would be different if there are (and I don't think we use them currently). The normalization of all complex infinities is essentially provided by the C9X proposal's cproj function (which uses Inf + 0.0im as the representation of ComplexInf). However, the C9X proposal acknowledges that there are use cases for both topologies and thus avoids recommending one or the other. To quote: Two topologies are commonly used in complex mathematics: the complex plane with its continuum of infinities and the Riemann sphere with its single infinity. The complex plane is better suited for transcendental functions, the Riemann sphere for algebraic functions. The complex types with their multiplicity of infinities provide a useful (though imperfect) model for the complex plane. The cproj function helps model the Riemann sphere by mapping all infinities to one, and should be used just before any operation, especially comparisons, that might give spurious results for any of the other infinities. In this proposal, the complex plane is clearly the more fundamental space. The Riemann sphere cannot be worked with directly; one can only work with the shadow of the Riemann sphere projected back onto the complex plane, one particular projection of which is provided by cproj. In contrast, it would be very Julian to represent the working topology with types rather than explicitly forcing users to choose between the ordinary complex plane and the Riemann sphere. The Riemann sphere would be a first-class working algebra instead of being accessible indirectly via cproj. The type model is also more composable, should there ever be a need to intermix the two topologies. The situation for zeros is less dire since they compare equal with ==, whereas the various representations of complex infinity do not. Thanks, Jiahao Chen Staff Research Scientist MIT Computer Science and Artificial Intelligence Laboratory
Re: [julia-users] Complex infinity
Instead of defining a new type, it should be faster to normalize complex
numbers after every operation. The current, IEEE-based complex numbers not only
have various representations of infinity, but also for zero, since they
distinguish between +0.0 and -0.0. Thus, the normalization would map all of the
four complex zeros to zero (e.g. +0.0+0.0im), and would also choose a single
normalized representation of infinity (maybe inf+inf*im, where inf=1.0/0.0).
I would probably be necessary to redefine most complex number arithmetic to
take these special cases into account.
If you do this for complex numbers, then you may also want to do it for real
numbers, just for good measure. Currently, there is +0.0 and -0.0 (and +inf and
-inf), and a real number corresponding to an extended complex number should
probably not distinguish between these either.
-erik
On Jan 10, 2015, at 14:37 , Jiahao Chen jia...@mit.edu wrote:
You might be interested in Issue #5234, where I tried to catalogue issues
arising from nonfinite floating-point computations in the ordinary complex
plane (without closure on the Riemann sphere).
Being able to switch between complex and extended complex (closed at
ComplexInf) would be interesting. I vaguely recall discussing complex
infinity with Mike Innes at some point, but I don't recall if we got anywhere.
One solution is to wrap ordinary complex numbers in a new type representing
finite extended complex numbers, and create yet another new type representing
complex infinity. Define each function over the exted For each function, and
redefine each function you want to use by either wrapping a finite answer in
the new type, or otherwise checking if the answer is infinite and returning
another new type instead which represents complex infinity.
abstract ExtendedComplex{T:Real}
immutable ComplexInf{T:Real} : ExtendedComplex{T} end
immutable FiniteExtendedComplex{T:Real} : ExtendedComplex{T}
z :: Complex{T}
end
function ExtendedComplex{T:Real}(r::T, i::T)
if isfinite(r) isfinite(i))
return ExtendedComplex(Complex(r, i))
else
return ComplexInf{T}()
end
end
import Base./
function /{T}(z:: ComplexInf{T}, w:: FiniteExtendedComplex{T})
return z
end
function /{T}(z:: FiniteExtendedComplex{T}, w:: ComplexInf{T})
return zero(z)
end
function /{T}(z:: FiniteExtendedComplex{T}, w:: FiniteExtendedComplex{T})
a = z.z/w.z
return isfinite(a) ? ExtendedComplex(a) : ComplexInf{T}
end
With these definitions you can compute as follows:
julia FiniteExtendedComplex(2.0, 0.5)/FiniteExtendedComplex(0.5, 2.0)
FiniteExtendedComplex{Float64}(0.47058823529411764 - 0.8823529411764706im)
julia FiniteExtendedComplex(2.0, 0.3)/FiniteExtendedComplex(0.0, 0.0)
ComplexInf{Float64}()
julia ComplexInf{Float64}()/FiniteExtendedComplex(0.0, 0.0)
ComplexInf{Float64}()
julia ComplexInf{Float64}()/ComplexInf{Float64}() #Don't know how to do ∞/∞
ERROR: `/` has no method matching /(::ComplexInf{Float64},
::ComplexInf{Float64})
Thanks,
Jiahao Chen
Staff Research Scientist
MIT Computer Science and Artificial Intelligence Laboratory
On Sat, Jan 10, 2015 at 8:55 AM, Ed Scheinerman
edward.scheiner...@gmail.com wrote:
Is there a way to have a single complex infinity? This may come at the cost
of computational efficiency I suppose, but I can think of situations where
all of the following give the same result:
julia (1+1im)/0
Inf + Inf*im
julia 1im/0
NaN + Inf*im
julia 1/0 + im
Inf + 1.0im
It would be nice (sometimes) if these were all the same ComplexInf, say.
Perhaps there's an extended complex numbers module for this sort of work?
--
Erik Schnetter schnet...@gmail.com
http://www.perimeterinstitute.ca/personal/eschnetter/
My email is as private as my paper mail. I therefore support encrypting
and signing email messages. Get my PGP key from http://pgp.mit.edu/.
signature.asc
Description: Message signed with OpenPGP using GPGMail
Re: [julia-users] Complex infinity
Whoops, that version didn't work. Try this version instead:
abstract ExtendedComplex{T:Real}
immutable ComplexInf{T:Real} : ExtendedComplex{T} end
immutable FiniteExtendedComplex{T:Real} : ExtendedComplex{T}
z :: Complex{T}
end
function FiniteExtendedComplex{T:Real}(r::T, i::T)
if isfinite(r) isfinite(i)
return FiniteExtendedComplex(Complex(r, i))
else
return ComplexInf{T}()
end
end
import Base./
function /{T}(z:: ComplexInf{T}, w:: FiniteExtendedComplex{T})
return z
end
function /{T}(z:: FiniteExtendedComplex{T}, w:: ComplexInf{T})
return zero(z)
end
function /{T}(z:: FiniteExtendedComplex{T}, w:: FiniteExtendedComplex{T})
a = z.z/w.z
return isfinite(a) ? FiniteExtendedComplex(a) : ComplexInf{T}()
end
Re: [julia-users] Complex infinity
You might be interested in Issue #5234
https://github.com/JuliaLang/julia/issues/5234, where I tried to
catalogue issues arising from nonfinite floating-point computations in the
ordinary complex plane (without closure on the Riemann sphere).
Being able to switch between complex and extended complex (closed at
ComplexInf) would be interesting. I vaguely recall discussing complex
infinity with Mike Innes at some point, but I don't recall if we got
anywhere.
One solution is to wrap ordinary complex numbers in a new type representing
finite extended complex numbers, and create yet another new type
representing complex infinity. Define each function over the exted For each
function, and redefine each function you want to use by either wrapping a
finite answer in the new type, or otherwise checking if the answer is
infinite and returning another new type instead which represents complex
infinity.
abstract ExtendedComplex{T:Real}
immutable ComplexInf{T:Real} : ExtendedComplex{T} end
immutable FiniteExtendedComplex{T:Real} : ExtendedComplex{T}
z :: Complex{T}
end
function ExtendedComplex{T:Real}(r::T, i::T)
if isfinite(r) isfinite(i))
return ExtendedComplex(Complex(r, i))
else
return ComplexInf{T}()
end
end
import Base./
function /{T}(z:: ComplexInf{T}, w:: FiniteExtendedComplex{T})
return z
end
function /{T}(z:: FiniteExtendedComplex{T}, w:: ComplexInf{T})
return zero(z)
end
function /{T}(z:: FiniteExtendedComplex{T}, w:: FiniteExtendedComplex{T})
a = z.z/w.z
return isfinite(a) ? ExtendedComplex(a) : ComplexInf{T}
end
With these definitions you can compute as follows:
julia FiniteExtendedComplex(2.0, 0.5)/FiniteExtendedComplex(0.5, 2.0)
FiniteExtendedComplex{Float64}(0.47058823529411764 - 0.8823529411764706im)
julia FiniteExtendedComplex(2.0, 0.3)/FiniteExtendedComplex(0.0, 0.0)
ComplexInf{Float64}()
julia ComplexInf{Float64}()/FiniteExtendedComplex(0.0, 0.0)
ComplexInf{Float64}()
julia ComplexInf{Float64}()/ComplexInf{Float64}() #Don't know how to do ∞/∞
ERROR: `/` has no method matching /(::ComplexInf{Float64},
::ComplexInf{Float64})
Thanks,
Jiahao Chen
Staff Research Scientist
MIT Computer Science and Artificial Intelligence Laboratory
On Sat, Jan 10, 2015 at 8:55 AM, Ed Scheinerman
edward.scheiner...@gmail.com wrote:
Is there a way to have a single complex infinity? This may come at the
cost of computational efficiency I suppose, but I can think of situations
where all of the following give the same result:
julia (1+1im)/0
Inf + Inf*im
julia 1im/0
NaN + Inf*im
julia 1/0 + im
Inf + 1.0im
It would be nice (sometimes) if these were all the same ComplexInf, say.
Perhaps there's an extended complex numbers module for this sort of work?
Re: [julia-users] Complex infinity
What I was imagining was a new type, called e.g. ExtendedComplex (or RiemannComplex) that substitutes for the current Complex type. I didn't know about cproj, but I was imagining something equivalent to RiemannComplex calls cproj after every complex number operation. If a function can return two different types, this is called type instability. It requires the compiler (at least currently, but I know of no plans to change this) to return values in heap-allocated data structures instead of in registers. This is very expensive. See e.g. http://julia.readthedocs.org/en/latest/manual/performance-tips/ for a discussion on type instability. Without this problem, using a separate type to denote infinity would indeed not only be a very elegant but also efficient solution. -erik On Jan 10, 2015, at 17:44 , Jiahao Chen jia...@mit.edu wrote: Julia's support for complex arithmetic is written purely in Julia, so I don't see how my suggestion to implement a new type is necessarily slower than hard-coding in special values indicating complex infinity. The code branches that one would need would simply be put in different places. Are there LLVM intrinsics and/or hardware support for complex arithmetic? The story would be different if there are (and I don't think we use them currently). The normalization of all complex infinities is essentially provided by the C9X proposal's cproj function (which uses Inf + 0.0im as the representation of ComplexInf). However, the C9X proposal acknowledges that there are use cases for both topologies and thus avoids recommending one or the other. To quote: Two topologies are commonly used in complex mathematics: the complex plane with its continuum of infinities and the Riemann sphere with its single infinity. The complex plane is better suited for transcendental functions, the Riemann sphere for algebraic functions. The complex types with their multiplicity of infinities provide a useful (though imperfect) model for the complex plane. The cproj function helps model the Riemann sphere by mapping all infinities to one, and should be used just before any operation, especially comparisons, that might give spurious results for any of the other infinities. In this proposal, the complex plane is clearly the more fundamental space. The Riemann sphere cannot be worked with directly; one can only work with the shadow of the Riemann sphere projected back onto the complex plane, one particular projection of which is provided by cproj. In contrast, it would be very Julian to represent the working topology with types rather than explicitly forcing users to choose between the ordinary complex plane and the Riemann sphere. The Riemann sphere would be a first-class working algebra instead of being accessible indirectly via cproj. The type model is also more composable, should there ever be a need to intermix the two topologies. The situation for zeros is less dire since they compare equal with ==, whereas the various representations of complex infinity do not. Thanks, Jiahao Chen Staff Research Scientist MIT Computer Science and Artificial Intelligence Laboratory -- Erik Schnetter schnet...@gmail.com http://www.perimeterinstitute.ca/personal/eschnetter/ My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://pgp.mit.edu/. signature.asc Description: Message signed with OpenPGP using GPGMail
Re: [julia-users] Complex infinity
This may lead to problems with nans. For example, 1.0/(0.0, 0.0) is (inf, nan), and 1.0/(inf,nan) is (0.0, nan); multiply by (1, im) yields (nan, nan), and all information is lost. With a proper treatment, this would be (0.0, 0.0) instead. -erik On Jan 10, 2015, at 18:07 , Jiahao Chen jia...@mit.edu wrote: Another possibility that occurs to me is to redefine the == equality comparison. Currently this is defined as ==(z::Complex, w::Complex) = (real(z) == real(w)) (imag(z) == imag(w)) but for some purposes it may be sufficient to redefine this as ==(z::Complex, w::Complex) = if isinf(z) isinf(w) return true else return (real(z) == real(w)) (imag(z) == imag(w)) end to model the topology of the Riemann sphere. With this redefinition: julia Complex(Inf, 0.0) == Complex(Inf, NaN) true julia Complex(Inf, 0.0) == Complex(Inf, 2.0) true -- Erik Schnetter schnet...@gmail.com http://www.perimeterinstitute.ca/personal/eschnetter/ My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://pgp.mit.edu/. signature.asc Description: Message signed with OpenPGP using GPGMail
Re: [julia-users] Complex infinity
On Sat, Jan 10, 2015 at 6:12 PM, Erik Schnetter schnet...@gmail.com wrote: This may lead to problems with nans. For example, 1.0/(0.0, 0.0) is (inf, nan), and 1.0/(inf,nan) is (0.0, nan); multiply by (1, im) yields (nan, nan), and all information is lost. With a proper treatment, this would be (0.0, 0.0) instead. We already have this problem. julia 1.0/Complex(0.0, 0.0) NaN + NaN*im (essentially #9531) I agree that (Inf, NaN) should be treated as a representation of complex infinity. (Ref: #5234)
Re: [julia-users] Complex infinity
Ah, I see your concern about the singleton ComplexInf type. This sounds just like the very same design issue behind NA in DataArrays and Nullables in 0.4.
Re: [julia-users] Complex infinity
Another possibility that occurs to me is to redefine the == equality comparison. Currently this is defined as ==(z::Complex, w::Complex) = (real(z) == real(w)) (imag(z) == imag(w)) but for some purposes it may be sufficient to redefine this as ==(z::Complex, w::Complex) = if isinf(z) isinf(w) return true else return (real(z) == real(w)) (imag(z) == imag(w)) end to model the topology of the Riemann sphere. With this redefinition: julia Complex(Inf, 0.0) == Complex(Inf, NaN) true julia Complex(Inf, 0.0) == Complex(Inf, 2.0) true
Re: [julia-users] Complex infinity
On Jan 10, 2015, at 18:15 , Jiahao Chen jia...@mit.edu wrote: Ah, I see your concern about the singleton ComplexInf type. This sounds just like the very same design issue behind NA in DataArrays and Nullables in 0.4. Yes, it is. But no matter how efficient this will be -- this approach will always require at least one additional bit of storage to describe the type, whereas a representation purely in terms of two reals won't. This will always be slower. -erik -- Erik Schnetter schnet...@gmail.com http://www.perimeterinstitute.ca/personal/eschnetter/ My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://pgp.mit.edu/. signature.asc Description: Message signed with OpenPGP using GPGMail
Re: [julia-users] Complex infinity
On Jan 10, 2015, at 18:22 , Jiahao Chen jia...@mit.edu wrote: On Sat, Jan 10, 2015 at 6:12 PM, Erik Schnetter schnet...@gmail.com wrote: This may lead to problems with nans. For example, 1.0/(0.0, 0.0) is (inf, nan), and 1.0/(inf,nan) is (0.0, nan); multiply by (1, im) yields (nan, nan), and all information is lost. With a proper treatment, this would be (0.0, 0.0) instead. We already have this problem. julia 1.0/Complex(0.0, 0.0) NaN + NaN*im (essentially #9531) I agree that (Inf, NaN) should be treated as a representation of complex infinity. (Ref: #5234) This will either require normalizing the result of complex operations, or normalizing the input. -erik -- Erik Schnetter schnet...@gmail.com http://www.perimeterinstitute.ca/personal/eschnetter/ My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://pgp.mit.edu/. signature.asc Description: Message signed with OpenPGP using GPGMail
