I implemented many of the ideas in the various comments in to a patch which is attached. The case for int64 is not ideal since it warns when there is potential loss of accuracy, but should just calculate an accurate result. The code is tricky enough that maybe Dan's suggestion of using mean() would be better. First, fix mean() to handle every single case and then let median call mean.
One possibility is to use the algorithm that Dan suggested in comment #11. One could use permute ahead of the algorithm in order to guarantee that the correct dimension was lined up the way the algorithm expects. After calculating the median one would then use ipermute to reverse the change.
(file #44882)

Dan: Right. Just thought that this might be of interest,
and, of course, it is a "signalprocessing" technique.

The Mohanty algorithm is actually for "running" median, or another common terminology, "moving median". It keeps a sorted linked list for which it can easily drop a sample when the window moves past that sample, and add a sorted sample in the list when that sample moves into the window. This would be something good for the signalprocessing toolbox if a moving median doesn't exist already. (Of course, Octave scripting of linked lists isn't very efficient, so written in C++ is preferred.)
Octave's current median() with sorting is pretty efficient. The mean() function definitely needs some work, though. There probably is not a lot of coding, but because of its importance the routine needs attention to detail which could take time.

I added comment #6 just for information. A median algorithm for
very large datasets (as, for example, the LIGO data)
is fully discussed in:
https://dccbackup.ligo.org/public/0027/T030168/000/T03016800.pdf
This reference includes the Ccode for the algorithm.
Also, the current Octave mean() just uses sum()/n which is
not robust nor efficient. Improving this might be a good
project.

BTW, we should at BIST tests for all these cases.
Also, we know mean() probably has the issue we are trying to avoid:
I.e., the large magnitude int64 numbers have a loss of precision in mean() so the result doesn't match the formula of Comment #11.

KT, the following seems to work for int64:
Notice there is one small issue with the fact I had to use the colon indexing, meaning that I'm assuming dim==2. That means you'll have to have an ifelse statement based on "dim" with slightly different formulas for dim=1 versus dim=2. That's unfortunate. I wonder if we should devise an indexing type of function that essentially extracts a submatrix from a matrix the way so many functions can act on "dim", e.g., nth_element(A,k,dim), sum(A,dim), mean(dim).
Also, if this approach proves more efficient than all the min/maxing of the uint64 formula in Comment #10, perhaps use of logical arrays could be used for unit64 as well.

KT, here's a test of the unsigned 64bit integer using min/max functions:
I tested what happens when the unsigned "wrap" occurs. Well, Octave is programmed to not wrap. I was curious what happens with wrapping because sometimes one can devise a formula that allows integer wrapping so long as it is guaranteed to unwrap. Since Octave truncates rather than wraps, there's no tricks that can be done.
OK, so let's think about an int64 formula then...

That's fine with me.
On the matter of two numbers, for what it is worth, it seems the floating point conversion (then sum, then divide) is the fastest. Consider
So the additional conversion to uint32, rounding and divisionbybitshift comes out slower than converting to float and doing the math. Even when x and y vectors are not converted, it is still on the order of converting to float first. I would think that in assembly code that the unsigned int conversion might be "free", i.e., just a case of a different machine instruction that treats the higher register bits as an extension from 8 bits to 32 bits. I suppose what the floating point instructions have going for them is inherent rounding.

Even though it makes for clear code to call mean(), we should not do it. This is the mean of only two numbers and the overhead involved in calling an Octave function is significant. For core Octave we do prioritize speed and efficiency to some degree over readability of code. So the calculation should be be done in median.m, and probably needs an exception case for (u)int64 values. I can review any amended changeset.

