Introduction to Algorithms (CLRS), Excises 5.1-2

CLRS_5_1-2.jl

#=
Introduction to Algorithms (CLRS), Excises 5.1-2

Problem Definition
- Describe an implementation of the procedure RANDOM(a,b) that only makes calls to RANDOM(0,1).
- What is the expected running time of your procedure, as a function of a and b?
- RANDOM(0, 1) returns either 0 or 1 with equal probability.


I came up with a lovely but not efficient solution.

At first glance, a naive approach is to halve candidates recursively by RANDOM(0, 1); however, this is wrong.
https://walkccc.github.io/CLRS/Chap05/5.1/#51-2-star

The counterexamples are when b-a+1 is odd.
https://github.com/walkccc/CLRS/issues/253

One of the correct solutions is to construct a bit array using RANDOM(0, 1).
https://sites.math.rutgers.edu/~ajl213/CLRS/Ch5.pdf

I reconsider the first wrong approach.
The problem is that the probabilities are not kept correctly with odd candidates.
My solution is doubling the list of candidates, e.g.;
[1, 2, 3] --> [ 1, 1, 2, 2, 3, 3 ]
We now can use the first approach while overcoming the problem!
The expected running time seems to be O((b-a)lg(b-a)).
I roughly explain as;
- at each iteration, candidates are almost halved -> O(lg(b-a)) iterations in total
- each iteration requires O(b-a), due to the doubling operation
I did not check the details, but you know, it is not efficient :D
I could not find a similar solution to the problem 5.1-2 on the web.
Any comments are welcome! (e.g., counterexamples, formal proof, etc)


Note that this is my first julia-lang experience, it may not be sophisticated.
=#


using Random

function random_ab(a::Int, b::Int)
    function f(lst::Array{Int})
        # finish
        if lst[begin] == lst[end]
            return first(lst)
        end
        half = length(lst)
        # doubling candidates, e.g., [1, 2, 3] -> [1, 1, 2, 2, 3, 3]
        doubled_lst = reshape(permutedims(hcat(lst, lst)), half*2)
        if rand((0, 1)) == 0  # RANDOM(0, 1)
            return f(doubled_lst[begin:half])
        else
            return f(doubled_lst[half+1:end])
        end
    end

    a <= b ? f(Array(a:b)) : error("a:$a should be <= b:$b")
end


# simple experiment
a = 1
b = 3
repeat = 1000
cnts = zeros(Int, b - a + 1)
for _ in 1:repeat
    res = random_ab(a, b)
    cnts[res - a + 1] += 1
end

for i in a:b
    freq = cnts[i-a+1] / repeat
    println("num=$i, freq=$(freq)")
end