I reviewed and amended Michael's patch and pushed it here (https://hg.savannah.gnu.org/hgweb/octave/rev/25d3e8e49d5c). This will be a part of the 5.0 release in early January. Marking as fixed and closing report.

@Michael: All fair comments. I saw that the loop would be executed infrequently, but didn't continue to calculate how infrequently that might be.
The worst case is unbounded. But, what I had in mind was column 3 of Table 1, "Maximal number of remainders per integer" which is Infinity for the Java algorithm which is why I wasn't interested in it.
Total time = time_to_create_random_integer + time_map_integer_onto_range.
Running the openbsd algorithm in randi_rej.cpp with
If I change the code to
, recompile, and rerun the 20 point benchmark I get
So the random number generation is about 43% of the total time and the mapping is 57% of the time. Since it is so balanced, I think there are opportunities for improvement in both halves. If I run the nearlydivisionless algorithm on a range where it works, imax=9, the results are impressive.
But, since it didn't work for all values of imax I can't advocate for it.

Rik, I have to comment on a number of points in your series of comments.
First, I would say that there are just two main variants of a single algorithm for deriving unbiased random integers in an arbitrary range from a source of unbiased random integers in a given range. Both use rejection: For instance, when you have a generator for integers in [0,65535] (i.e., 16 bit) and you want to generate integers in [0,9], you first generate a source integer, which you reject and regenerate iteratively if it should fall within [65530,65535]. So now you have an unbiased integer r within [0,65529], and you get your output by doing either mod(r,10) or floor(r/6553) (when you are doing integer arithmetic, you can drop the floor). What is discussed in the arxivmanuscript is how you can do that while keeping the number of integer divisions low. So actually all the implementations have unbounded worstcase behaviour, because they all use rejection.
I think you misunderstood my solution. Yes, there is a loop, but it does not loop over the elements. Instead, I generate first more numbers that I expect to need, do everything vectorized, and only when at the end I see that the number of generated numbers is not large enough because I had to reject too many, I try again. If have the comment "should practically always be true" in the code, because I generate use an excess number that is 10 times the expected standard deviation of the number of rejections. In octave, normcdf(10) evaluates to exactly 1, while 1normcdf(8) is 1e15. Thus the probability for the loop having to be evaluated a second time is probably about 1e20. That's what I mean by "practically always".
Yes, a compiled function will always be faster than a vectorized mfile, mainly because of the fact that with vectorization, you have to repeatedly retrieve and store back your intermediate vectors to main memory, while in a loop you have only to store the final value, everything else happens in registers. Further, when you directly use randi53(), you obviate the cast to double and subsequent division that happen in rand(), and you do not have to undo the division. I think that these are the main reasons why your code is faster, not so much that the arxivalgorithms are better (because they are actually the same).
So how fast is my code? I see
About a factor of 2.2 slower than the simple biased mfile when one in three values has to be rejected, dropping to a factor 1.6 when practically nothing is rejected. The slowingdown with rejection cannot be avoided, and a factor of 1.6 for an mfile is not bad, I think. So as long as it is about addressing the bias, and as long as nobody finds a bug in my implementation, I still think it can be used.
Thinking about what to do when you really want to beat Matlab: I think you meant that rand should keep returning doubles within [0,1). Yes, that goes without saying. For this you need randi53. But what would be nice to expose would be something that indeed returns a uint64 with the full range of [0,2^641], because this is in most likelihood what the Mersenne twister internally returns. And further, for beating Matlab with something "new" I would propose to change the underlying PRNG to something of the PCG family http://www.pcgrandom.org/, which are both among the fastest and "randomest" PRNGs available today. I see that Melissa O'Neill even has a very recent post about just the problem we are discussing here: http://www.pcgrandom.org/posts/boundedrands.html

Going forward, I have other bug reports that I want to sqaush so I won't be working on this. My suggestion would be implement the functionality in liboctave in octrand.cc or a combination of octrand.cc and randmtzig.cc. The interpreter then needs a way to access the new functionality. I would probably put that in libinterp/corefcn/rand.cc. There is code in rand.cc, or fill_matrix() in data.cc, for extracting the dimensions from a variable number of input arguments which could be borrowed for randi input processing.
For the moment I would prefer to keep rand returning values within [0,2^53). The base data type in Octave is a double with only 53 bits of precision and I fear introducing more bias during the conversion from 64 to 53 bits. Also, minimal change necessary to solve the problem seems best from a software engineering perspective.
There is still a good GSOC project here to convert to using C++ standard libraries for random numbers. Maybe simultaneously with that project we could investigate the impact of returning a uint64_t from an internal API.