Nice article on mean computation... This is why I wondered if it were better to put this case of 64bit average in the mean() script and call mean() inside median() instead. Could there be an internal library routine that is stable that making a builtin mean() makes sense?
Still, here for median() it is only two numbers we are averaging, plus the formula in the reference pertains to float values, not integer "wrap" effects. Consider that the formula from the reference becomes
which is the same formula as in Comment #4.
KT, keep thinking about the formulas; it shouldn't be too difficult to implement conditionals with scripts. They slow things down, but I'd prefer being correct than have a warning message. For example, the formula of Comment #4 could be implemented as
retval = min (x, dim) + (max (x, dim)  min (x, dim))/2;
Maybe there is a way to generalize the formula from the reference for integers. The issue with integers is the loss of the fractional portion when doing the intermediate mean in the recursion. Perhaps there is a way to accumulate those lost bits as some sort of residual and then at the very end of the computation adjust the final computation with the residual divided by N.

Numerically stable mean value algorithms have been around
for a long time. The one shown in the URL below dates from
the 1960's:
https://diego.assencio.com/?index=c34d06f4f4de2375658ed41f70177d59
The original references should be easy to find.

Good catch with the (u)int64 data type.
I think the second version of comment #2 will deal with averaging two numbers for the even case with (u)int832, single, double, and logical. For (u)int64, we might have to introduce an exception
The formula of comment #3 and comment #4 is hard to implement regarding multidimensional arrays. If you found a more elegant way, I would be very happy ;)

"For uint64, it's just meanXY = X + (YX)/2."
Actually, need to ensure that we pick Y >= X or there is wraparound that way.

This is one of those I'd classify as "think about for a while". I did notice the "logical" case, both within median() and mean(). However, I din't think why it is there, just that mean() with the followup check within median() should still work fine. Of course, in the alternative changeset there is an extra call as well, to cast().
Yes, mean() is an extra function call of overhead, as mean() in turn calls sum(). For large data this overhead should be diminished. There's nothing recursive or looplike about these routines, so the greater the input matrix size, the more CPU goes toward the inner matrix operations written in C++.
I'm fine with dropping the "native" and then use cast(). I was thinking more along the lines of there being just one location where the sophistication of mean() is implemented. The one thing to be careful of, though, is something we've run into with variable indexing, i.e., a double has 52 bit mantissa or significand. Hence, there is this 12bit range where int64 numbers lose resolution when computed with double values.
The general case of mean() means this is complicated. We are summing together unknown N numbers, so we need ceil(log2(N)) bits overhead to ensure there is no arithmetic overflow. But in the case of median() it is just two numbers that involve mean. So, should we consider the special cases of median of int64 and uint64?
For int64, given two numbers X and Y, if X and Y are of different sign, then meanXY = (X+Y)/2. If X and Y are the same sign, then meanXY = (2X  X + Y)/2 = X + (YX)/2, or something close ensuring rounding is appropriate.
For uint64, it's just meanXY = X + (YX)/2.
Keep in mind that this needs to be for vector/matrix data, so the test above would involve logical arrays in some way.

Matlab R2018a returns "140" in both cases of the original comment.
The patch file #44857 of comment #1 works
But I think it makes things more complicated and maybe slower for large data?! Basically, mean() would perform the computation with double precision and finally convert the result to the input data type. Why not doing the very same without additional function calls in one line and tidy up the nasty "logical" exception:
I thought about other solutions from the "real plain dump" integer world (like C/C++) but these would not apply for Octave... This I am afraid to stay to this solution.
What do you think about this Dan Sebald and Etienne de Foras?

The question is whether Matlab does the arithmetic using the intrinsic variable type. The median() routine is supposed to return the same variable type, but that doesn't necessarily mean the average of the center values should be computed using that type. Why should it, I guess? In fact, mean() appears to do its arithmetic not using native, but instead converts it afterward:
So, attached a patch to change the behavior for median. Basically it uses mean() instead of sum() with the assumption that the proper averaging behavior is already baked into mean(). I added a couple more tests as well.
(file #44857)

If I define:
u= [1 140 140 255];
u8=uint8(u);
then median(u)=140 (correct)
but median(u8)=128 (incorrect, should be 140)
The bug disappear with odd size data, the median with even size data is computed as (140+140)/2, but with saturation to 255 due to uint8 type and then divided by two >127 (or 128 with rounding).
We may use the type int32 or double for the sum and division by 2 and then convert back to uint8
(Thanks for GNU Octave !)
