numinteg_brob.r

# numerical integration methods for brob numbers

#########################################################################
# Bolstad2:::sintegral
#
###########################################################
sintegral.brob.parallel <- function(fx, sL, sH, Nsteps=100, parallel=FALSE)
{
  # rewritten from sintegral from Bolstad2 package
  sek <- seq(sL,sH,length=2*Nsteps+1)
  sek.l <- length(sek)
  h <- sek[2] - sek[1]
  if(parallel == TRUE)
  {
    require(parallel)
    if(length(grep("windows", sessionInfo())) > 0) cores <- 1 else cores <- detectCores()
    inity <- mclapply(seq_along(1:sek.l), function(x) fx(sek[x]), mc.silent = FALSE, mc.cores = cores)
  } else
  {
   inity <- lapply(seq_along(1:sek.l), function(x) fx(sek[x]))
  }
  summe <- (h/3)*( sum( list2vec.brob(inity[2 * (1:Nsteps) - 1]) ) +
                               sum( 4*list2vec.brob(inity[2 * (1:Nsteps)]) ) +
                               sum( list2vec.brob(inity[2 * (1:Nsteps) + 1]) )
                              )
return(summe)
}
#########################################################################


#########################################################################
# pracma:::romberg
#
################################################################
romberg.brob <- function (f, a, b, maxit = 25, tol = 1e-12, ...) 
{
  stopifnot(is.numeric(a), is.numeric(b), length(a) == 1, length(b) == 1)
  tol <- abs(tol)
  if (a == b) 
    return(list(value = 0, iter = 0, rel.error = 0))
  if (a > b) 
    return(-1 * romberg(f, b, a, tol = tol, ...))
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  eps <- .Machine$double.eps
  
  #infinity check boundaries
  #if (!is.finite(f(a))) 
  #  a <- a + eps * sign(b - a)
  #if (!is.finite(f(b))) 
  #b <- b - eps * sign(b - a)
  
  #I <- matrix(0, nrow = maxit + 1, ncol = maxit + 1)
  # matrix as brob list
  # matrix[row, col] = brob.list[[row]][col]
  I.brob <- vector(mode="list", length=maxit + 1)
  I.brob <- lapply(I.brob, function(x) x <- as.brob(rep(0, maxit + 1)))

  n <- 1
  iter <- 0
  err <- 1
  
  while (err > tol && iter < maxit) {
    iter <- iter + 1
    n <- 2 * n
    h <- (b - a)/n
    S <- f(a)
    for (i in 1:(n - 1))
    {
      xi <- a + h * i
      S <- S + 2 * f(xi)
    }
    S <- (S + f(b)) * h/2

    #I[iter + 1, 1] <- S
    I.brob[[iter +1]][1] <- S
    
    for (k in 2:(iter + 1)) {
      j <- 2 + iter - k
      #I[j, k] <- (4^(k - 1) * I[j + 1, k - 1] - I[j, k - 1])/(4^(k - 1) - 1)
      I.brob[[j]][k] <- ( 4^(k - 1) * I.brob[[j + 1]][k - 1] - I.brob[[j]][k - 1] )/ ( 4^(k - 1) - 1 )
    }
    err <- abs(I.brob[[1]][iter + 1] - I.brob[[2]][iter] )
  }
  
  #while (err > tol && iter < maxit) {
  #  iter <- iter + 1
  #  n <- 2 * n
  #  h <- (b - a)/n
  #  S <- f(a)
  #  for (i in 1:(n - 1)) {
  #    xi <- a + h * i
  #    S <- S + 2 * f(xi)
  #  }
  #  S <- (S + f(b)) * h/2
  #  I[iter + 1, 1] <- S
  #  for (k in 2:(iter + 1)) {
  #    j <- 2 + iter - k
  #    I[j, k] <- (4^(k - 1) * I[j + 1, k - 1] - I[j, k - 1])/(4^(k - 1) - 1)
  #  }
  #  err <- abs(I[1, iter + 1] - I[2, iter])
  #}
  
  if (iter == maxit) 
    warning("Maximum number of iterations has been reached.")
return(list(value = I.brob[[1]][iter + 1], iter = iter, rel.error = err))
}
#########################################################################


#########################################################################
# pracma:::cotes
#
# tweaked 'cotes' from 'pracma' to work with brob
# TODO parallel computing version
#
########################################################
cotes.brob <- function (f, a, b, n = 20, nodes = 5, ...) 
{
  if (nodes < 2 || nodes > 8) 
    stop("Number of nodes, 'nodes', must be between 2 and 8.")
  if (n < nodes) 
    stop("Argument 'n' must be greater or equal to the number of nodes.")
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  N <- (nodes - 1) * ceiling(n/(nodes - 1))
  N1 <- N + 1
  x <- linspace(a, b, N1)
  h <- x[2] - x[1]
  
  #g <- f(x)
  g <- list2vec.brob( lapply(x,function(x) f(x)) )
  
  if (length(g) != N1) {
    #g <- numeric(N1)
    g <- as.brob(rep(0,N1))
    for (i in 1:N1) g[i] <- f(x[i])
  }
  endpts <- g[1] + g[N1]
  switch(nodes - 1,
         
         Q <- (h/2) * (endpts + 2 * sum(g[2:N]))
         , 
         
         Q <- (h/3) * (endpts +
                         4 * sum(g[seq(2, N, 2)]) +
                         2 * sum(g[seq(3, N, 2)]))
         ,
         
         Q <- (3 * h/8) * (endpts +
                             3 * sum(g[seq(2, N, 3)] +
                                       g[seq(3, N, 3)]) +
                             2 * sum(g[seq(4, N, 3)]))
         ,
         
         Q <- (4 * h/90) * (7 * endpts +
                              32 * sum(g[seq(2, N, 4)]) +
                              12 * sum(g[seq(3, N, 4)]) +
                              32 * sum(g[seq(4, N, 4)]) +
                              14 * sum(g[seq(5, N, 4)]))
         , 
         
         Q <- (5 * h / 288) * (19 * endpts +
                                 75 * sum(g[seq(2, N, 5)] +
                                            g[seq(5, N, 5)]) +
                                 50 * sum(g[seq(3, N, 5)] +
                                            g[seq(4, N, 5)]) +
                                 38 * sum(g[seq(6, N,5)]))
         ,
         
         Q <- (6 * h / 840) * (41 * endpts +
                                 216 * sum(g[seq(2,N, 6)] +
                                             g[seq(6, N, 6)]) +
                                 27 * sum(g[seq(3, N,6)] +
                                            g[seq(5, N, 6)]) +
                                 272 * sum(g[seq(4, N, 6)]) +
                                 82 * sum(g[seq(7, N, 6)]))
         ,
         
         Q <- (7 * h / 17280) * (751 * endpts +
                                   3577 * sum(g[seq(2, N, 7)] +
                                                g[seq(7,N, 7)]) +
                                   1323 * sum(g[seq(3, N, 7)] +
                                                g[seq(6,N, 7)]) +
                                   2989 * sum(g[seq(4, N, 7)] +
                                                g[seq(5,N, 7)]) +
                                   1502 * sum(g[seq(8, N, 7)]))
  )
  
return(Q)
}
#########################################################################


