bugGNU Octave - Bugs: bug #61968, pow2 with two arguments is not...

 
 

bug #61968: pow2 with two arguments is not Matlab compatible

Submitter:  Kai Torben Ohlhus <siko1056>
Submitted:  Tue 01 Feb 2022 05:41:45 AM UTC
   
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  siko1056
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 09 Feb 2022 12:42:44 AM UTC, comment #26: 

That file was deleted in the repository. Updating an installation with `make install` basically copies the new files to the installation prefix. It doesn't delete files that might be already installed in that prefix.
Consider wiping that prefix before installing a newer version of Octave. Or delete that file manually.

Markus Mützel <mmuetzel>
Group administrator
Tue 08 Feb 2022 11:59:57 PM UTC, comment #25: 

Is this warning when starting Octave 8 related to this bug report? If it is intended, can it be suppressed?


warning: function .../octave/8.0.0/m/specfun/pow2.m shadows a built-in function


Anonymous
Mon 07 Feb 2022 07:20:49 AM UTC, comment #24: 

I'm not so certain. That document also says §26.8(4):

> The contents of these headers are the same as the Standard C library headers <math.h> and <stdlib.h> respectively, with the following changes:


I'm not sure what of the following that means:
- The functions as standardized for C99 are imported (with C linkage).
- Functions with the same name and signature as standardized in C99 are defined in the `std` namespace (with C++ linkage).
- Both of the above.

Given that the (less ambiguous) standard documentation for C++14 and newer explicitly lists those functions in the `std` namespace, I'd guess that it is probably the last. That'd mean that those functions should exist in the `std` namespace according to that document.

All of this discussion is quite academic though because (like we both already wrote) those functions aren't implemented in the `std` namespace for modern compilers. That might be an oversight on the part of the compiler vendors. Or it might be a conscious decision on their part to not implement those functions (deviating from the standard). I'd understand that decision because these functions don't "feel like C++". But I didn't find any information on that...

Markus Mützel <mmuetzel>
Group administrator
Mon 07 Feb 2022 12:52:06 AM UTC, comment #23: 

Thank you for the additional csets Markus and Rik 🙂

Regarding comment #21: xelem should be okay here too.  In the first branch "e" is a dense matrix and in the second a sparse matrix.  Maybe I was overcautious when facing Sparse(rrors).

Regarding comment #22: Your understanding of the standard is correct and as written in comment #19 (2.) There is no `std::ldexpf` and `std::ldexpl` to find.

http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3690.pdf

§ 26.8 C library (8) [p. 993ff.]

In addition to the double versions of the math functions in <cmath>, C++ adds float and long double overloaded versions of these functions, with the same semantics.


§ 26.8 C library (8)

The added signatures are:
...
float ldexp(float, int);
...
long double ldexp(long double, int);
...


Kai Torben Ohlhus <siko1056>
Group Member
Sun 06 Feb 2022 04:52:41 PM UTC, comment #22: 

Regarding `std::ldexpf` and `std::ldexpl`: I can't find those in the explicit list of functions that should be defined when including `<cmath>`. See 26.8 of the C++11 standard (pages 993 ff.):
http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3690.pdf
(They might be implicitly included by one of the paragraphs though?)

Neither for C++14 (pages 1004 ff.):
http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2014/n4296.pdf

It does appear in C++17 though (p. 1018):
http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2017/n4713.pdf

And in the draft for C++20 (p. 1211):
https://isocpp.org/files/papers/N4860.pdf

Might be an error on cppreference or I don't understand the standard documentation...

However, none of this matters because we don't need those functions anyway.

Markus Mützel <mmuetzel>
Group administrator
Sun 06 Feb 2022 03:17:55 PM UTC, comment #21: 

Thanks for clarifying the difference between xelem and elem. I might have confused that with the difference between operator() and elem.
I still don't understand why one of the branches would need to be handled differently from all other branches in this respect.
I changed it to use `xelem` in that branch together with a few other changes here:
https://hg.savannah.gnu.org/hgweb/octave/rev/555676697eef

