bugGNU Octave - Bugs: bug #31974, Interpreter fails to create...

 
 

bug #31974: Interpreter fails to create complex numbers with Inf complex parts

Submitted by:  Rik <rik5>
Submitted on:  Sat 25 Dec 2010 03:53:34 PM UTC  
 
Category: InterpreterSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: Need InfoAssigned to: None
Originator Name: Open/Closed: Open
Release: devOperating System: GNU/Linux

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Thu 06 Jan 2011 07:37:24 PM UTC, comment #27:

BTW, the other bug report is https://savannah.gnu.org/bugs/?32053

As with that report, I think we should move any further discussion about this issue to the maintainers list since that will maybe get more people involved in the discussion.

John W. Eaton <jwe>
Project Administrator
Thu 06 Jan 2011 07:28:38 PM UTC, comment #26:

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?

John W. Eaton <jwe>
Project Administrator
Thu 06 Jan 2011 07:00:25 PM UTC, comment #25:

If we're betting, I think lazily creating and destroying the real and imaginary parts would be a very easy optimization for them.

Perhaps you're right though that Matlab's coders wouldn't have stumbled onto this, but I guess they have had 40+ years... maybe they didn't even intentionally do it and their JIT is just that good(TM) and comes up with it on its own ;)

Judd Storrs <judd>
Thu 06 Jan 2011 06:46:47 AM UTC, comment #24:

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

John W. Eaton <jwe>
Project Administrator
Wed 05 Jan 2011 08:05:17 PM UTC, comment #23:

After I tested this in an old Matlab version (R2007b), I realized that R2007b doesn't behave the way R2010a is reported to below:

But below it is claimed

So, Matlab seems to have changed since the version I have access to.

FWIW this is what I see in the old Matlab version:

So, after mxGetPr is called, the real part exists. It may of course be that Matlab creates that field on demand either before mxGetPr() or before handing off any value to MEX.

As for the rest, this is what I see in R2007b, but I doubt it is helpful because of the difference I noted above.

Judd Storrs <judd>
Wed 05 Jan 2011 07:06:29 PM UTC, comment #22:

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..."

John W. Eaton <jwe>
Project Administrator
Wed 05 Jan 2011 07:04:12 PM UTC, comment #21:

Judd, I think that Matlab always allocates storage for the real part. If you have access to Matlab, try this mex file:

and call it with

etc. Can you come up with a way to make Matlab not allocate the real part of a numeric value? That would be interesting, because it would mean that implemented pure imaginary values. But I don't think they have done that. Instead, I think there are some limited number of places where values with zero real parts are treated as if they are pure imaginary. Unfortunately, I think this special treatment also applies to computed values, not just ones explicitly defined with the special values "i" or "j". For example, what happens with

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

This topic of having a pure imaginary type has come up before. If it is really important to have this behavior for pure imaginary values, then I think we need to have it consistently in Fortran, C++, C, etc., not just in some portions of programs like Octave or Matlab. Do we really want that kind of inconsistent behavior?

John W. Eaton <jwe>
Project Administrator
Mon 03 Jan 2011 04:14:52 PM UTC, comment #20:

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.

Judd Storrs <judd>
Mon 03 Jan 2011 08:13:32 AM UTC, comment #19:

Matlab gives me:

>> f (Inf, i)


ans =

0 + Infi

ans =

0 - Infi

Søren Hauberg <hauberg>
Mon 03 Jan 2011 07:02:17 AM UTC, comment #18:

What does Matlab do for

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.

John W. Eaton <jwe>
Project Administrator
Sun 02 Jan 2011 11:57:07 AM UTC, comment #17:

It is also possible that Matlab changed the behavior of realcomplex just to make this thing work. I would advise to avoid doing such a general (and kinda stupid) change and only implement the constant folding optimization as suggested by JWE, i.e. a+bi -> complex(a, b) where a, b are scalar constants, because that's what most people mean. Don't change anything else.
It does mean a sort of inconsistency, but not quite, because we can say that the expression is simply parsed as a constant.

Jaroslav Hajek <highegg>
Thu 30 Dec 2010 03:45:43 PM UTC, comment #16:

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...

John W. Eaton <jwe>
Project Administrator
Thu 30 Dec 2010 03:32:39 PM UTC, comment #15:

I was a bit unsure about the test from JWE (one of vectors needs a transposition), but I get:

>> [Inf; 0] * [i, 0]


ans =

NaN + Infi NaN + NaNi
0 0

>> [Inf, 0] * [i; 0]


ans =

NaN + Infi

>> [Inf, 0] .* [i, 0]