######################################################################### 
# pracma:::integral
#
# TODO: add "Clenshaw"
#############################################################################
integral.brob <- function (fun, xmin, xmax, method = c("Kronrod", #"Clenshaw", 
                           "Simpson"), no_intervals = 8, random = FALSE, reltol = 1e-08, 
                           abstol = 0, ...) 
{
  stopifnot(is.numeric(xmin), length(xmin) == 1, is.numeric(xmax), 
            length(xmax) == 1)
  no_intervals <- max(1, floor(no_intervals))
  fun <- match.fun(fun)
  f <- function(x) fun(x, ...)
  if (length(f(xmin)) > 1 || length(f(xmax)) > 1) {
    stop("Function 'fun' is array-valued! Use 'quadv'.\n")
  }
  if (length(f(c(xmin, xmax))) != 2) {
    cat("Warning: Function 'fun' is not vectorized!\n")
    f = Vectorize(f)
  }
  if (xmin == xmax) 
    return(0)
  method <- match.arg(method)
  tol <- if (abstol > 0) 
    min(reltol, abstol)
  else reltol
  if (is.infinite(xmin) || is.infinite(xmax)) {
    cat("For infinite domains Gauss integration is applied!\n")
    Q <- quadinf(f, xmin, xmax, tol = tol)$Q
    return(Q)
  }
  if (random) {
    xs <- c(xmin, (xmax - xmin) * sort(runif(no_intervals - 1)) + xmin, xmax)
  } else
    {
    xs <- linspace(xmin, xmax, no_intervals + 1)
    }
  Q <- 0
  if (method == "Kronrod") {
    for (i in 1:no_intervals) {
      Q = Q + quadgk.brob(f, xs[i], xs[i + 1], tol = tol)
    }
  }
  else if (method == "Clenshaw") {
    stop("Method 'Clenshaw' not implemented for very large numbers.")
  #  for (i in 1:no_intervals) {
  #    Q = Q + quadcc(f, xs[i], xs[i + 1], tol = tol)
  #  }
  }
  else if (method == "Simpson") {
    for (i in 1:no_intervals) {
      Q = Q + simpadpt.brob(f, xs[i], xs[i + 1], tol = tol)
    }
  }
  else {
    stop("Unknown method; not available as integration routine.")
  }
  return(Q)
}
#########################################################################


#########################################################################
# pracma:::quadgk
#
#####################################################################
quadgk.brob <- function (f, a, b, tol = .Machine$double.eps^0.5, ...) 
{
  stopifnot(is.numeric(a), length(a) == 1, is.numeric(b), length(b) == 
              1)
  eps <- .Machine$double.eps
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  if (a == b) 
    return(0)
  else if (a > b) 
    return(-1 * quadgk(f, b, a, tol = tol))
  n15 <- c(-0.991455371120813, -0.949107912342758, -0.864864423359769, 
           -0.741531185599394, -0.586087235467691, -0.405845151377397, 
           -0.207784955007899, 0, 0.207784955007899, 0.405845151377397, 
           0.586087235467691, 0.741531185599394, 0.864864423359769, 
           0.949107912342758, 0.991455371120813)
  n7 <- c(-0.949107912342758, -0.741531185599394, -0.405845151377397, 
          0, 0.405845151377397, 0.741531185599394, 0.949107912342758)
  w15 <- c(0.0229353220105292, 0.0630920926299785, 0.10479001032225, 
           0.140653259715526, 0.169004726639268, 0.190350578064785, 
           0.204432940075299, 0.209482141084728, 0.204432940075299, 
           0.190350578064785, 0.169004726639268, 0.140653259715526, 
           0.10479001032225, 0.0630920926299785, 0.0229353220105292)
  w7 <- c(0.12948496616887, 0.279705391489277, 0.381830050505119, 
          0.417959183673469, 0.381830050505119, 0.279705391489277, 
          0.12948496616887)
  
  # pracma:::.gkadpt
  .gkadpt <- function(f, a, b, tol = tol) {
    x15 <- 0.5 * ((b - a) * n15 + b + a)
    x7 <- 0.5 * ((b - a) * n7 + b + a)
    Q7 <- sum(w7 * f(x7)) * (b - a)/2
    Q15 <- sum(w15 * f(x15)) * (b - a)/2
    
    #infinity check boundaries
    #if (!is.finite(Q7) || !is.finite(Q15)) {
    #  warning("Infinite or NA function value encountered.")
    #  return(Q15)
    #} else
    if (abs(Q15 - Q7) < tol) {
      return(Q15)
    }
    else if (abs(b - a) < 16 * eps) {
      warning("Minimum step size reached; singularity possible.")
      return(Q2)
    }
    Q2 <- .gkadpt(f, (a + b)/2, b, tol = tol)
    Q1 <- .gkadpt(f, a, (a + b)/2, tol = tol)
    return(Q1 + Q2)
  }
  .gkadpt(f, a, b, tol = tol)
}
#########################################################################


