See discussion and voting by Octave Maintainers on this thread: https://octave.discourse.group/t/should-octave-support-special-case-code-for-val-i-where-val-inf-nan-for-matlab-compatibility/4341

The conclusion is that Octave will not implement this feature because it introduces an inconsistency between scalar multiplication and vector/matrix multiplication which is deemed confusing.  In addition, the current behavior is only a problem with NaN/Inf construction of complex imaginary values and that has a very low incidence rate which can be worked around by calling complex (xxx, Inf).

Marking as "Won't Fix" and closing report.

update: checked in v7.2.0, v8.0.0 hg id 8b75954a4670, and matlab 2022a.
no change from comment #29

also, for easier copy/paste next time:

1+Inf*i
str2double('1 + Inf*i')
complex(1,Inf)
1+Inf*j
1 + Inf*j
1 + Infj
str2num('1+Inf*j')
123*j
NaN*j
Inf * 4j
x = Inf;x*j
j = complex(0,1);Inf*j

Just updating that the outputs from Comment #29 are still the same in Octave 5.2.0 (checked in Win-10)

just updating that the outputs from Comment #29 are still the same in Octave 4.2.0 and Matlab 2016b

I have distilled the comments on this bug report and made those tests. I will add them to data.cc (under the function complex) but maybe they should be somewhere else?

%!error <undefined> 1+Infj
%!error <undefined> 1+Infi

%!test <31974>
%! assert (Inf + Inf*i, complex (Inf, Inf))
%!
%! assert (1 + Inf*i, complex (1, Inf))
%! assert (1 + Inf*j, complex (1, Inf))
%!
%! ## whitespace should not affect parsing
%! assert (1+Inf*i, complex (1, Inf))
%! assert (1+Inf*j, complex (1, Inf))
%!
%! assert (NaN*j, complex (0, NaN))
%!
%! assert (Inf * 4j, complex (0, Inf))

%!test <31974>
%! x = Inf;
%! assert (x * j, complex (0, Inf))
%! j = complex (0, 1);
%! assert (Inf * j, complex (0, Inf))

%!test <31974>
%! exp = complex (zeros (2, 2), Inf (2, 2));
%! assert (Inf (2, 2) * j, exp)
%! assert (Inf (2, 2) .* j, exp)
%! assert (Inf * (ones (2, 2) * j), exp)
%! assert (Inf (2, 2) .* (ones (2, 2) * j), exp)

%!test <31974>
%! assert ([Inf; 0] * [i, 0], complex ([NaN NaN; 0 0], [Inf NaN; 0 0]))
%! assert ([Inf, 0] * [i; 0], complex (NaN, Inf))
%! assert ([Inf, 0] .* [i, 0], complex ([0 0], [Inf 0]))

%!test <31974>
%! m = @(x, y) x * y;
%! d = @(x, y) x / y;
%! assert (m (Inf, i), complex (0, +Inf))
%! assert (d (Inf, i), complex (0, -Inf))

I decided to sort through "Need Info"s and saw this old bug as still open.  I noticed the original complaint still produces the same result in Octave 4.0.2, different from Matlab 2016a.

The referenced related bug 32053 was closed as fixed. So, is this a Need Info?, should it become a Won't fix? A 'Matlab is wrong'...? I couldn't find a maintainers list thread to follow, but it's 5 years old so...

For posterity, here's are some simpler tests I culled from below with current output:

>> 1+Inf*i
Octave:  NaN + Infi
Matlab:  1 + Infi

>> str2double('1 + Inf*i')
O: 1 + Infi
M: 1 + Infi

>> complex(1,Inf)
O: 1 + Infi
M: 1 + Infi

>> 1+Inf*j
O: NaN + Infi
M: 1+Infi

>> 1 + Inf*j
O: NaN + Infi
M: 1+Infi

>> 1 + Infj
O: error Infj undefined
M: error: Infj undefined

>> str2num('1+Inf*j')
O: NaN + Infi
M: 1 + Infi

>> 123*j
O: 0 + 123i
M: 0 + 123i

>> NaN*j
O: NaN + NaNi
M: NaN + NaNi

>> Inf * 4j
O: NaN + Infi
M: 0 + Infi

