@jwe: regarding comment #9, do you know of a fancy way to use templates to reduce the number of files that instantiate operators? Or, is it just the case that every file in libinterp/operators that defines operators between two data types needs to be modified to generate a warning?

I checked in a major rewrite of the indexing section based on your work and my own. See http://hg.savannah.gnu.org/hgweb/octave/rev/3d5f953e2ef6.

I rewrote large parts of 8.1 and 8.1.1 (and a very small modification to 8.2). I removed the unnecessarily complicated stuff and tried to write the rest both more accurately, thus the level of complexity may seem about the same. Further, I gave more basic examples. The diff is attached.
Further, I think that the documentation would gain by deleting the last part of 8.1.1 starting from "A second use of indexing...". First, I do not see the welldefined first use of indexing, but then the given examples are either not about indexing at all (the point with sqrt, which is rather about using vectors and matrices as fundamental data types rather than scalars, which is probably the reason why anybody would turn to Octave in the first case, or in any case the first thing he would learn  discussed already in 19.1), or should rather be done by a repmat in the second example. Preallocation is already done in 19.6. I do not see that it is specific to indexing.
Last point: From inspecting this part of the documentation, I did not really get the conventions with @samp, @math, @var, @code. I did it as good as I could, but feel free to make these more consistent if there are conventions regarding their use.
(file #43982)

The manual chapter for this is in the file doc/interpreter/expr.txi if you navigate from the top of the Mercurial repository.
Mercurial changesets are the best submissions. There is some documentation on the Wiki here: https://wiki.octave.org/Mercurial.
Otherwise, a unified diff will work (diff u).

I tried adding a warning message for double arrays multiplied by double arrays only. The patch is attached as mm.diff. It seems to work, but as I understand it, the full solution would require modifying every file in libinterp/operators that involves arraybyarray operators. That is a lot of files because there seems to be one file for each twoway combination of arrays that can be operated on. For example, there is float arraybydouble array, complex double arraybypermutation matrix, diagonalmatrixbydoublematrix, etc., etc. That is a lot of files, unless there is a easier way to do this.
Octave is already using a bunch of macros in these files. Maybe clever use of templates could help? Or maybe I have to buckle down and write a Perl script to do all the modification.
(file #43974)

Aha. I did not find that documented anywhere on the Mathworks page, they just present the cases of conventional and socalled linear indexing, which is the special case of indexing into an Ndimensional object by just a single dimension.
By the way, what is the preferred way to propose patches for the online documentation? Is the source for the documentation to be downloaded somewhere?

In Matlab R2016a:
So Octave is behaving compatible to Matlab with that respect.

"and is thus unnecessarily complex". Yes, the language is really difficult to decipher. It was written ages ago and really not updated. Anything that could be done to improve it would be a help.
As for whether this is an extension over Matlab, I'm not sure. I don't believe it is, but someone could test with this code
The original array has only 2 columns. Asking for column 3 should either provoke an error from Matlab, or return 5 as Octave does.

With respect to the behaviour of mtimes and *: The last paragraph in my initial bug report was nonsense: as you say, what I describe can be done with the double colon trick without any explizit reshaping to singleton dimension and so on. So I do not see a use case for the present behaviour (apart from saving 5 characters for the parentheses, colons and the comma), and thus it is probably best to do input validation and raise an error when any argument has more than two dimensions.
As regards the documentation: yes, in 8.1.1 there is some documentation, but specifically the action of the colon in this respect is at best implicitly defined. If I find time, I will think about a rewording of this section.
But: isn't it the case that the second paragraph in 8.1.1 that you reproduced below treats cases that cannot validly appear and is thus unnecessarily complex? What I mean is the case "nidx < nd but one of the indices idx(i) is outside...". If i<length(nidx), that is, if this excessive idxentry concerns any but the last specified dimension, it errors out in any case? So to me it seems that we can define the action described in this paragraph as: "If nidx < nd, then the result is an implicit reshape that folds the last nd  nidx + 1 dimensions into one." This includes the case that idx(end) is within the size of the corresponding dimension, in which case it would be equivalent to the present definition via addition of singleton dimensions (meaning: idx=[idx ones(1,ndnidx)]), as well as if the last entry is larger than the corresponding dimension. If any other entry is larger than its dimension, it would give an error.
By the way, this is presumably again an extension over what Matlab allows, I would guess Matlab would allow only nidx==nd or nidx==1, in which case the answer is A(:)(idx). This should then probably also be documented.

From the test in comment #3, and because I think it is good practice, Octave should implement some sort of input validation for matrix multiplication. Should it be just a warning, or a full error as Matlab does? I've added jwe to the CC list.
I'm leaning towards a full error, as any programmer who really, really wants to condense an ND matrix to a 2D matrix can do so themselves with indexing:

Is something wrong with Savannah? Michael Godfrey posted a response to this report, which I will include manually, but it doesn't seem to have been received to this bug report.

Do you have access to Matlab? Does Matlab produce a warning for any of the code below?

We certainly don't want to get rid of it. As for documentation, see the section 8.1.1 "Advanced Indexing". There is documentation there, but it is not an easy concept to explain. The text is
In this situation, it is the last bit which takes precedence. By indexing a 3D object with only 2 indices (:,:) all of the dimensions >= 2 are folded in to the second dimension. So rand (3,3,3) becomes a 3x9 2D matrix.
If you can think of language that would be clearer to explain this it would be useful.

I do not know whether I am searching incorrectly, but it seems to me that all that the documentation has to say about mtimes is "Return the matrix multiplication product of inputs." But what if the arguments are not matrices, that is, have more than two dimensions? By testing it seems to me that
is exactly equivalent (in terms of results, but also under which conditions it errors out) to
(I do not know whether or where this syntax is documented), which in turn is equivalent to
Even the error message is the same: for
all three give "error: operator *: nonconformant arguments (op1 is 1x6, op2 is 3x9)".
I have been using Octave now eight years nearly daily and did not know about that feature, so this should be documented at least in section "8.3 Arithmetic operators", but perhaps one could add a subsection under "8 Expressions" called something like "Advanced trickery", where such things (also including indexing by two or more colons) is presented. Further, perhaps the error checking should be done a bit earlier, so that in the abovementioned example the actual dimensions of the input arguments are reported. And finally, perhaps it would be good to emit a warning whenever (length(size(a))>2)(length(size(b))>2), because this is probably often unintended.
However, I do not want to argue that under the abovementioned condition this should raise an error, as the behaviour is actually useful: for instance, it allows me to have vectors organized in more dimensions, and to transform them all by a simple matrix multiplication from the left with a subsequent reshape, as opposed to have to reshape each argument and introduce singleton dimensions, do an elementwise multiplication with broadcasting, sum over the relevant dimension, and do a final reshape.