#########################################################################
# pracma:::simpadpt
#
#####################################################
simpadpt.brob <- function (f, a, b, tol = 1e-06, ...) 
{
  stopifnot(is.numeric(a), length(a) == 1, is.numeric(b), length(b) == 
              1)
  eps <- .Machine$double.eps
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  if (a == b) 
    return(0)
  else if (a > b) 
    return(-1 * simpadpt.brob(f, b, a, tol = tol))
  h <- 1/8 * (b - a)
  x <- c(a, a + h, a + 2 * h, (a + b)/2, b - 2 * h, b - h, b)
  y <- c(f(a), f(a + h), f(a + 2 * h), f((a + b)/2), f(b - 2 * h), f(b - h), f(b))
  
  #infinity check boundaries
  #if (!is.finite(y[1])) 
  #  y[1] <- f(a + eps * (b - a))
  #if (!is.finite(y[7])) 
  #  y[7] <- f(b - eps * (b - a))
  
  hmin <- eps * (b - a)/1024
  Q1 <- .simpadpt.brob(f, x[1], x[3], y[1], y[2], y[3], tol, hmin)
  Q2 <- .simpadpt.brob(f, x[3], x[5], y[3], y[4], y[5], tol, hmin)
  Q3 <- .simpadpt.brob(f, x[5], x[7], y[5], y[6], y[7], tol, hmin)
  return(Q1 + Q2 + Q3)
}

# pracma:::.simpadpt
.simpadpt.brob <- function (f, a, b, fa, fc, fb, tol, hmin) 
{
  h <- b - a
  g <- (a + b)/2
  d <- (a + g)/2
  fd <- f(d)
  e <- (g + b)/2
  fe <- f(e)

  if(is.list(fa)) fa <- list2vec.brob(fa)
  if(is.list(fc)) fc <- list2vec.brob(fc)
  if(is.list(fb)) fb <- list2vec.brob(fb)
  
  Q1 <- (h/6) * (fa + 4 * fc + fb)
  Q2 <- (h/12) * (fa + 4 * fd + 2 * fc + 4 * fe + fb)
  Q <- Q2 + (Q2 - Q1)/15
  
  #infinity check boundaries
  #if (!is.finite(Q)) {
  #  warning("Infinite or NA function value encountered.")
  #  return(Q)
  #} else
  if (abs(Q2 - Q) <= tol) {
    return(Q)
  } else if (abs(h) < hmin || g == a || g == b) {
    warning("Minimum step size reached; singularity possible.")
    return(Q)
  }
  Q4 <- .simpadpt.brob(f, a, g, fa, fd, fc, tol, hmin)
  Q5 <- .simpadpt.brob(f, g, b, fc, fe, fb, tol, hmin)
  return(Q4 + Q5)
}
#########################################################################

#########################################################################
# pracma:::quadv
#
# seems to be just a vectorized version, no need to add this
# if one wants a vectorized version that works on multiple functions
# in parallel, one can write a short wrapper around an integral
# function
#
# otherwise it would require to simulate matrices for brob objects
#
#########################################################################


#########################################################################
# pracma:::quadinf
#
# seems to work without modification for brob
#
######################################################
quadinf.brob <- function (f, xa, xb, tol = 1e-12, ...) 
{
  stopifnot(is.numeric(xa), length(xa) == 1, is.numeric(xb), 
            length(xb) == 1, is.numeric(tol), length(tol) == 1)
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  if (xa == xb) {
    return(list(Q = 0, reltol = 0))
  } else
  if (xa > xb) {
    stop("For the integration limits 'xa <= xb' is required.")
  }
  u <- -1
  v <- 1
  if (is.finite(xa) && is.finite(xb)) {
    duv <- (xb - xa)/(v - u)
    ff <- function(y) duv * f(xa + duv * (y - u))
  } else
  if (is.finite(xa) && is.infinite(xb)) {
    ff <- function(y) {
        duv <- (v - u)/(1 - y)^2
        z <- (y - u)/(v - u)
        duv * f(xa + z/(1 - z))
        }
  } else
  if (is.infinite(xa) && is.finite(xb)) {
    xa <- -xb
    ff <- function(y) {
      duv = (v - u)/(1 - y)^2
      z <- (y - u)/(v - u)
      duv * f(-(xa + z/(1 - z)))
    }
  } else
  if (is.infinite(xa) && is.infinite(xb)) {
    ff <- function(y) {
      z <- pi * (y - u)/(v - u) - 0.5
      duv <- (1/cos(z)^2) * pi/(v - u)
      duv * f(tan(z))
    }
  } else
  {
    stop("Other integration domains will be treated later.")
  }
  cc <- .quadinf_pre.brob()
  xx <- cc$nodes
  ww <- cc$weights
  h <- 1/2
  x <- xx[[1]]
  w <- ww[[1]]
  s <- w[7] * ff(x[7])
  for (j in 1:6) {
    s <- s + w[j] * (ff(x[j]) + ff(-x[j]))
  }
  Q <- s * h
  for (k in 2:7) {
    x <- xx[[k]]
    w <- ww[[k]]
    s <- 0
    for (j in 1:length(w)) {
      s <- s + w[j] * (ff(x[j]) + ff(-x[j]))
    }
    h <- h/2
    newQ <- s * h + Q/2
    delta <- abs(newQ - Q)
    Q <- newQ
    if (delta < tol) 
      break
  }
  return(list(Q = Q, relerr = delta, niter = k))
}