>> x = Inf;x*j
O: Nan + Infi
M: 0 + Infi

>> j = complex(0,1);Inf*j
O: NaN + Infi
M: 0 + Infi

and many more...

If you really want to find out, then I think you might be able to cast an mxArray* to char* and look at the bytes, then look at the values returned by mxGetPr and mxGetPI and compare.  That and a little more digging might uncover something about the layout of the mxArray object and when things are allocated.  I don't really care to do that myself as those details don't really matter to me.

What does matter is whether we should treat complex values with zero real part as pure imaginary values so we don't generate NaNs when multiplying a complex value with zero real part by Inf.  Apparently this surprises some people because they don't realize that there are no pure imaginary values, just real and complex.  The MathWorks apparently made this choice for Matlab, for all cases, not just the ones Jaroslav has suggested.  So should we be fully compatible, partly compatible, or not at all with regard to this feature?

I'm getting feedback on another bug report about complex numbers that is arguing against automatic narrowing of complex values that have a zero imaginary part (which is compatible behavior).  So if we narrow there, why not (appear to) narrow for the opposite case of zero real part?  Or why do either, even if it means incompatibility, if not ignoring the zero part is the right thing to do?

Whatever we do, I'm sure someone will not be happy about it.

I don't see the point of just being partly compatible here, and I don't really like the idea that internal computations can have different results from ones that are implemented in the interpreter operators themselves.  But whatever we decide, I don't plan to make any changes in this until after 3.4 is released, so can we stop wasting time on this issue now?

Yes, it's possible that Matlab creates the real part when asked, but I'd bet that is not how things really work.  In any case, the point is that the behavior at the scripting language level is a special case that is not present in the complex type in Fortran or C++, so there will be inconsistencies between what happens for the scripting language and operations in the BLAS (for example).  That was the point I was trying to make with the example of

[Inf, 0] * [i; 0]

Oops, in my previous comment, I wrote

Do you still see the special treatment of the zero imaginary part here, or do you get NaN + Infi?

That should have been "special treatment of the zero REAL part..."

It seems like this may result from the differences in complex number representation in Matlab and in Octave?  Matlab stores the real and imaginary components separately whereas octave interleaves them.

If the components are stored separately, it would be trivial to get Matlab's behavior based on the states of the pointers to the real and imaginary memory chunks.  i.e. NaN*j would create a complex number that only has memory for the imaginary field allocated.  In some cases when Matlab passes the value off to the blas or whatever in some cases it will be forced to first allocate both components of the complex value, losing the internal state and resulting in identical behavior as octave.

So for the "1+Inf*j" Matlab could simply do "(1,unallocated)" + "(unallocated,Inf)" by letting allocated fields simply replace unallocated fields.  This would work without a special parser rule and copy-on-write means there is essentially no penalty for behaving this way.

Approaching this from the interleaved representation (used in octave) it's more difficult because both the real and imaginary parts are allocated simultaneously as a single region of memory that cannot be separated--i.e. octave can't hide state inside allocation patterns.

One approach would be to add a state variable to the complex values that labels the status of the (real,imaginary) components to obtain the semantics of the separate real/imag representation.  Another option would be to add a pure-imaginary array similar to the pure-real arrays we already have.

What does Matlab do for

function f (x, y)
  x * y
  x / y
end

f (Inf, i)

If it produces 0+Infi and 0-Infi, then this is not limited to constant folding, so we have to implement the special case in the operators that might be affected.  I would guess this includes , ., /, ./, \, and .\ for mixed real and complex cases for scalar by scalar and scalar by array operations for both double and single values.  Or, we can use a class derived from std::complex that implements the special cases we want.  Using a derived class like this would probably require fewer changes to Octave since we already have Complex and FloatComplex typedefs, but then it would mean that anywhere the class is used, we get the special behavior, and I'm not sure that's really what we want.

In any case, I think this change should wait until after 3.4 is released.

Sorry, I was thinking inner product, but your tests show that the outer product does the same.  The key thing is that then the operation is done by the BLAS library which doesn't have the special case code.

Also, my intro text was mangled by the rich markup feature...