It would be fun to have bragging rights to an algorithm which is both faster than Matlab and unbiased so I tried implementing the "new" algorithm in the paper in randi_new.cpp. The author outlined the algorithm, and also provided a code implementation for L=64. Because we have L=53 I could not use his code directly and instead coded the algorithm as I understood it. What I found was that it worked for small values of imax, but hung in an essentially infinite loop for large values. That may be my misunderstanding, but I wasn't super motivated to find out. The algorithm relies on 128bit unsigned data types, which g++ does support, but they are not part of the C++ standard and so there would need to be some extra work in the configure script to determine which compiler was being used, whether 128bit primitives were supported, and a fallback mechanism if they were not.

For testing, I used a script tst_randi.m which is shown below, and attached.
It generates one million values using a particular method and then uses Michael's test even/odd test. It prints out both the absolute counts and the percentages (which should be 50%:50%).
As a warmup, running tst_randi with the existing algorithm (randi) produces
So, there is a clear bias in the results using the existing algorithm. For some performance benchmarkings,
So, roughly, 30 milliseconds. Using the rejection method,
Very roughly, 45 milliseconds, or about 50% slower. Quoting from the paper,
The worst case is when 2^L mod s is large which is at (2^L / 2) + 1. For our case, this is (2^53 / 2) + 1 = 4503599627370497. Testing this idea proves to be correct.
Maybe this averages 80 milliseconds. But of course, the real reason to use thismethod is because it is unbiased. Running tst_randi shows
Also, it is unlikely that users are choosing pathological values for imax frequently.
(file #44968)

I did some work on this problem over the weekend just to hack something up. First, I'm attaching the academic paper that Michael Godfrey pointed to (Fast_Random_Integer_Generation_in_an_Interval.pdf). This outlines two methods used by existing libraries: 1) rejection (OpenBSD, GNU libc, others), 2) Java. Of these, rejection is used by many and it has guaranteed worst case behavior. The Java library has better average performance, but unbounded worst case behavior, which I'm not interested in adopting. Finally, the author proposes a third method of his own devising which I will call "new" method. The paper is from May, 2018 so this is very "new" compared to the existing methods.
I appreciate Michael Leitner's mfile approach, but it does contain a dountil loop. Because Octave is interpreted, and we don't have a JIT implementation for loops, any looping behavior involves a big performance hit. In situations where looping can't be avoided the suggestion is to move to C++ which is what I have done.
Second attachment is a small change to expose the function randi53 (53bit random integer) from liboctave/numeric/randmtzig.cc. See randi53.diff. Because of some vagaries of the build system to use this you will need to patch the Octave source tree and then install Octave so that 'mkoctfile' has the new header file randmtzig.h. Instructions
After that you can work with the two octfile samples I made: randi_rej.cpp and randi_new.cpp. These are not completethere is no input validation, documentation, support for specifying number of values, BIST tests, etc.but they do allow you to set imin and imax and they return a column vector [1e6,1] of values.
To compile them use
More to follow, but I can only attach 4 files at a time so I have to continue in a new comment.
(file #44964, file #44965, file #44966, file #44967)

See the attached patch (a diff against the current randi.m). I think it is correct, at least it passes the test on the ratio of evens and odds. This could be also added in the tests, but I do not know how to write a test on stochastic functions. Of course, that it does pass this specific test does not prove anything, so I invite also other people to have a look on it. Essentially it just tries to create integers between 0 and flintmax1 (2^531 for IEEE754 doubles), divides them by another integer K so that the requested range is covered (for exact integer arithmetic this should result in exactly uniform probabilities within the requested range) and rejects all above the range. It initially generates a buffer of random numbers that is in all likelihood larger than what is requested, so that it is only insignificantly slower than the previous implementation for cases where not many numbers have to be dropped (for small ranges, that is, where the previous implementation was already nearly correct). In these cases, the sequence of returned numbers is also the same, provided not too many numbers are requested.
A potential issue: is IEEE754 double division predicted to be exact on exactly representable integers, when the result is again an integer? If division is multiplication by the inverse, then probably not, because the inverse is in general not exactly representable. In this case, it would probably be better to do these manipulations in uint64, which should be required to be exact.
As a side note: I would be very much in favour of exposing the fundamental random number generator not via rand(), which returns doubles in the range [0,1], but by a function returning uniformly distributed uint64, which they internally indeed are. This would save operations in the present problem which have only the purpose to undo this conversion to scaled doubles, further you can get the scaled doubles by a twoline mfile, and finally this would then document what rand() really returns (is exact 0 and/or exact 1 included, if so, with which probability, and is indeed every double between 0 and 1 returned or only those that are rounded uint64 divided by 2^64 and so on).
(file #44961)

Thanks Michael for the observable example. I've changed the release field to dev. This should be fixed for the 5.0 release due at the end of this year.

dev Octave (~4 sept):
...and Matlab r2018b prerelease:

I think the test case should be changed to:
Thanks, rahnema1

So at least it's a Matlabcompatible bug

Matlab:

Thanks for illuminating!
Octave
Julia seems to get it right
Note: there seems to be abug in fit which fits 3 bins instead of 2
Numpy aslo gets it right
From my view this is now a serious bug, we should fix it asap.

Under normal operation, the bias is really tiny, so that you will never be able to see it with such tests. The deviations you see are just what you expect for a binomial distribution. But as we know how the code is working, we can easily construct a pathological problem:
This just generates a bunch of random integers, takes the larger half of them, and computes how many of them are even or odd. The result should be that on average exactly half of them are odd, but here you see that it is rather two thirds. So this is a significant deviation, and when you blindly assume that randi is doing its job correctly, you could generate capital nonsense.
This works because the magic number 6004799503160661 is bitmax()*2/3, so when we assume that rand() itself works correctly, that is, that every possible double between 0 and 1 appears with the correct probability (not with equal probability, because there is exactly the same number of floatingpoint numbers between 0.25 and 0.5 as between 0.5 and 1), then in the upper half of numbers either one or two doubles get assigned to one and the same integer.
It would be interesting to see how Matlab performs here, or Julia and numpy. If most of the competition does it correctly, then I would just quickly hack up the rejection part, which would not be too difficult.

The code below to try to quantify the bias is pointless, sorry. Verifying correctness of a PRNG is a very complicated thing and nobody knows how to do it right. You need to look at this issue from first principles: compute the probability for each output value.
Here's a quick explanation of the problem with a extreme example:
Assume you have a 2bit RNG (generates values in the range [0,3]. You want to generate integers in the range [1,3]. You use Octave's method:
You are now mapping 4 possible values into 3 bins. 2 of these bins will have a 25% chance of being hit, 1 bin will have a 50% chance. You can do the mapping however else you want, the outcome will always have this bias. The only solution is to reject 1 in 4 results of the RNG.
With a 53bit RNG as Octave uses, this effect is much, much smaller, but it is still there.
I hope this example clarifies the problem.
I suggested earlier to "throwing away randmtzig.cc and replacing that with stuff from the Standard Library", but what I really meant was to replace the `randi` Mfile with a call to the Standard Library. The random number generator you use is OK, and it might even be possible to use `std::uniform_int_distribution` with it. I would recommend replacing that also, just because it's less code to maintain, but if you need to keep the old behaviour intact, maybe make the fewest changes that fix `randi`.

Rik,
Sounds good to me. This is a good Summer project task for a start.
The C++ code should in any case be selected based on an option
setting. As you said, the old code will need to be supported for
compatibility reasons. That also means that the sooner the
old code is frozen the better. This also supports (my) view
that "fixing" the bias in the current code is not a good idea.

@Michael: I think this might make a good Google Summer of Code project. This task isn't large enough to be a standalone project, but it could be part of a larger project to change over the core to use the C++ libraries for random number generation. I read through the paper you sent and I was ambivalent. It's quite nicely worked out, and there are code samples, but if we move in that direction then we also are back in the business of maintaining a code base for random numbers. There is a simplicity, and reduced maintenance burden, to just calling library code and letting the library maintainer's worry about such things. Of course, a true engineering decision would look at a number of variables like accuracy, performance, maintenance burden, etc. That's why I think it might make a good GSOC project for someone to explore these concerns and then implement a solution after consultation with the Octave Maintainers.

Rik,
Discussion of "correctness" of random number generators
has a long history and has not been particularly helpful.
The idea of a "perfect" random number generator appears
to be wishful thinking.
However, using the C++ implementation for this case seems
a reasonable choice combined with the option of backward
compatibility as you point out.

Rik,
This deserves a careful look. The paper at:
https://arxiv.org/pdf/1805.10941.pdf
may also be helpful.

For the sake of comparisson I run Rik's code in Julia 1.0.0, I get thefollowing deviations
max 0.04% in bin 2
and a second time
max 0.056% in bin 5
In numpy 1.14.2 I get
max 0.076% at bin 7
and again
max 0.08% at bin 9
Both are similar to Octave's deviations.
I appreciate the comments from the OP. If he can just illustrate the value of his suggestion with a working example, then I am pretty somebody will be happy to improve the method.
It seems there is no agreement on what "serious scientific software" is, at least Julia and numpy are considered pretty serious.
Julia code:
Numpy code:

I hear the assertion that the Octave method is incorrect, but I haven't been able to determine that myself. But, assuming it is incorrect, is it inadequate? That's why I want a numerical estimate of the bias of the current system. The last argument, not reinventing the wheel, is the most compelling. Octave only recently began requiring C++11 so it is now possible to contemplate using the STL for this case. Long term, we would like to move all of the random number generation to the standard library, but we also have to support legacy random number generators (as Matlab does), so the code is always going to be rather complicated so that programmers can reproduce simulations they have written years ago.

There is a correct way to generate random integers in a given range, and there are incorrect ways. The method that I linked is the correct way. Other ways are incorrect. The bias might be small, but it's there.
GNU's C++ Standard Library does it right also:
No serious scientific software would use the method used by Octave. This is why I'm reporting it as a bug. Take it or leave it.
Note that random floatingpoint numbers are totally different from random integers, and should not be used as a comparison. And they should not be used as a basis either.
This "large project" is just throwing away randmtzig.cc and replacing that with stuff from the Standard Library. Random number generation is complicated enough that one should not attempt to reinvent the wheel.

the rand from Octave is implemented in liboctave/numeric/randmtzig.cc. As far as I can tell, there is no 64bit integer that is cast to double. Similarly, there is no "division + floor" in randi. There is multiplication, and floor is used to remove the fractional part. Because I can't verify some of the first assertions, I can't go on to accept the conclusions that follow. The Stack Overflow response has this comment
So it's not clear, even to the original author, that Octave's current system doesn't work well enough. Before any large project is undertaken it makes sense to have an idea of the costs and the benefits. It seems to me that the benefits need to be quantified. How large is the bias of the existing system (if it even exists) that is implemented in Octave?
Also, relevant, are we any worse than the competition? According to The Mathworks, the algorithm Octave implements is the one suggested for random double numbers (http://www.mathworks.com/help/matlab/ref/rand.html)
For random integers, The Mathworks suggestion is to use randi(), but I don't know if they implement something special there or not.
As an attempt at quantifying the bias, I used
So, as much as .05% deviation, but this could be statistical fluctuations. Running a second time I get
This time it is a .06% deviation, but it is not in the same bin as before so it doesn't look like systematic bias.

This has nothing to do with the C `rand`. Octave's `randi` is an Mfile, you can type `edit randi` to see it. It uses Octave's `rand`. Line 106:
`rand` uses the Mersenne Twister algorithm, which outputs a 64bit integer, cast to a double. Even if this cast is done correctly, leading to a perfectly uniform distribution of values, this double float value doesn't sample integer values equidistantly, and the output range likely doesn't divide the input range evenly either. Thus, the division + floor in `randi` will lead to a bias in the output distribution.
In short, a correct `randi` is implemented using the integer output of the Mersenne Twister, and using the rejection scheme as shown in the linked Stack Overflow answer.

Could you point to the code in Octave where this division occurs? Octave is not using rand() or random() from the stdlib for its implementation of rand() at the interpreter level. Thus, I'm not sure the analysis applies.

The randi function uses an incorrect algorithm that leads to biased results. It is implemented in terms of `rand`. But `rand` is created by dividing the output of a 64bit random number generator by 2^641, scaling this value up and rounding down doesn't lead to a uniform distribution, as the original integer intervals are most likely not evenly divisible into the output interval, and thus some output values are slightly more likely than other output values.
This is maybe only a small bias, but it is there.
See here for a correct algorithm: https://stackoverflow.com/questions/2509679/howtogeneratearandomintegernumberfromwithinarange/6852396#6852396