# pracma:::.quadinf_pre
# can be used without any modification, use it to be independet from 'pracma'
.quadinf_pre.brob <- function () 
{
  nodes <- list(c(0.999999999999957, 0.999999988875665, 0.999977477192462, 
                  0.997514856457224, 0.951367964072747, 0.674271492248436, 
                  0),
                c(0.999999999952856, 0.999999204737115, 0.999688264028353, 
                  0.987040560507377, 0.859569058689897, 0.377209738164034), 
                c(0.999999999998232, 0.999999999142705, 0.999999892781612, 
                  0.99999531604122, 0.999909384695144, 0.999065196455786, 
                  0.994055506631402, 0.973966868195677, 0.914879263264575, 
                  0.7806074389832, 0.539146705387968, 0.194357003324935), 
                c(0.999999999999708, 0.999999999990394, 0.999999999789733, 
                  0.999999996787199, 0.999999964239081, 0.999999698894153, 
                  0.999998017140595, 0.999989482014818, 0.999953871005628, 
                  0.999828822072875, 0.999451434435275, 0.998454208767698, 
                  0.996108665437509, 0.991126992441699, 0.981454826677335, 
                  0.964112164223547, 0.935160857521985, 0.88989140278426, 
                  0.823317005506402, 0.731018034792561, 0.610273657500639, 
                  0.461253543939586, 0.287879932742716, 0.0979238852878323), 
                c(0.999999999999886, 0.999999999999271, 0.999999999995825, 
                  0.999999999978455, 0.999999999899278, 0.999999999570788, 
                  0.999999998323362, 0.999999993964134, 0.999999979874503, 
                  0.999999937554078, 0.999999818893713, 0.999999507005719, 
                  0.999998735471866, 0.99999693244919, 0.999992937876663, 
                  0.999984519902271, 0.999967593067943, 0.999935019925082, 
                  0.99987486504878, 0.999767971599561, 0.999584750351518, 
                  0.999281111921792, 0.998793534298806, 0.998033336315434, 
                  0.996880318128192, 0.995176026155327, 0.992716997196827, 
                  0.989248431090134, 0.984458831167431, 0.977976235186665, 
                  0.969366732896917, 0.958136022710214, 0.943734786052757, 
                  0.925568634068613, 0.903013281513574, 0.875435397630409, 
                  0.842219246350757, 0.802798741343241, 0.75669390863373, 
                  0.703550005147142, 0.643176758985205, 0.575584490635152, 
                  0.501013389379309, 0.419952111278447, 0.333142264577638, 
                  0.241566319538884, 0.146417984290588, 0.0490559673050779), 
                c(0.99999999999993, 0.999999999999817, 0.999999999999537, 
                  0.999999999998861, 0.999999999997274, 0.999999999993646, 
                  0.999999999985569, 0.999999999968033, 0.999999999930888, 
                  0.999999999854056, 0.999999999698754, 0.999999999391777, 
                  0.999999998797983, 0.999999997673237, 0.999999995585634, 
                  0.999999991786456, 0.999999985003076, 0.999999973113236, 
                  0.999999952642664, 0.999999918004795, 0.999999860371215, 
                  0.999999766023332, 0.99999961398855, 0.999999372707335, 
                  0.99999899541069, 0.999998413810965, 0.999997529623805, 
                  0.999996203347166, 0.999994239627617, 0.999991368448345, 
                  0.999987221282001, 0.999981301270121, 0.999972946425232, 
                  0.999961284807857, 0.999945180614459, 0.999923170129289, 
                  0.999893386547593, 0.999853472773111, 0.999800481431138, 
                  0.999730761519808, 0.9996398313456, 0.999522237651217, 
                  0.999371401140938, 0.999179448934886, 0.998937034833512, 
                  0.998633148640677, 0.998254916171996, 0.997787391958906, 
                  0.997213347043469, 0.996513054640254, 0.995664076816953, 
                  0.994641055712511, 0.993415513169264, 0.991955663002678, 
                  0.990226240467528, 0.988188353800743, 0.985799363025283, 
                  0.983012791481101, 0.979778275800616, 0.976041560256577, 
                  0.971744541565487, 0.966825370312356, 0.961218615151116, 
                  0.954855495805023, 0.947664190615153, 0.939570223933275, 
                  0.930496937997153, 0.920366053031953, 0.90909831816302, 
                  0.896614254280076, 0.882834988244669, 0.867683175775646, 
                  0.851084007987849, 0.832966293919411, 0.813263608502974, 
                  0.791915492376142, 0.768868686768246, 0.744078383547347, 
                  0.717509467487324, 0.689137725061668, 0.65895099174335, 
                  0.626950208051043, 0.593150353591953, 0.557581228260778, 
                  0.52028805069123, 0.481331846116905, 0.440789599033901, 
                  0.398754150467238, 0.355333825165075, 0.310651780552846, 
                  0.264845076583448, 0.218063473469712, 0.170467972382011, 
                  0.122229122201558, 0.0735251229856713, 0.0245397635746492), 
                c(0.999999999999945, 0.999999999999911, 0.999999999999856, 
                  0.999999999999769, 0.999999999999632, 0.999999999999419, 
                  0.999999999999088, 0.99999999999858, 0.999999999997803, 
                  0.999999999996623, 0.999999999994845, 0.999999999992181, 
                  0.999999999988217, 0.999999999982353, 0.999999999973737, 
                  0.99999999996115, 0.999999999942877, 0.999999999916506, 
                  0.99999999987867, 0.999999999824698, 0.999999999748146, 
                  0.999999999640173, 0.999999999488718, 0.999999999277422, 
                  0.99999999898421, 0.999999998579458, 0.999999998023619, 
                  0.999999997264174, 0.999999996231727, 0.999999994835051, 
                  0.999999992954804, 0.999999990435633, 0.999999987076265, 
                  0.999999982617173, 0.999999976725267, 0.99999996897499, 
                  0.999999958825101, 0.999999945590273, 0.999999928406514, 
                  0.999999906189268, 0.999999877582855, 0.999999840899736, 
                  0.999999794047876, 0.999999734444233, 0.999999658912159, 
                  0.999999563560232, 0.999999443639729, 0.999999293377668, 
                  0.999999105781996, 0.999998872415162, 0.999998583131984, 
                  0.999998225777311, 0.999997785838635, 0.999997246048426, 
                  0.999996585930571, 0.999995781284928, 0.999994803603649, 
                  0.999993619412539, 0.999992189530434, 0.999990468239214, 
                  0.999988402356833, 0.999985930205475, 0.999982980466756, 
                  0.999979470915751, 0.999975307025492, 0.999970380433589, 
                  0.999964567262626, 0.999957726286078, 0.999949696931689, 
                  0.999940297114479, 0.999929320891884, 0.999916535933951, 
                  0.999901680802003, 0.999884462029766, 0.999864551001635, 
                  0.999841580623482, 0.999815141782273, 0.999784779591651, 
                  0.999749989421649, 0.999710212711751, 0.999664832567665, 
                  0.999613169143347, 0.999554474811097, 0.999487929123823, 
                  0.999412633574952, 0.99932760616283, 0.999231775767899, 
                  0.99912397635239, 0.999002940993729, 0.998867295764364, 
                  0.998715553472209, 0.998546107277413, 0.998357224202664, 
                  0.998147038555748, 0.997913545284577, 0.997654593286378, 
                  0.997367878694235, 0.997050938165609, 0.996701142198913, 
                  0.996315688505661, 0.9958915954671, 0.995425695705623, 
                  0.994914629802639, 0.994354840195921, 0.993742565290762, 
                  0.993073833820584, 0.992344459493918, 0.991550035965892, 
                  0.990685932173615, 0.989747288075987, 0.988729010839602, 
                  0.987625771513511, 0.986432002236608, 0.985141894022398, 
                  0.983749395166731, 0.982248210324941, 0.980631800305449, 
                  0.978893382627494, 0.977025932891058, 0.975022187007306, 
                  0.972874644337964, 0.9705755717919, 0.968117008926856, 
                  0.96549077410362, 0.962688471739078, 0.959701500703337, 
                  0.956521063904581, 0.953138179103394, 0.949543690995936, 
                  0.945728284602632, 0.941682499995769, 0.937396748395719, 
                  0.93286132966124, 0.928066451194555, 0.923002248276554, 
                  0.917658805841549, 0.912026181694487, 0.906094431166404, 
                  0.89985363319617, 0.893293917818211, 0.886405495026979, 
                  0.879178684979389, 0.871603949486378, 0.863671924734105, 
                  0.855373455164271, 0.846699628431464, 0.8376418113436, 
                  0.828191686679357, 0.818341290764098, 0.808083051673338, 
                  0.797409827920312, 0.78631494747182, 0.774792246924435, 
                  0.762836110661392, 0.75044150979924, 0.737604040722825, 
                  0.724319962997407, 0.710586236438006, 0.696400557108459, 
                  0.681761392016427, 0.666668012265738, 0.651120524424304, 
                  0.635119899864422, 0.618668001832728, 0.601767610009633, 
                  0.584422442322664, 0.566637173785019, 0.548417452139792, 
                  0.529769910101775, 0.510702174002558, 0.491222868660811, 
                  0.471341618317998, 0.451069043500452, 0.430416753691437, 
                  0.409397335721529, 0.388024337812118, 0.366312249234904, 
                  0.344276475579705, 0.321933309653369, 0.29929989806396, 
                  0.276394203576179, 0.253234963356, 0.229841643254361, 
                  0.206234388311029, 0.182433969690289, 0.158461728289299, 
                  0.134339515287672, 0.110089629932628, 0.085734754877651, 
                  0.06129788941366, 0.0368022809500251, 0.0122713551180822))
  weights = list(c(1.35817842745391e-12, 2.1431204556943e-07, 
                  0.000266200513752717, 0.0183431669899278, 0.230022394514789, 
                  0.965976579412301, 1.5707963267949), c(1.16311658142558e-09, 
                  1.19837013631707e-05, 0.00290251774790131, 0.0763857435708323, 
                  0.531078275428054, 1.38961475924726), c(4.93785387766319e-11, 
                  1.86872822687364e-08, 1.82633205937107e-06, 6.24825592407441e-05, 
                  0.000949946804283469, 0.00774260102606424, 0.0391750054936008, 
                  0.137422107733168, 0.360461418469344, 0.737437848361548, 
                  1.19346302584916, 1.52328371863471), c(8.6759314149796e-12, 
                  2.52163479185301e-10, 4.87600609742406e-09, 6.58351851271834e-08, 
                  6.47775660359297e-07, 4.82371820326155e-06, 2.81101643279401e-05, 
                  0.0001320523412561, 0.000513393824067903, 0.00169087399814264, 
                  0.00481629814392846, 0.012083543599158, 0.027133510013712, 
                  0.0552896837422406, 0.103432154223333, 0.179324412110728, 
                  0.290240679312454, 0.440833236273858, 0.630405135164744, 
                  0.85017285645662, 1.08163498549007, 1.29747575042498, 
                  1.46601442671697, 1.55877335553333), c(3.48419376702611e-12, 
                  2.09893354045115e-11, 1.13060553474947e-10, 5.4828357797095e-10, 
                  2.40917732564759e-09, 9.64988889610896e-09, 3.54347771714219e-08, 
                  1.19924427829028e-07, 3.75954118623606e-07, 1.09688351259013e-06, 
                  2.99166158781388e-06, 7.65957585252031e-06, 1.84818135998792e-05, 
                  4.21831838417576e-05, 9.13908174907101e-05, 0.000188564429767003, 
                  0.000371666936216778, 0.000701859515684242, 0.00127332794470824, 
                  0.00222508270647864, 0.00375425097743183, 0.00613003763208303, 
                  0.00970722373939169, 0.0149378350960501, 0.0223794710636485, 
                  0.032698732726609, 0.0466682080548466, 0.0651555334325362, 
                  0.0891031392409415, 0.119497411288696, 0.157326203484366, 
                  0.203523998858602, 0.258904639514054, 0.324082539611529, 
                  0.399384741525717, 0.484758091214755, 0.579678103087788, 
                  0.683068516344264, 0.793242700820517, 0.907879379154895, 
                  1.02404493311181, 1.13827224337631, 1.24670120745186, 
                  1.34527888476625, 1.4300083548723, 1.49722622254104, 
                  1.54388111617696, 1.56778143130722), c(2.18359220992336e-12, 
                  5.51823694681749e-12, 1.35425129123363e-11, 3.23044643332524e-11, 
                  7.49673975738182e-11, 1.69394577894116e-10, 3.72995018430528e-10, 
                  8.00997844797297e-10, 1.67888976821619e-09, 3.43718567446501e-09, 
                  6.8784610955899e-09, 1.3464645522302e-08, 2.57995682295359e-08, 
                  4.84209501980724e-08, 8.90713951402424e-08, 1.60693945790762e-07, 
                  2.84499236591598e-07, 4.94582887027542e-07, 8.44737563848599e-07, 
                  1.41830671554939e-06, 2.34216672085281e-06, 3.80619832646449e-06, 
                  6.0899100320949e-06, 9.59819412837847e-06, 1.49085140318706e-05, 
                  2.28321181090361e-05, 3.44921247593432e-05, 5.14214974476588e-05, 
                  7.56839965862015e-05, 0.000110021128466667, 0.000158027884007012, 
                  0.000224359652050085, 0.000314972091860212, 0.000437394956159117, 
                  0.000601039879911474, 0.000817541013324695, 0.00110112611345194, 
                  0.00146901435994298, 0.00194183577598437, 0.00254406576752917, 
                  0.00330446699403483, 0.00425652959901786, 0.005438899797624, 
                  0.006895785969066, 0.00867733074953918, 0.0108399371682559, 
                  0.0134465366052857, 0.0165667862542476, 0.0202771838175001, 
                  0.0246610873147533, 0.0298086281173101, 0.0358165056041964, 
                  0.0427876521577257, 0.0508307575725705, 0.0600596423586363, 
                  0.070592469906867, 0.0825507881107017, 0.0960583918651895, 
                  0.111239998988745, 0.128219733631201, 0.147119413257857, 
                  0.168056637948269, 0.191142684133427, 0.216480209117296, 
                  0.24416077786984, 0.274262229689068, 0.306845909417917, 
                  0.341953795923017, 0.379605569386652, 0.419795668445015, 
                  0.462490398055368, 0.507625158831908, 0.555101878003633, 
                  0.604786730578404, 0.656508246131628, 0.710055901205469, 
                  0.765179298908956, 0.821588035266965, 0.878952345552782, 
                  0.936904612745668, 0.995041804046133, 1.05292887995527, 
                  1.11010319396534, 1.16607986993243, 1.22035810957936, 
                  1.27242834553786, 1.32178011744377, 1.3679105116809, 
                  1.41033297144626, 1.44858625496132, 1.48224329788554, 
                  1.51091972307417, 1.5342817381543, 1.55205316984541, 
                  1.56402140377323, 1.57004202927959), c(1.72376440360427e-12, 
                  2.76085876713983e-12, 4.38886518997793e-12, 6.92552801526814e-12, 
                  1.08491618343371e-11, 1.6874519343915e-11, 2.60619605028053e-11, 
                  3.99733895199303e-11, 6.08933870643807e-11, 9.21405642265189e-11, 
                  1.38502875258341e-10, 2.06842035390292e-10, 3.06927020787233e-10, 
                  4.52575761041538e-10, 6.63208634701621e-10, 9.65948518580993e-10, 
                  1.39844143125654e-09, 2.01262102188672e-09, 2.87970134644711e-09, 
                  4.09675791396852e-09, 5.79534957309661e-09, 8.15274648285765e-09, 
                  1.14064655548505e-08, 1.58729780909682e-08, 2.19716489170893e-08, 
                  3.02551964648995e-08, 4.14482335585054e-08, 5.64957638716706e-08, 
                  7.66238739748176e-08, 1.03415280407456e-07, 1.3890286996036e-07, 
                  1.85684913700251e-07, 2.47066245112497e-07, 3.27230373305176e-07, 
                  4.31448255871654e-07, 5.66330284020508e-07, 7.40128934909086e-07, 
                  9.63100521176592e-07, 1.24793551211775e-06, 1.61026800945077e-06, 
                  2.06927612573647e-06, 2.64838622540434e-06, 3.37609523480551e-06, 
                  4.28692649399982e-06, 5.42253589182403e-06, 6.83298627742185e-06, 
                  8.57820935370791e-06, 1.07296754068523e-05, 1.33722922845324e-05, 
                  1.66065559765083e-05, 2.05509759449342e-05, 2.53447989688501e-05, 
                  3.11510556774483e-05, 3.81599541202458e-05, 4.65926446304946e-05, 
                  5.67053798539329e-05, 6.87940931134964e-05, 8.31994172400365e-05, 
                  0.000100312164601124, 0.00012057928729057, 0.000144510334290946, 
                  0.000172684419885929, 0.000205757714679988, 0.000244471467286892, 
                  0.000289660561088772, 0.000342262606462977, 0.000403327564549541, 
                  0.000474027894018167, 0.000555669207425785, 0.00064970141867429, 
                  0.000757730357827359, 0.00088152982417297, 0.00102305404297454, 
                  0.00118445048589028, 0.0013680730096097, 0.00157649526191056, 
                  0.00181252429912886, 0.00207921435400867, 0.00237988068810044, 
                  0.00271811345835023, 0.00309779152330082, 0.00352309611044299, 
                  0.00399852426273353, 0.00452890197915629, 0.00511939696145378, 
                  0.00577553087680655, 0.00650319104428258, 0.00730864145131325, 
                  0.0081985330052612, 0.0091799129243109, 0.0102602331714108, 
                  0.0114473578347998, 0.0127495693587312, 0.0141755735283305, 
                  0.0157345031130598, 0.0174359200739775, 0.0192898162408456, 
                  0.0213066123661261, 0.0234971554639885, 0.0258727143436176, 
                  0.0284449732473342, 0.0312260235053317, 0.0342283531201757, 
                  0.037464834195629, 0.0409487081258673, 0.0446935684627655, 
                  0.0487133413807014, 0.0530222636602872, 0.0576348581146547, 
                  0.0625659063844551, 0.067830419030658, 0.0734436028576388, 
                  0.0794208254030082, 0.0857775765352682, 0.0925294271057785, 
                  0.0996919846077886, 0.107280845802556, 0.115311546280937, 
                  0.123799506938413, 0.132759977352446, 0.142207976063402, 
                  0.1521582277742, 0.162625097499384, 0.173622521711631, 
                  0.185163936552775, 0.197262203197437, 0.209929530480208, 
                  0.223177394922185, 0.237016458319419, 0.25145648308451, 
                  0.266506245563135, 0.282173447579596, 0.298464626499445, 
                  0.315385064132678, 0.332938694837748, 0.351128013224425, 
                  0.369953981892114, 0.389415939679212, 0.409511510938194, 
                  0.430236516389789, 0.451584886147545, 0.473548575540614, 
                  0.496117484397316, 0.519279380484246, 0.543019827824843, 
                  0.567322120646739, 0.592167223728207, 0.617533719929909, 
                  0.643397765708253, 0.669733055410291, 0.696510795146547, 
                  0.72369968702683, 0.751265924524357, 0.779173199704798, 
                  0.807382723018763, 0.835853256308267, 0.864541159619669, 
                  0.893400452347196, 0.922382889152455, 0.951438051016353, 
                  0.980513451680855, 1.00955465962943, 1.03850543563739, 
                  1.0673078857975, 1.09590262979297, 1.1242289840506, 1.15222515926255, 
                  1.17982847161733, 1.20697556693131, 1.23360265672195, 
                  1.25964576511667, 1.28504098534673, 1.30972474443744, 
                  1.33363407457576, 1.35670688951561, 1.37888226427314, 
                  1.40010071626944, 1.42030448599969, 1.43943781524641, 
                  1.45744722081255, 1.47428176172808, 1.4898932978833, 
                  1.50423673806368, 1.51727027540505, 1.52895560835458, 
                  1.53925814531188, 1.54814719123556, 1.55559611463166, 
                  1.56158249349181, 1.56608823891746, 1.56909969535167, 
                  1.57060771653828))
  return(list(nodes = nodes, weights = weights))
}
#########################################################################


