I don't think there's a problem with high-dimensional data as long as you are using a 64-bit octave_idx_type which just about everyone is these days.  If you want to use 32-bit types on a Raspberry Pi you can do so, but that's not really a platform for high performance numerical computing.

The reason is that the number of elements on the squashed dimension is

sz = size (data);
nel = prod (sz(squashed_dimensions));

In order to overflow the octave_idx_type which is int64_t the number of elements would have to exceed 2^63 - 1.  If this was a double array then each element would take 8 bytes and this would be 2^66 bytes.  I don't know any machine that has that type of memory.  Hence, there would have been an out-of-memory error when the attempt was made to create the array well before any squashing could happen.

On question 1, I found it easier to understand the reshape process if the dimensions were at the end.  When indexing, you can use the magic colon ':' to collapse all remaining dimensions and this is what the code reminded me of.  For example,

x = reshape (1:24, 2,3,4)
x =

ans(:,:,1) =

   1   3   5
   2   4   6

ans(:,:,2) =

    7    9   11
    8   10   12

ans(:,:,3) =

   13   15   17
   14   16   18

ans(:,:,4) =

   19   21   23
   20   22   24

# Now preserve index 1 (rows) and collapse dimensions 2 and 3 into final dimension (2 = columns)
x([1 2], :)
ans =

    1    3    5    7    9   11   13   15   17   19   21   23
    2    4    6    8   10   12   14   16   18   20   22   24

# Equivalently, use reshape
reshape (x, 2, 3*4)
ans =

    1    3    5    7    9   11   13   15   17   19   21   23
    2    4    6    8   10   12   14   16   18   20   22   24

When the dimensions are not at the end the reshape vector looks a bit awkward.  From median, the code is

      ## Reshape all vecdims into dim1
      num_dim = prod (szx(dim));
      szx(dim) = [];
      szx = [num_dim, ones(1, numel(dim)-1), szx];
      x = reshape (x, szx);

I find the need to insert a vector of ones awkward.  But, it isn't a super big deal and that's why I left it alone in median because it looked like it was going to involve a fair number of changes to calls to sum() etc. where the dimension would no longer be 1 but a variable 'dim'.  I'm okay, though, if we want to rationalize this.

On question 2, that's weird that there is a performance penalty for using true/false.  Was the performance penalty in the assignment 'flag = true' or was it in the evaluation of 'if (flag)'?  I can kind of see that 'flag = true' might be slower because true() is a function, whereas a number like 0 or 1 is directly understood by the parser.  Still, I'm surprised it makes that much of a difference.

I checked in a changeset that updates the code to conform to best practices for Octave core.  See http://hg.savannah.gnu.org/hgweb/octave/rev/6a2638cbea96.

There are a number of separate bug reports that address pieces of bringing the functions named above up to full matlab compatibility.  primary concerns are:
- proper empty input handling
- compatible vector dimension handling
- "all" dimension option handling
- nanflag handling (option to ignore nan values)
- compatible outtype handling

The statistics package has shadowed these functions for octave core and very recently all four functions were made fully matlab compatible.  mean had been ported from statistics to core previously, but compatibility was still incomplete and there was forked development.

I've worked with the statistics package to bring the functions into a compatible state, including avoiding any mean regressions, and the plan is now to merge those changes back into core octave.  Moving forward then statistics can just incorporate any changes to core and shadow them for older versions of octave (targeting v9.1)

In addition, these versions of mean and median should solve bug #54567 regarding large int overflow and double precision limits.

The functions need some header and other text rework to match core function style.  Once I do that cleanup I'll attach those versions to this report and push a patch for review after i verify they don't seriously break any calling functions in core. I'll also try to go through the pile of related bug reports and update/close/dependency-mark them accordingly.

Note: i've done some optimization and speed testing, but the extra input handling for the functions does add overhead. some advice on refactoring to buy back some of that performance would be appreciated.