ans =

0 + Infi 0

Søren Hauberg <hauberg>
Thu 30 Dec 2010 03:20:46 PM UTC, comment #14:

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

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?

John W. Eaton <jwe>
Project Administrator
Thu 30 Dec 2010 02:25:04 PM UTC, comment #13:

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


ans =

0 + Infi 0 + Infi
0 + Infi 0 + Infi

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


ans =

0 + Infi 0 + Infi
0 + Infi 0 + Infi

Søren Hauberg <hauberg>
Thu 30 Dec 2010 01:56:14 PM UTC, comment #12:

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.

David Bateman <dbateman>
Project Member
Thu 30 Dec 2010 11:48:11 AM UTC, comment #11:

I get

>> Inf(2,2) * j


ans =

0 + Infi 0 + Infi
0 + Infi 0 + Infi

>> Inf(2,2) .* j


ans =

0 + Infi 0 + Infi
0 + Infi 0 + Infi

>> Inf * j(2,2)

??? Error using ==> j
Too many input arguments.

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

??? Error using ==> j
Too many input arguments.

>> Inf * sparse(j(2,2))

??? Error using ==> j
Too many input arguments.

>> Inf * sparse(1)


ans =

(1,1) Inf

Søren Hauberg <hauberg>
Thu 30 Dec 2010 11:37:42 AM UTC, comment #10:

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.

David Bateman <dbateman>
Project Member
Thu 30 Dec 2010 12:35:12 AM UTC, comment #9:

In Matlab I see:

>> x = Inf;
>> x * j


ans =

0 + Infi

>> j = complex (0, 1);
>> Inf * j


ans =

0 + Infi

>> fcn * j


ans =

0 + Infi

Søren Hauberg <hauberg>
Wed 29 Dec 2010 11:11:29 PM UTC, comment #8:

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.

David Bateman <dbateman>
Project Member
Wed 29 Dec 2010 08:42:43 PM UTC, comment #7:

I get

>> Inf * 4j


ans =

0 + Infi

Søren Hauberg <hauberg>
Wed 29 Dec 2010 08:18:44 PM UTC, comment #6:

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.

David Bateman <dbateman>
Project Member
Wed 29 Dec 2010 07:49:44 PM UTC, comment #5:

In Matlab I get:

>> NaN*j


ans =

0 + NaNi

Søren Hauberg <hauberg>
Wed 29 Dec 2010 07:21:36 PM UTC, comment #4:

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.

Rik <rik5>
Project Administrator
Tue 28 Dec 2010 11:30:57 PM UTC, comment #3:

Using Matlab 7.10.0.499 (R2010a) I see

>> 1+Inf*j


ans =

1.0000 + Infi

>> 1+Infj

??? Undefined function or variable 'Infj'.

>> 1 + Inf*j


ans =

1.0000 + Infi

>> 1 + Infj

??? Undefined function or variable 'Infj'.

>> str2num('1+Inf*j')


ans =

1.0000 + Infi

>> 1 + Inf*i


ans =

1.0000 + Infi

>> str2double ('1 + Inf*i')


ans =

1.0000 + Infi

>> complex(1, Inf)


ans =

1.0000 + Infi

Søren Hauberg <hauberg>
Tue 28 Dec 2010 10:18:43 PM UTC, comment #2:

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.

Rik <rik5>
Project Administrator
Tue 28 Dec 2010 08:29:42 PM UTC, comment #1:

Just "Inf i" is sufficient to show this behavior. However, I'm not sure this is a bug. Read carefully what you've written and how Octave will treat it. That is "1 + Infi". Octave will parse this as "(1) + (Inf * (i))". The "(i)" equals "0 + i", which when multiplied by "Inf" has "Inf * 0" in the real part and guess what "Inf * 0" is "NaN"

What does Matlab do in this case? I suspect it does the same as Octave. If it does I propose to mark this bug as invalid.

D.

David Bateman <dbateman>
Project Member
Sat 25 Dec 2010 03:53:34 PM UTC, original submission:

The code below shows the problem.

Rik <rik5>
Project Administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by judd (Posted a comment)
  • -unavailable- added by highegg (Posted a comment)
  • -unavailable- added by jwe (Posted a comment)
  • -unavailable- added by hauberg (Posted a comment)
  • -unavailable- added by dbateman (Posted a comment)
  • -unavailable- added by rik5 (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follows 1 latest change.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 28 Dec 2010 08:29:42 PM UTCdbatemanStatusNone=>Need Info

    Back to the top


    Powered by Savane 3.1-cleanup