#########################################################################
# pracma:::quad
#
################################################################
quad.brob <- function (f, xa, xb, tol = .Machine$double.eps^0.5,
                       TRACE = FALSE, digs=6, ...) 
{
  stopifnot(is.numeric(xa), length(xa) == 1, is.finite(xa), 
            is.numeric(xb), length(xb) == 1, is.finite(xb))
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  if (xa == xb) 
  {
    return(xb - xa)
  } else
  if (xa > xb) {
    tmp <- xa
    xa <- xb
    xb <- tmp
    rev_p <- TRUE
  } else
  {
    rev_p <- FALSE
  }  
  eps <- .Machine$double.eps

  #infinity check boundaries
  #if (!is.finite(f(xa))) 
  #  xa <- xa + 2 * eps
  #if (!is.finite(f(xb))) 
  #  xb <- xb - 2 * eps
  
  Q <- .adaptsim.brob(f, xa, xb, tol, TRACE, digs)
  if (rev_p) Q <- -1 * Q
return(Q)
}

#pracma:::.adaptsim
.adaptsim.brob  <- function (f, xa, xb, tol = tol, TRACE = TRACE, digs=6) 
{
  x <- c(xa, (xa + xb)/2, xb)
  y <- c(f(xa), f((xa + xb)/2), f(xb))
  fa <- y[1]
  fm <- y[2]
  fb <- y[3]
  yy <- f(xa + c(0.9501, 0.2311, 0.6068, 0.486, 0.8913) * (xb - xa))
  ab <- (xb - xa)/8 * ( sum(list2vec.brob(y)) + sum(yy))
  if (ab == 0) ab <- xb - xa
  ab <- ab * tol/.Machine$double.eps
  Q <- .adaptsimstp.brob(f, xa, xb, fa, fm, fb, ab, TRACE, digs)
return(Q)
}