Markus Mützel <mmuetzel>
Group administrator
Sun 06 Feb 2022 02:38:18 PM UTC, comment #20: 
Kai Torben Ohlhus <siko1056>
Group Member
Sun 06 Feb 2022 02:34:15 PM UTC, comment #19: 

Okay, I fixed the remaining issues and I am sorry for the long ldexp(f) discussion.  At some point I started to confuse myself.  In short std::ldexp works:

1. C11 has ldexp and ldexpf and a macro to overload.
2. C++11 only has std::ldexp with proper overload for float and "long double".  No need to call something else and https://en.cppreference.com/w/cpp/numeric/math/ldexp should be corrected to not contain those overloads with "f" and "l" suffix 😓

> I was thinking the only difference between `elem` and `xelem` is that the latter doesn't check for out-of-range indices. By that, it's slightly more efficient than `elem`. If you are sure that you'll not exceed the range of the array, you could use `xelem`. Is there really a difference with respect to "unique data"?


The difference between elem and xelem is best described in the code https://octave.org/doxygen/6.2.0/df/dd0/Array_8h_source.html#l00499


T& elem (octave_idx_type n)
{
  make_unique ();
  return xelem (n);
}


The feature with range checking is "checkelem()".

Kai Torben Ohlhus <siko1056>
Group Member
Sun 06 Feb 2022 09:12:07 AM UTC, comment #18: 

Re comment# 16: IIUC, `::ldexpf` corresponds to the function defined in C99. `std::ldexpf` would correspond to the function defined in C++11. But it looks like the latter isn't implemented on the platforms I checked.

Markus Mützel <mmuetzel>
Group administrator
Sun 06 Feb 2022 09:03:19 AM UTC, comment #17: 

`std::ldexpf` not existing in in GCC nor clang might be a compiler bug, a documentation issue, or an (intentional) deviation from the standard. I didn't find any further information on that.
But honestly, I don't know why anyone would even need to use `std::ldexpf`. Instead, it would be clearer and more in the spirit of C++ to just use `std::ldexp` (and potentially cast the input argument to `float` if that is really what one would like to do).

