Here is a version that attempts to properly account for the difference in probability for the last "hi" bin. Admittedly, this probability is very small, so I don't expect much practical impact.
function R = randi64(maxval,sz)
%RANDI64 Generate random numbers between 1 and a given 'uint64' value.
% Total number of random numbers desired
N = prod(sz);
% We work on generating values between 0 and maxval-1, and then add 1 at
% the end to make the indices one-based.
zmaxval = maxval - 1;
if zmaxval < uint64(2^52-1)
% By default, randi generates doubles with mantissas of 52 bits, so we
% should be able to safely generate all integers between zero and
% 2^52-1, and then cast these to 64-bit integers.
R = randi([0 zmaxval], N, 1);
R = uint64(R);
% The logic here is that we split the random number into two 32-bit
% parts, referred to as "hi" and "lo".
nbits = 32;
zmaxval_hi = idivide(zmaxval,uint64(2^nbits),'floor'); % number of hi bins
zmaxval_lo1 = mod(zmaxval,2^nbits); % number of lo bins for last hi bin
zmaxval_lo0 = 2^nbits-1; % number of lo bins for every hi bin but the last
% We need to handle the last high bin specially because its probability
% is not equal to the other bins. We figure out its probability and
% separate the random numbers into those that are in the last bin (lo1)
% and those that aren't (lo0).
p_last_hi = (zmaxval_lo1 + 1) / maxval;
last_bin = rand(N,1) <= p_last_hi; % True if in last bin
R_hi = zeros(N,1);
R_lo = zeros(N,1);
% First, let's deal with the last "hi" bin. The "lo" values in this bin
% don't go all the way up to 2^32. It's generally a very small
% probability of being in the very last bin!
n1 = sum(last_bin);
if n1 > 0
R_hi(last_bin) = zmaxval_hi;
R_lo(last_bin) = randi([0,zmaxval_lo1],n1,1);
% The rest is generated normally, i.e., the hi bin is between 0 and the
% next to last bin, and the lo bin can take all 2^32 values.
n0 = N-n1;
if n0 > 0
R_hi(~last_bin) = randi([0,zmaxval_hi-1],n0,1);
R_lo(~last_bin) = randi([0,zmaxval_lo0],n0,1);
% Now combine the hi and lo bits.
R = uint64(R_hi)*(2^nbits) + uint64(R_lo);
R = 1+R;
R = reshape(R,[sz 1 1]);