#pracma:::.adaptsimstp
.adaptsimstp.brob <- function (f, xa, xb, fa, fm, fb, ab, TRACE, digs=6) 
{
  m <- (xa + xb)/2
  h <- (xb - xa)/4
  x <- c(xa + h, xb - h)
  y <- c(f(xa + h), f(xb - h))
  fml <- y[1]
  fmr <- y[2]
  
  #brob
  if(is.list(fa)) fa <- list2vec.brob(fa)
  if(is.list(fb)) fb <- list2vec.brob(fb)
  if(is.list(fm)) fm <- list2vec.brob(fm)
  if(is.list(fml)) fml <- list2vec.brob(fml)
  if(is.list(fmr)) fmr <- list2vec.brob(fmr)
  
  i1 <- h/1.5 * (fa + 4 * fm + fb)
  i2 <- h/3 * (fa + 4 * (fml + fmr) + 2 * fm + fb)
  i1 <- (16 * i2 - i1)/15
  if ((ab + (i1 - i2) == ab) || (m <= xa) || (xb <= m)) {
    if (((m <= xa) || (xb <= m))) 
      warning("Required tolerance may not be met.")
    Q <- i1
    if (TRACE) 
      cat("xa = ",signif(xa,digs),"\t\txb-xa = ", signif(xb-xa,digs),"\t\tQ = ", ifelse(Q@positive,"+","-"),"exp(",Q@x,")\n",sep="")
  } else
  {
    Q <- .adaptsimstp.brob(f, xa, m, fa, fml, fm, ab, TRACE) + 
         .adaptsimstp.brob(f, m, xb, fm, fmr, fb, ab, TRACE)
  }
return(Q)
}
#########################################################################


