Thank you, Guillaume, for applying the principle to the helper function. Actually, last night I was also active. I attach a version of _unite_shared_vertices_.m that is now perhaps a bit more efficient because the faces are updated only once, and I have now a tolerance criterion that should be the tightest possible (it is now again vectorvalued). And further, it is stable in the sense that the order of vertices as they are returned is the same as they were provided. This was not the case with my previous version, but it was the case with Markus' initial version, and this was a condition for his deletion of facevertexcdata to work. If you add the lines at the end
then you should be able to use this as a dropin replacement for the previous version. More efficient would be to return J as it is and select the corresponding elements from facevertexcdata (then it would be more flexible if in the future this function is used for other purposes), or as you did it to pass it in and out and select from it in the loop.
Rik, sorry, I am doing this here on Octave 4.0.3 (greetings from Debian Stable) by pulling in the relevant .mfiles from the web repository. I do not have a full mercurial setup here so that I could prepare a changeset.
(file #46213)

I attach a version of Michael's patch where the changes now take place in _unite_shared_vertices_.m. As discussed here, the tolerance is now different than before and I don't know what is best; it passes all tests.
(file #46212)

@Guillaume or Michael: Can you work the patch in to a complete changeset? The speedup is very impressive.

Answers to your comments:
 yes, then it would be best to apply my patch just to _unite_shared_vertices_.m (and compute J according to its previous definition), because it indeed is the vectorization Markus asked for in the FIXME. Thus the speed improvement would apply also to reducepatch.
 I was initially thinking that marching_cubes works sequentially, thus that every cube can affect every other one. Then the eps of the largest vertex overall would indeed apply. But yes, seemingly it works in parallel, so only the largest vertex in the cube defines the achievable precision. So yes, it would become a vector, but I think the previous choice was too optimistic (when one vertex results very close to zero, its position is still affected by the finite precision of the other vertices of the cube it resulted from). But there is no one true solution  you will always either keep returning vertices as distinct that with perfect arithmetic would be just one, or you will merge vertices that actually should be separated. So I would say, it is good enough.
 You want to drop all faces that have at least two equal vertices. And if neither f(1)==f(2) nor f(2)==f(3) nor f(3)==f(1), then necessarily all three are distinct. If you knew that the faces were sorted, then you would need only two comparisons (as in the last lines of _unite_shared_vertices_.m).

Thanks Michael, the speed improvement is impressive! To be fair, the speed issue was acknowledged by the authors in the FIXME note in _unite_shared_vertices_.m and it was probably not so important for small examples.
A few comments/questions:
 _unite_shared_vertices_.m probably still has to be kept as it is also used by reducepatch.m (patch #8912).
 Concerning tolerance, this was discussed in bug #46946, see Marco's comment #8. This would make the tolerance a vector instead of a scalar in the loop?
 The one line I don't understand is this one. Why is it enough to only compare with a [2 3 1] reordering?

This is a great opportunity to use the builtin profiler.
If I then use profexplore I see
So, Michael is correct that _marching_cube_ is way more efficient than _unite_shared_vertices_.
Michael's patch definitely decreases the time required. This isn't my area of expertise so someone else should evaluate the mathematical logic (5 eps, for example as the tolerance).

Yes, it is embarrassing if the main problem is handled by an efficient algorithm like marching cubes, and the postprocessing is extremely inefficient  to give numbers, on my computer the four lines with n = 64 take 15 seconds, while adding "noshare" to the invocation of isosurface cuts this down to 0.08 sec.
Your idea is obviously correct. Below I attach a patch that incorporates this idea. I also reworked the actual deletion of the vertices, which are now only seven lines (excluding comments) and thus probably do not warrant an outsourcing to _unite_shared_vertices_.m any more.
The example runs now in 0.10 sec. This could perhaps be decreased still by some 10 percent, along the lines of the FIXME in the code. But no need for C++.
(file #46178)

isosurface becomes very slow when the size of the input data increases:
Octave takes 11.34s while Matlab takes 0.04s.
If you increase n to 128, Octave takes 200s and Matlab 0.25s
And for n = 256, Matlab takes 2s while Octave was still running after an hour.
Most of the time is spent in _unite_shared_vertices_ where the loop over vertices is known to be slow (see bug #46946 and patch #8912).
One way would be to move _unite_shared_vertices_ into C++. Otherwise I wonder if one can implement 'unique with tolerance' from sorting the vertices:
I also notice that the help text of isosurface still mentions the use of 'unique':
