What does it mean to sample from a hypergeometric distribution? According to https://crypto.stackexchange.com/a/8800/53007:

Start with the entire domain [M] and range [N]. Call y←N/2 our range
gap. Now using our key k we generate some pseudorandom coins and give
them to our HGD sampling routine along with y, M, and N. This gives us
an x≤y that describes the number of points of our order-preserving
function less than y.

What does the HGD sampling does exactly? It gives me a "good" or a "bad" ball at random? How can this give me any information about the number of points less than y?
For me sampling from an interval according to a distribution would give me a random number in any distribution possible. I don't see how it gives me the information cited there.
What does it mean to sample from the hypergeometric distribution HGD?
 A: When people talk of "sampling from a distribution", they usually refer to obtaining a value of an appropriately distributed random variable.
A (real) random variable assigns a number to each outcome of a sample space. An example is the function that assigns $0$ to the outcome "heads" and 1 to the outcome "tails" in the random experiment of flipping a fair coin. This variable follows a uniform distribution in $\{0,1\}$: each value happens with probability $1/2$. If you want to "sample from this distribution", you can just flip your fair coin and write down the number corresponding to the observed outcome.
Imagine you now want to sample from a different distribution: $0$ happens with probability $1/4$ and $1$ happens with probability $3/4$. You could set up an urn with 1 red ball and 3 blue balls, and pick one at random. If you don't have this equipment nearby, what you can do is flip your coin twice. If you get two heads, you write $0$. Otherwise, you write $1$. You are now sampling from the desired distribution.
Notice that what you're doing in order to sample from the second distribution is performing a perfectly uniform random experiment, but assigning numbers to outcomes so that they are distributed as you wish.
How can we sample from a hypergeometric distribution? Ideally, you would have $M$ red balls and $N-M$ blue balls handy, and could afford to spend the time blindly choosing $y$ among them and counting how many are red. Since this is rarely the case, what you can do is flip your coin a number of times, and devise a way to assign a number from $0$ to $M$ to each outcome, so that the number $x$ appears
$$
\frac{{y\choose x} {N-y \choose M-x}}{N\choose M}\times 100
$$
percent of the time. Thus, even if you never took balls out of an urn and counted the red ones, it's just as if you had. That is, if you repeat the experiment many times, the percentage of times you will get the number $x$ will get closer and closer to the number of times you would get $x$ red balls if you did the "true" experiment. That's how you "sample from a hypergeometric distribution".
In practice, one can generate a pseudorandom number uniformly distributed in $[0,1]$ using a computer, and then apply a transformation so that the result is distributed as desired. I don't know what's the precise transformation for the hypergeometric distribution, but I suppose you can find it somewhere.
A: There are already numerical routines that sample efficiently numbers distributed according to a hypergeometric distribution.
One can also generate samples of the hypergeometric distribution by  sampling from the uniform distributions in $(0,1)$. I will outline an algorithm to do so at the end.
Overview: Suppose there are $M$ white balls and $N$ black balls. We choose randomly  without replacement $k$ balls, meaning that each of the $\binom{M+N}{k}$ possible samples have the same chance of occurrence. We count the number of white balls $X$ that appear in our selection.  $X$ is a random variable whose distribution is the hypergeometric distribution with parameters M,N,k:
$$\mathbb{P}(X=j)=\frac{\binom{M}{j} \binom{N}{k-j}}{\binom{M+N}{k}},\qquad\max(0,M+N+k)\leq j\leq \min(M,k).$$
As I mentioned before, one can use  independent samples from the  uniform $(0,1)$ distribution to construct a sample of the hypergeometric distribution.
Suppose we have sample $k$ uniform independently random variables  $U_1,\ldots, U_k$. Here is an algorithm to sample a random variable from the hypergeometric distribution with parameters $M$, $N$ and $k$.
    For j=1 to k: 
     If U_j less.equal M/(N+M), 
             set B_j=w
             reset  M as M-1
     else,  
             set B_j=b
             reset N as N-1.
     continue
    
   Xhyper = sum(B==w)

where Xhyper = sum(B=w) means $X=\sum^k_{j=1}\mathbb{1}(B_j=b)$.
There are already well establish numerical recipes to produce pseudo-random numbers form the $0-1$ uniform distribution. Here is an R implementation of the algorithm above:
    ### Method 1: using samples from uniform 0-1 distribution.
    hp.geo.fun <- function(M,N,k){
      u <- runif(k)  # sample k iid 0-1 unifrom distributed numbers 
      b <- vector(mode = "character", length = k) # this records the k balls drawn
      for(j in 1:k){
        if(u[j] <= M/(M+N)){
          b[j] <- "w"
          M <- max(M-1,0)
        } else{
          b[j]  <- "b"
          N<- max(N-1,0)
        }
      }
      sum(b=="w")
    }
    ### Here we generate a sample following the algorithm described above
    M <- 10
    N <- 4
    k <- 7
    pop <- c(rep("w",M),rep("b",N))
    nsample <- 10000
    hp.geo.samp1 <- sapply(1:nsample,function(x){hp.geo.fun(M,N,k)})
    hist(hp.geo.samp1)

Efficient algorithms for sampling with/without replacement have been already implemented in many computer languages. In R for example, the  function sample that allows from sampling say k objects from a collection of n objects with/or without replacement and with certain probabilities of being chosen for each of the n objects. Also, efficient implementations of the  hypergeometric distribution are available in many languages. In R for example, there is the function rhyper based on  Kachitvichyanukul, V. and Schmeiser, B. (1985). Computer generation of hypergeometric random variates. Journal of Statistical Computation and Simulation, 22, 127–145.
#########

M <- 10
N <-  4
k <- 7
pop <- c(rep("w",M),rep("b",N))
nsample <- 10000

### Method 2: based on the sample function
hp.geo.samp.2 <- sapply(1:nsample, 
   function(x){sum(sample(pop,n,replace = FALSE)=="w")})
hist(hp.geo.samp.2)

### Method 3: based on the function rhyper
hp.geo.samp.3 <- rhyper(10000, M, N, k)
hist(hp.geo.samp.3)

I hope this clarify things for you.