In mixed complex/real element-by-element operations (either explicit .*, ./, etc., operations, or scalar by matrix operations) we can intercept the operation and pretend that a complex value with 0 real part is a pure imaginary number and skip the 0 * Inf operation.  This should be done in the function that implements the operator.  This may appear to do the "right" thing in some cases, but will be inconsistent with anything that doesn't use our special functions for these mixed ops.  For example, I think the test you want is

[Inf, 0] * [i, 0]

This operation should use the BLAS, so the Inf*i won't be subject to a special definition of a mixed real by complex multiplication operator.

I'm not convinced this is something we should "fix", but if we do, and we are doing it for Matlab compatibility, I think we need to determine exactly how deeply this special case should be handled.  Should it just be done in the operator function(s)?  Should we replace <complex> with our own class that has the special cases?  Do we only have to change the behavior of multiplication or also division?  Other operators?

As j(2,2) doesn't work in matlab one final test. What does

Inf (ones(2,2) j)
Inf(2,2) .* (ones(2,2) * j)

give in matlab? I'm desperately hoping there is a NaN in the real part, as then we'd only have to check if a scalar is purely imaginary.

D.

Frankly, what the hell was mathworks thinking when they implemented this exception to IEEE754..... What a mess..

Does this behavior also hold for matrices or is it particular to scalars? What about sparse matrices? or even scalars stored in sparse matrices? What does

Inf(2,2) * j
Inf(2,2) .* j
Inf * j(2,2)
Inf(2,2) .* j(2,2)
Inf * sparse(j(2,2))
Inf * sparse(1)

give in matllab? If the same behavior holds for matrices, then we are going to take a serious performance hit to special case the imaginary by infinity or NaN case.

D.

I had a look at coding this as it seemed to me that everything could be treated in the function

oct-parse.yy (static tree_expression fold (tree_binary_expression ))

but the problem is harder than I thought.. We need to treat not just "Inf j" but also "Inf 1j" as the parser treats the first as two identifiers and the second as an identifier times a constant. However identifier can be overloaded by user functions and variables, so

Inf = pi;
Inf * j

should return "0 + pi * j" so a simple test of the identifier in the parser can work. A solution for the above I looked at however is stymied by the case

for i = 1 : 10
  d(i) = Inf * j;
  Inf = pi;'t
endfor

as d(1) should be "0 + Inf j" but d(2:10) should be "0 + pi j".

Also another test case for matlab is

x = Inf;
x * j

Does this return "0 + j"  or "NaN + j"? And what about

j = Complex (0, 1);
Inf * j

does it return "0 + j" or "NaN + j"? And finally what about if you create a user function

function y = fcn ()
  y = Inf;
end

and then run

fcn * j

does it return "0 + j" or "NaN + j"? If the above three cases all return "0 + j" then the problem should be treated in ov.cc (do_binary_op) rather then in the parser and the problems above would go away.

D.

One final test for matlab I propose is

Inf * 4j

If that equal "0 + Inf" then the special case in the parser might be simpler as then we probably only need to special case the cases

Inf * constant
NaN * constant
constant * Inf
constant * NaN

where "constant" is already defined in the parser. As I don't know much yacc I hope John steps forward and looks at this.

D.

So Matlab seems to have a special case for Inf*[complex unit].  One additional test is to see whether Matlab behaves correctly for XXX*[complex unit] where XXX is any value or whether 'Inf' is actually special.  My guess is that the parser does something like this.

XXX*j  is translated to function call complex (0, XXX)

rather than translating to XXX * complex (0, 1) as Octave does.

Final test case:
NaN*j

I think this will yield 0 + NaNi in Matlab

I know what Octave is doing, I'm just not sure it's the best interface to present to the user.  There are a number of possible ways to create a number with infinite complex magnitude, but most of them fail with Octave.  Whoever can run Matlab to check its behavior should also try some of these other corner cases.

1+Inf*j
1+Infj
1 + Inf*j // to see whether spaces affect parsing
1 + Infj  // ditto
str2num("1+Inf*j")

The code below shows the problem.

1 + Inf*i
ans = NaN + Infi

str2double ("1 + Inf*i")
ans =    1 + Infi

complex(1, Inf)
ans =    1 + Infi