@Kai: Why do you prefer the C functions `ldexp` and `ldexpf` over the (easier to use) C++ function `std::ldexp`? After all, Octave is mainly written in C++. And your proposed implementation of `Fexp2` is also written in C++ (with some functions mixed in from C for some reason I'd like to understand).

> pow2(f,e) with two arguments calls ldexp.  This is documented Matlab behavior (comment #0): discussions forbidden 😅


As I'm understanding their documentation, they are describing the function by comparing it to `ldexp`. They don't actually write that the function is implemented with `ldexp`...
But that's just a minor detail.
It might still be interesting to check how much of a performance difference there is between using `std::exp2` or `std::ldexp` for different compiler or STL implementations.

> Somehow in Octave the C symbols leak (intentionally?) into the C++ "octave::" namespace.  I got hard errors using std::exp2, etc.


Tbh, I don't understand what you are trying to say...


I was thinking the only difference between `elem` and `xelem` is that the latter doesn't check for out-of-range indices. By that, it's slightly more efficient than `elem`. If you are sure that you'll not exceed the range of the array, you could use `xelem`. Is there really a difference with respect to "unique data"?

Markus Mützel <mmuetzel>
Group administrator
Sat 05 Feb 2022 04:23:08 PM UTC, comment #16: 

Hmm, ::ldexpf works in C++ with GCC.  Is it a bug that ldexpf is not imported into the std:: namespace or is cppreference wrong in stating that it should be?

John W. Eaton <jwe>
Group administrator
Sat 05 Feb 2022 04:00:08 PM UTC, comment #15: 

In C++ with function overloading, if you have included <cmath> then you should get the correct float or double function if you use the name ldexp depending on the type of the argument.  In C99, I guess function overloading is done with a magic macro.  In C++ ldexp should not be a macro but should refer to actual overloaded function names.

John W. Eaton <jwe>
Group administrator
Sat 05 Feb 2022 10:28:24 AM UTC, comment #14: 

Thanks again for the review @mmuetzel.

comment #12 and comment #13: Copyright, BIST pattern, and warning ID will fixed.  Thanks for the hints.

comment #12 and comment #9:

> On many implementations, ldexp is less efficient than multiplication or division by a power of two using arithmetic operators.


There must be a distinction.  pow2(x) with single argument calls exp2(x) or pow(2,x) and is equal to "2.^x".  For (non-complex) double and float (single) input we are even better now and pow2(x) calls the more specialized exp2(x) function explicitly.

pow2(f,e) with two arguments calls ldexp.  This is documented Matlab behavior (comment #0): discussions forbidden 😅

comment #9: More off-topic: Why are we using C functions (like ::exp) in many cases instead of their corresponding C++ functions (like std::exp)?

Somehow in Octave the C symbols leak (intentionally?) into the C++ "octave::" namespace.  I got hard errors using std::exp2, etc.

comment #12: "xelem" vs. "elem", etc., I reviewed the usages.  Sometimes it is safer to call "elem" to ensure a unique data instance.  Input data structures created outside the pow2 function are beyond the functions control (thinking of race conditions in the future).

comment #12: Why not use `std::ldexp` in both cases (single and double precision).

That is true, for C++ ldexpf does not exist and double/float overloads named "ldexp" exist.

C no error: https://godbolt.org/z/GoTYGrWE5
C++ error: https://godbolt.org/z/n4r7xxrfc

The wiki reference is a little misleading:

https://en.cppreference.com/w/cpp/numeric/math/ldexp
https://en.cppreference.com/w/c/numeric/math/ldexp

As I am calling C functions, I try avoid implicit float->double typecasts and call the respective function:

https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/flt-32/s_scalbnf.c;h=d521d9a3035d572c43454ac17ea5477c8f2b18d7;hb=refs/heads/master

https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_scalbn.c;h=4491227f3e3b5cf431564146b4aadc9cc206339e;hb=refs/heads/master

Do you know it is safe to ignore this float/double distinction in C?

Kai Torben Ohlhus <siko1056>
Group Member
Fri 04 Feb 2022 10:20:00 AM UTC, comment #13: 

More nit: Afaict, we usually don't claim copyright for the past if a new file is added.

Markus Mützel <mmuetzel>
Group administrator
Fri 04 Feb 2022 10:13:22 AM UTC, comment #12: 

I haven't tested anything. Just a few remarks from reading the patch:

> On many implementations, ldexp is less efficient than multiplication or division by a power of two using arithmetic operators.


Did you check if that has a measurable difference?

What is the difference between `pow3` and `pow2` in your test?


+%! f = str2num (char ({'1/2','pi/4','-3/4','1/2','1-eps/2','1/2'}'));


That coding pattern looks odd. Why not use something like this directly?


+%! f = [1/2, pi/4, -3/4, 1/2, 1-eps()/2, 1/2];



+      map_2_xldexp (y, f, e, ldexpf);


Why do you use that pattern of passing a function handle? Why not use `std::ldexp` in both cases (single and double precision). The standard defines overloads for both anyway:
https://en.cppreference.com/w/cpp/numeric/math/ldexp


+    warning_with_id ("pow2:imag_ignored", "pow2: Imaginary part is ignored.");


Warning ids in core Octave should start with "Octave:" and use "-" instead of "_".


+          int ee = static_cast<int> (e.elem (f.ridx (i), col - 1));


You are using `xelem` in most of the other cases. Here, too?

Markus Mützel <mmuetzel>
Group administrator
Fri 04 Feb 2022 09:26:21 AM UTC, comment #11: 

Thanks for the input.  I overhauled the patch (file #52801) to only include a single file in libinterp/corefcn.  I think this pattern is much better to keep the implementation at a single place, than spread about all libinterp.

This way I can also use this code in my own project as single file 🙂

Regarding comment #8.  I also find the current implementation of pow2() too simple to be true.  If "2.^x" was really required, I find this statement more expressive than pow2(x).  One does not need an extra function for this.

Unfortunately, pow2() does not do what Matlab intended it to do.  Further it is inefficient for sparse input, as "2.^e" fills-up merciless to a dense NxN matrix, even though the final "f .* 2 .^ e" is very sparse 😓


>> e = rand(500); tic; y = pow3(e); toc, tic; y = pow2(e); tocElapsed time is 0.021579 seconds.
Elapsed time is 0.00504899 seconds.
>> N = 1e3; f = sprand (N, N, 1/N); e = sparse(N,N); tic; y = pow3(f,e); toc, tic; y = pow2(f,e); toc
Elapsed time is 0.0445991 seconds.
Elapsed time is 0.000132799 seconds.


A final (a little exotic) issue I care about is of numerical stability nature of exp2(x) vs. pow(2,x): https://godbolt.org/z/zs3s9dqfs

Again, the great advantage of the patch is it finally makes pow2() do what Matlab is doing for the two argument call.

If no further objections, I would push the patch next week.

Kai Torben Ohlhus <siko1056>
Group Member
Tue 01 Feb 2022 05:01:42 PM UTC, comment #10: 

Oops. Meant to write `isfloat ()` (not `is_float ()`).

Markus Mützel <mmuetzel>
Group administrator
Tue 01 Feb 2022 04:56:38 PM UTC, comment #9: 

Some more comments on the patch:

+  if (! args(0).is_single_type () && ! args(0).is_double_type ())


Why not `! args(0).is_float ()`?


+  if (! args(1).is_single_type () && ! args(1).is_double_type ())


Same.


+      FloatNDArray e = args(0).float_array_value ();


Typo? Should be `args(1)` instead?

IIUC, `ov.exp2 ()` and `Fexp2 (ov)` would do two different operations. (The former essentially calls the C function from math.h `exp2`, the latter `ldexp`.) Is that what we want? Or do we just need to live with that difference...

Also from https://en.cppreference.com/w/c/numeric/math/ldexp

> On many implementations, ldexp is less efficient than multiplication or division by a power of two using arithmetic operators.


More off-topic: Why are we using C functions (like ::exp) in many cases instead of their corresponding C++ functions (like std::exp)?

Markus Mützel <mmuetzel>
Group administrator
Tue 01 Feb 2022 04:21:16 PM UTC, comment #8: 

"As nice side effect the function will move from m-code to C++ code."

It's fine to make this change, but I'm not sure that defining more functions in C++ is a good thing.  I understand there are performance issues, but I'd rather have fewer things defined in C++ unless they are simple wrappers around existing library code.  Another option would be to provide wrappers for ldexp and scalbn and using those functions to define pow2.

John W. Eaton <jwe>
Group administrator
Tue 01 Feb 2022 04:16:13 PM UTC, comment #7: 

I've never like that whitespace is significant inside the square brackets used for defining matrices.  But we have to have it for compatibility.

John W. Eaton <jwe>
Group administrator
Tue 01 Feb 2022 12:27:14 PM UTC, comment #6: 

To answer my own question from comment #4:

>> e = [1 + 1i, 1]
e =
   1 + 1i   1 + 0i
>> e = [1+ 1i, 1]
e =
   1 + 1i   1 + 0i
>> e = [1 +1i, 1]
e =
   1 + 0i   0 + 1i   1 + 0i
>> e = [1+1i, 1]
e =
   1 + 1i   1 + 0i


I always found the parsing rules around spaces in arrays/vectors pretty confusing. (They do make sense when you take the time to think about it. But I'm lazy.😉)
That's why I usually try to avoid them (unless next to a comma).

Markus Mützel <mmuetzel>
Group administrator
Tue 01 Feb 2022 12:26:30 PM UTC, comment #5: 

Thanks for the idea, I will find a proper solution.

Regarding the question of comment #4, it should be a complex and a real number, thus 2 elements in a complex row vector.


e =

   1 + 1i   1 + 0i


Kai Torben Ohlhus <siko1056>
Group Member
Tue 01 Feb 2022 10:05:30 AM UTC, comment #4: 

I understand your personal motivation for having this as soon as possible. But I'm not sure if it is a good idea to have this change in behavior for a minor release either.
Maybe you could use an .oct file implementation for your project in the meantime?

More on topic: While scrolling through your patch, I noticed this test:

+%!test
+%! f = 1;
+%! e = [1 + 1i, 1];
+%! z = [2, 2];
+%! assert (pow2 (f, e), z);


Haven't tested yet. But how many elements does `e` have after that assignment? 2 or 3?

Markus Mützel <mmuetzel>
Group administrator
Tue 01 Feb 2022 09:56:02 AM UTC, comment #3: 

The operation turned out to be more complex than I initially thought.  An almost ready to apply patch attached (file #52788).

Maybe dev is better for now.  The dilemma before an release: now or in another year? 😅  One project of mine will use it... maybe for Octave 7.2.0? 😇

Kai Torben Ohlhus <siko1056>
Group Member
Tue 01 Feb 2022 09:47:22 AM UTC, comment #2: 

IIUC, this isn't a regression. And Octave's current behavior looks reasonable to me as well (albeit deviating from Matlab).
Is there a reason this is targeting Octave 7?

Markus Mützel <mmuetzel>
Group administrator
Tue 01 Feb 2022 05:47:31 AM UTC, comment #1: 

One addendum for TWO argument call and non-integer exponents in Matlab R2021a:


>> f = 2; e = [2, 2.1, 2.2, 2.4, 2.5, 2.8]; y1 = pow2(f,e), y2 = f.*2.^e
y1 = 8     8     8     8     8     8
y2 = 8.0000    8.5742    9.1896   10.5561   11.3137   13.9288


Kai Torben Ohlhus <siko1056>
Group Member
Tue 01 Feb 2022 05:41:45 AM UTC, original submission:  

From the Matlab documentation https://www.mathworks.com/help/matlab/ref/pow2.html:

> The syntax Y = pow2(X,E) corresponds to the ANSI® C function ldexp() and the IEEE® floating-point standard function scalbn(). The result Y is computed quickly by simply adding E to the floating-point exponent of X.


Octave so far simply computes https://hg.savannah.gnu.org/hgweb/octave/file/66456820ff59/scripts/specfun/pow2.m#l57


Y = X .* 2 .^ E


Some test cases for verification (fails in Octave 7):


Yexpect = str2num (char ({'1','pi','-3','eps','realmax','realmin'}'))
X = str2num (char ({'1/2','pi/4','-3/4','1/2','1-eps/2','1/2'}'));
E = [1 2 2 -51 1024 -1021]';
Y = pow2 (X, E)
isequal (Y, Yexpect)  % true == 1


Especially the complex part is ignored in the TWO argument call.  Sparse and Single type are preserved in the output.


>> f = [1+1i, 1]; e = 2; y1 = pow2(f,e), y2 = f.*2.^e
Warning: Imaginary part is ignored.
y1 = 4     4
y2 = 4.0000 + 4.0000i   4.0000 + 0.0000i

>> f = 1; e = [1 + 1i, 1 + 1i]; y1 = pow2(f,e), y2 = f.*2.^e
Warning: Imaginary part is ignored.
y1 = 2     2
y2 = 1.5385 + 1.2779i   1.5385 + 1.2779i


I am working on a patch.  As nice side effect the function will move from m-code to C++ code.

Kai Torben Ohlhus <siko1056>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #52801:  bug61968_v2.patch added by siko1056 (13KiB - text/x-patch)
file #52788:  bug61968.patch added by siko1056 (14KiB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Improved the code.)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by siko1056 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-02-07 siko1056 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-02-07 siko1056 Carbon-Copy- Added rik5
    2022-02-06 siko1056 StatusPatch Submitted Ready For Test
    2022-02-04 siko1056 StatusIn Progress Patch Submitted
    2022-02-04 siko1056 Attached File- Added bug61968_v2.patch, #52801
    2022-02-01 siko1056 Release7.0.90 dev
    2022-02-01 siko1056 Attached File- Added bug61968.patch, #52788

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code