Replace generation of random reals in `reservoir_sampling_aexpj` with non-biased PRNG

When MR !2528 (merged) is merged, implement this function:

!> @brief Return a PRN in the range [a, b].
function random_number_in_range(a, b, seed) result(scaled_rn)
   real(real64), intent(in) :: a, b  !< Lower and upper limits
   integer, intent(inout) :: seed
   real(real64) :: scaled_rn

   integer(int64) :: rn_int
   
   rn_int = xorshift64(seed)
   ! CHECK: Assume max value returned by to_double() is (1.0_real64 - M_EPSILON)
   scaled_rn = a + (b - a) * to_double(rn_int) / (1.0_real64 - M_EPSILON)

end function

and replace the calls in quickrnd_oct_m/reservoir_sampling_aexpj with random_number_in_range