Thanks to all for their investigation, and also for reporting this!

This is indeed a fascinating problem. The reason why the higher-dim output ndarrays have dims that start with 1, then are powers of the original size, is because the pdl struct has a def_dims array, size 6. Straight after that is the def_dimincs (stride-size at each dimension), also size 6 (the API has code to allocate memory if >6 dims, clearly that's not getting called here). Both are being used, even though there are more than 6 dims. Reading the array of dims will therefore give data that starts out correct, then looks crackers. It's almost certainly a bug in the rangeb code, since higher-dim ndarrays surely get tested elsewhere and the API ought to be solid.

Further discussion etc on the GitHub issue, but I thought I'd give an insight here into what's going on.