Here, I am trying to implement Benjamini

Wikipedia provides the write up

The Original paper published in 1995

Perl has a module: Multtest on cpan, it implemented what implemented in an R module: p.adjust

Implementation in R is as following:

        i <- lp:1L
        o <- order(p, decreasing = TRUE)
        ro <- order(o)
        pmin(1, cummin(n/i * p[o]))[ro]

Given a sequence of example p-values in the original paper, R produces:

pval <- c(0.0001 ,0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, .03240, 0.4262, 0.5719, 0.6528, 0.7590, 1)
p.adjust (pval, "BH")
 [1] 0.00150000 0.00300000 0.00950000 0.03562500 0.05733333 0.05733333 0.05733333 0.05733333
 [9] 0.06885000 0.05733333 0.58118182 0.71487500 0.75323077 0.81321429 1.00000000

Therefore, FDR is control at 0.05. p(4) = 0.0356 < 0.05