Strange inconsistent in complex algebra

Mathmatica is not choosing a different branch than Julia. Since Julia uses floating point arithmatic, there is no branch choice here. 0.0 uses one branch, -0.0 uses a different one. This is unambiguously THE correct choice.

1 Like

It is also possible in practice — just define your own number type with the relevant methods. There are packages with custom number types, eg

2 Likes

Thank you @Tamas_Papp very much. The SaterIntegers.jl is what I need

You can try Mathematica first. I am confusing on the Julia output at the beginning is because that the Julia output does not the same as the Mathematica output.

The following is the Mathematica code. You can see the difference with Julia. From the Math or science point of view. This is the natural, convenient and effecient Complex function. Zero should be an exact zero. And Zero times any number should equals Zero obviously. This is the principle of nature.

Yeah, there’s a reason mathmatica is slow. It’s because processors don’t like this and mathmatica has to special case lots of stuff to make it look like they do. No language with these behaviors will be reasonably fast on modern processors.

1 Like

Edit: Somehow I didn’t see all the previous answers before I posted this.

Both answers are valid. A simpler version of your example is

julia> Complex(-1.0, 0.0)^(0.5)
0.0 + 1.0im

julia> Complex(-1.0, -0.0)^(0.5)
0.0 - 1.0im

Any algorithm that takes roots must choose which “branch” to take; this is arbitrary. What is curious in this example is that different roots are picked even though mathematically -0im equals 0im. However, the IEEE floating point standard (which Julia follows) encodes -0.0 differently than 0.0. Thus

julia> 0.0 == -0.0
true

julia> 0.0 === -0.0
false

Thus technically Complex(-1.0,0.0) is in the upper half of the complex plane while Complex(-1.0,-0.0) is in the lower half, which leads to the root-taking algorithm picking a different root.

No, 0*z == 0 is still true for any z except Inf or NaN.

1 Like

Appealing to grand philosophical statements like this usually has very low returns in the Julia community, for two reasons:

  1. You can define your own number type with your desired semantics. You never have to convince or coordinate with others about this, it will be a first-class citizen with near-zero overhead (other than the cost you pay for not relying on IEEE semantics, but that’s a given here).

  2. Compelling reasons are needed to make the semantics for the built-in default types deviate from what is common in scientific computing (ie pretty much what you see in C/Fortran).

So if you want a number type like Mathematica’s, just create a package that does it and you are all set.

4 Likes

That is a good point. But Inf or NaN is not a number. The problem in my example is 0*(-1.0) is became -0.0. My solution to this is define a Zero such that Zero*(any number)==Zero, and redefine complex() function to generate the Zero im. It will not affect the efficient and generate a nature output. So Julia can solve this problem.

I do not want to Julia developer to change 0::Int to meet my requirement. My claim is from the costemer or user point of view, not from the developer. I think most Mathematica users will meet similar problem when they want to try Julia.

When I call C or Fortan or Jave or Python program , I know this problem long time ago. But it does not matter for these Languages because I know they are not math freindly launguage. But Julia has a math-friendly syntax which is one of major advantages. I suggest Julia devloper to add a True zero in the core math packgage. It will not effect current users. And it is helpful in attracting many old Mathematica users.

1 Like

I am not sure if you are aware of false being a “strong zero” already:

help?> false
  Bool <: Integer

  Boolean type, containing the values true and false.

  Bool is a kind of number: false is numerically equal to 0 and true is 
  numerically equal to 1.
  Moreover, false acts as a multiplicative "strong zero":

  julia> false == 0
  true
  
  julia> true == 1
  true
  
  julia> 0 * NaN
  NaN
  
  julia> false * NaN
  0.0
1 Like

It is nice to know that. But I try the following code. It is not a strong zero, at least not so strong.

Just to clarify: it actually does hold that -1 * 0.0 == 0.0 (but not ===).

The way I see it, the thing that you take issue with is not with whether 0.0 is a multiplicative zero in the mathematical sense or not. Rather, it’s the following:

There are functions (log, fractional powers) for which f(x) != f(y) even when x == y.

To be fair to julia, this only happens when f is a multi-valued function with a branch cut and when x and y are exactly on the branch cut. To be fair to you, this can definitely take someone by surprise and then it’s an annoying gotcha. @Tamas_Papp shared an interesting paper about the design reasons and @stevengj explained one of the paper’s examples in a comment.

The reason why we get side-tracked about multiplicative zeros is because the branch-cut happens along \Im(z) = 0 but it’s otherwise not very relevant to the issue.

Is that a reasonable summary?

2 Likes

Yes, possibly a signless strong zero would make sense for your use case:

struct SignEater <: Real end
Base.:*(::SignEater, ::T) where {T <: Real} = zero(T)
Base.:*(::SignEater, ::Complex{T}) where {T <: Real} = Complex(zero(T),zero(T))
# NOTE: more methods are needed for ambiguities and I did not do that carefully
const 𝟎 = SignEater()
𝟎 * -Inf # 0.0
𝟎 * 1im  # 0 + 0im

(Incidentally, please post quoted code, not screenshots).

1 Like

Nice code. Thanks. I think it is helpful temporary.

Since -0.0 == +0.0 according to IEEE 754 rules, you still have 0*(-1.0) == 0.

What IEEE arithmetic is doing is keeping track of one extra bit of information alongside the real number 0 — from which side of ℝ did we approach zero? But it is still maintaining the rule that 0*z == 0.

I think most Mathematica users will meet similar problem when they want to try Julia.

Mathematica users will meet the same “problem” when they use machine numbers in Mathematica. Or when they switch to C, Fortran, Python, R, or any other popular language. Floating-point arithmetic rules take some getting used to!

3 Likes

As I said before, Julia is a Math friendly language. But the languages like C, Fortran, Python, R, or nearly all other languages are not. (expect Mathematica and Matlab). So I do not think other languages (especially C and fortran) can avoid this annoying gotcha conveniently. But Julia can, by just introducing a new strong Zero. The code is provided by @Tamas_Papp. Even in my primary julia code style, It just need to add fewer than ten lines code. Another way is define a new Power function, that is more easy to realize (three lines code).

Mathematica keep both machine number and Math number. But for Julia, add the Math number is not a good idea. It really will affect the performence. For general number, some other rules like 1::int +2.0::float= 3::int does lower performence of the Language and such rules are not quite useful in practice. The only exepection number is the math Zero. Introducing the strong math zero of rule like Zero1.0=Zero or Zero(-1.0)=Zero will not affect the performence. Moreover, it is quite helpful in practice. In my expeirence in scientific computation, although 90% problem do not need the strong math zero, but there does exist nearly 10% problem need this strong math zero.

My suggestion as an Julia user is introduce a strong Zero to core math package in Julia. (Not replace 0::int or 0.0::Float. The IEEE arithmetic rule does make sense.).

I can not agree on this point. Such rule is suitable for mature Language like c or Jave but not for Julia. Julia is too new for most people and the devoloper should pay more efforts on the user’s feedback. If adding one more line code can attract one more user, it worth to do it.

Anyway, thanks for the discussion and help me in dealing with this annoying gotcha in Julia style. This problem is solved already.

I’m pretty sure Matlab has the same behavior as Julia also.

Please consider that code just a proof of concept: as the note indicates, there is more work needed to flesh it out (also, I am not endorsing the use case, just wanted to show how easy it was).

There is no reason why it has to live in Base. A package can implement it just fine for those interested.

4 Likes

You are right. +1