==================================================================

 Exp <- function(n, r, N = 256, k = 1, wav = "d8", ...)
{
# Ver 0.4 Brani Vidakovic, 10/9/1996.
# this function needs S+Wavelets' procedure
# dwt. All other functions used are supplied.
        ret <- list()
        eps <- 1/(2 * N)
        ret$x <- seq(eps, 1 - eps, length = N)
        #-------define lik function------------------------
        fun <- function(x, n, r)
        {
                return(exp((r - 1) * log(x) + (n - r) * log(1 - x) - lgamma(r) - 
                        lgamma(n - r + 1) + lgamma(n + 1)))
        }
#---------------------------------------------------
        ret$g <- fun(ret$x, n, r)
        ret$g <- ret$g/sum(ret$g)
        ret$h <- qnorm(ret$x)^k
        ret$gd <- as.vector(dwt(ret$g, wavelet = wav, ...))
        ret$hd <- as.vector(dwt(ret$h, wavelet = wav, ...))
        #--------------------------------------------------
        ret$res <- print(sum(ret$gd * ret$hd), digits = 18)
        #---------------------------------
        return(ret)
}


================================================================
 Exp2d <- function(n, r, s, N, eps = 0.005, wav = "d4", lev = 4, ...)
{
# Ver 0.1 Brani Vidakovic, 9/22/1996.
# this function needs S+Wavelets' procedure
# dwt. All other functions used are supplied.
        ret <- list()
        ret$x <- seq(eps, 1 - eps, length = N)
        ret$y <- ret$x  #--------------------------------------------------
        ss <<- s
        nn <<- n
        rr <<- r
        fun <- function(u, v)
        {
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                vv <- u * as.numeric(u >= v) + v * as.numeric(v > u)
                uu <- u * as.numeric(u <= v) + v * as.numeric(v < u)
                Be <- exp(lgamma(rr) + lgamma(ss - rr) + lgamma(nn - ss + 1) - 
                        lgamma(nn + 1))
                logi <- (rr - 1) * log(uu) + (ss - rr - 1) * log(vv - uu) + (nn -
                        ss) * log(1 - vv)
                return(0.5/Be * exp(logi))
        }
        ret$gm <- outer(ret$x, ret$y, fun)
        ret$gm <- ret$gm/sum(ret$gm)
        #--------------------------------------------------
        ret$h <- qnorm(ret$x)
        ret$hm <- outer(ret$h, ret$h)
        #       ret$hdm <- as.matrix(dwt.2d(ret$hm, wavelet = wav, n.levels = lev, ...)
#               )
#       ret$gdm <- as.matrix(dwt.2d(ret$gm, wavelet = wav, n.levels = lev, ...)
#               )       #--------------------------------------------------
#       ret$res <- print(sum(ret$gdm * ret$hdm), digits = 18)
#---------------------------------
        return(ret)
}

==========================================