#########################################################################
# pracma:::trapz
#
#############################
trapz.brob <- function (x, y) 
{
  if (missing(y)) {
    if (length(x) == 0) 
      return(0)
    y <- x
    x <- seq(along = x)
  }
  if (length(x) == 0 && length(y) == 0) 
    return(0)
  if (!(is.numeric(x) || is.complex(x)) || !(attr(y,"class") == "brob") || is.complex(y)) #???????how handle complex check?????
    stop("Arguments 'x' and 'y' must be real or complex vectors.")
  m <- length(x)
  if (length(y) != m) 
    stop("Arguments 'x', 'y' must be vectors of the same length.")
  if (m <= 1) 
    return(0)
  xp <- c(x, x[m:1])
  #brob
  yp  <- cbrob(as.brob(numeric(m)), y[m:1])
  n <- 2 * m
  p1 <- sum(xp[1:(n - 1)] * yp[2:n]) + xp[n] * yp[1]
  p2 <- sum(xp[2:n] * yp[1:(n - 1)]) + xp[1] * yp[n]
  return(0.5 * (p1 - p2))
}
#########################################################################


#########################################################################
# pracma:::quadgr
#
#######################################################################
quadgr.brob <- function (f, a, b, tol = .Machine$double.eps^(1/2), ...) 
{
  stopifnot(is.numeric(a), length(a) == 1, is.numeric(b), length(b) == 1)
  fun <- match.fun(f)
  f <- function(x) fun(x, ...)
  if (length(f(a)) != 1 || length(f(b)) != 1) 
    stop("Function 'f' may be multi-valued; please check.")
  if (length(f(c(a, b))) != 2) 
    stop("Function 'f' needs to be vectorized (see 'Vectorize').")
  if (a == b) 
  {
    return(list(value = 0, rel.err = 0))
  } else
  if (a > b) {
    tmp <- a
    a <- b
    b <- tmp
    rev_p <- TRUE
  } else
  {
    rev_p <- FALSE
  }  
  if (is.infinite(a) || is.infinite(b)) {
    if (is.finite(a) && is.infinite(b)) {
      f1 <- function(t) f(a + t/(1 - t))/(1 - t)^2
      Q <- quadgr.brob(f1, 0, 1, tol = tol)
    }
    else if (is.infinite(a) && is.finite(b)) {
      f2 <- function(t) f(b + t/(1 + t))/(1 + t)^2
      Q <- quadgr.brob(f2, -1, 0, tol = tol)
    }
    else if (is.infinite(a) && is.infinite(b)) {
      f1 <- function(t) f(t/(1 - t))/(1 - t)^2
      f2 <- function(t) f(t/(1 + t))/(1 + t)^2
      Q1 <- quadgr.brob(f1, 0, 1, tol = tol/2)
      Q2 <- quadgr.brob(f2, -1, 0, tol = tol/2)
      Q <- list(value = Q1$value + Q2$value, rel.err = Q1$rel.err + 
                  Q2$rel.err)
    }
    if (rev_p) 
      Q$value <- -Q$value
    return(Q)
  }
  xq <- c(0.125233408511469, 0.36783149899818, 0.587317954286617, 
          0.769902674194305, 0.904117256370475, 0.981560634246719)
  wq <- c(0.249147045813403, 0.233492536538355, 0.203167426723066, 
          0.160078328543346, 0.106939325995318, 0.0471753363865118)
  xq <- matrix(c(xq, -xq), ncol = 1)
  wq <- c(wq, wq)
  nq <- length(xq)
  maxit <- 17
  
  #Q0 <- zeros(maxit, 1)
  #Q1 <- zeros(maxit, 1)
  #Q2 <- zeros(maxit, 1)
  Q0 <- as.brob( rep(0,maxit) )
  Q1 <- as.brob( rep(0,maxit) )
  Q2 <- as.brob( rep(0,maxit) )
  
  hh <- (b - a)/2
  x <- (a + b)/2 + hh * xq
  
  #brob
  fx <- lapply(as.list(x), function(x) fun(x))
  hhwq <- hh * wq
  Q0[1] <- scalarprod.brob(as.brob(hhwq), list2vec.brob(fx))
  #Q0[1] <- hh * wq %*% fun(x)
  
  for (k in 3:maxit) {
    hh <- hh/2
    x <- cbind(x + a, x + b)/2
    #
    A <- as.brob(numeric(12))
    for (i in 1:12) A[i] <- sum(f(x[i, ]))

    #brob
    hhwq <- hh * wq
    #Q0[k] <- hh*wq %*% A
    Q0[k] <- scalarprod.brob(hhwq, A)
    
    if (k >= 5) {
      Q1[k] <- .rich.brob(Q0, k)
      Q2[k] <- .rich.brob(Q1, k)
    } else
    if (k >= 3)
    {
      Q1[k] <- .rich.brob(Q0, k)
    }
    if (k >= 6) {
      Qv <- c(Q0[k], Q1[k], Q2[k])
      Qw <- c(Q0[k - 1], Q1[k - 1], Q2[k - 1])
    } else 
    if (k >= 4) {
      Qv <- c(Q0[k], Q1[k])
      Qw <- c(Q0[k - 1], Q1[k - 1])
    } else
    {
      Qv <- Q0[k]
      Qw <- Q0[k - 1]
    }
    
    #brob
    if(is.list(Qv)) Qv <- list2vec.brob(Qv)
    if(is.list(Qw)) Qw <- list2vec.brob(Qw)
    absQvminusQw <- abs(Qv - Qw)
    #print(absQvminusQw)
    err <- min(absQvminusQw)

    #j <- which.min(abs(Qv - Qw))
    j <- which(err@x %in% absQvminusQw@x)
    
    Q <- Qv[j]
    if ( err < tol || (!is.finite(is.numeric(Q)) & is.numeric(Q) != 0) ) break
  }
  if (rev_p) Q <- -Q
return(list(value = Q, rel.err = err))
}

#pracma:::.rich
#associated to quadgr - just works, only renamed here...
.rich.brob <- function (Q, k) 
{
  if (Q[k] != Q[k - 1]) {
    cc <- (Q[k - 1] - Q[k - 2])/(Q[k] - Q[k - 1]) - 1
  } else
  {
    cc <- 1
  }
  cc <- max(cc, 0.07)
  return(Q[k] + (Q[k] - Q[k - 1])/cc)
}
#########################################################################