library(arrangements)
.P <- function(m, n){
  sum(sapply(1L:m, function(i){
    sum(sapply(1L:min(m,n,i), function(k){
      npartitions(i, k)
    }))
  }))
}

.T <- function(alpha, a, b, kappa, i){
  if(length(kappa) == 0L || kappa[1L] == 0L) return(1)
  c <- -(i-1)/alpha + kappa[i] - 1
  d <- kappa[i]*alpha - i
  s <- seq_len(kappa[i]-1L)
  e <- d - s*alpha + vapply(s, function(i) sum(kappa >= i), integer(1L))
  g <- e + 1
  f <- kappa[seq_len(i-1)]*alpha - seq_len(i-1) - d
  h <- f + alpha
  l <- h*f
  prod(a+c)/prod(b+c) * prod((g-alpha)*e/(g*(e+alpha))) * prod((l-f)/(l+h))
}

.betaratio <- function(kappa, mu, k, alpha){
  t <- k - alpha*mu[k]
  s <- seq_len(k)
  u <- t + 1 - s + alpha*kappa[s]
  s <- seq_len(k-1L)
  v <- t - s + alpha*mu[s]
  s <- seq_len(mu[k]-1L)
  w <- vapply(s, function(i) sum(mu >= i), integer(1L)) - t - alpha*s
  alpha * prod(u/(u+alpha-1)) * prod((v+alpha)/v) * prod((w+alpha)/w)
}

HypergeoI <- function(m, alpha, a, b, n, x){
  summation <- function(i, z, j){
    for(kappai in seq_len(if(i==1L) j else min(kappa[i-1L],j))){
      kappa[i] <<- kappai
      t <- .T(alpha, a, b, kappa, i)
      z <- z * x * (n-i+1+alpha*(kappai-1L)) * t
      s <<- s + z
      if(j > kappai && i < n){
        summation(i+1L, z, j-kappai)
      }
    }
    kappa[i] <<- 0L
  }
  s <- 1
  kappa <- integer(0L)
  summation(1L, 1, m)
  s
}


Hypergeo <- function(m, alpha, a, b, x){
  if(all(x == x[1L])){
    return(HypergeoI(m, alpha, a, b, length(x), x[1L]))
  }
  #
  summation <- function(i, z, j){
    r <- if(i == 1L) j else min(kappa[i-1L],j)
    for(kappai in seq_len(r)){
      kappa[i] <<- kappai
      Nkappa <<- .Nkappa(kappa) - 1L
      z <- z * .T(alpha, a, b, kappa, i)
      if(Nkappa > 1L && (length(kappa) == 1L || kappa[2L] == 0L)){
        J[Nkappa,1L] <<- x[1L]*(1+alpha*(kappa[1L]-1L)) * J[Nkappa-1L,1L]
      }
      for(t in 2L:n) jack(0L, 1, 0L, t, kappa, Nkappa)
      s <<- s + z * J[Nkappa,n]
      if(j > kappai && i < n){
        summation(i+1L, z, j-kappai)
      }
    }
    kappa[i] <<- 0L
  }
  #
  jack <- function(k, beta, c, t, mu, Nmu){
    for(i in max(k,1L) : sum(mu > 0L)){
      if(length(mu)==i || mu[i] > mu[i+1L]){
        d <- Nmu
        gamma <- beta * .betaratio(kappa, mu, i, alpha)
        mu[i] <- mu[i]-1L
        Nmu <- .Nkappa(mu) - 1L
        if(mu[i] > 0L){
          jack(i, gamma, c+1L, t, mu, Nmu)
        }else{
          if(Nkappa > 1L){
            J[Nkappa,t] <<- J[Nkappa,t] + gamma*x[t]^(c+1L) *
              ifelse(any(mu>0L), J[Nmu,t-1L], 1)
          }
        }
        mu[i] <- mu[i] + 1L
        Nmu <- d
      }
    }
    if(k == 0L){
      if(Nkappa>1L) J[Nkappa,t] <<- J[Nkappa,t] + J[Nkappa,t-1L]
    }else{
      J[Nkappa,t] <<- J[Nkappa,t] + beta * x[t]^c * J[Nmu,t-1L]
    }
  }
  #
  n <- length(x)
  Pmn <- .P(m,n)
  s <- 1
  J <- matrix(0, Pmn, n)
  J[1L,] <- cumsum(x)
  kappa <- integer(0L)
  #
  D <- rep(NA_integer_, Pmn)
  Last <- t(as.integer(c(0,m,m)))
  fin <- 0L
  for(. in 1L:n){
    NewLast <- matrix(NA_integer_, nrow = 0L, ncol = 3L)
    for(i in 1L:nrow(Last)){
      manque <- Last[i,2L]
      l <- min(manque, Last[i,3L])
      if(l > 0L){
        D[Last[i,1L]+1L] <- fin + 1L
        s <- 1L:l
        NewLast <- rbind(NewLast, cbind(fin+s, manque-s, s))
        fin <- fin + l
      }
    }
    Last <- NewLast
  }
  .Nkappa <- function(kappa){
    kappa <- kappa[kappa>0L]
    if(length(kappa)==0L) return(1L)
    D[.Nkappa(head(kappa,-1L))] + tail(kappa,1L)
  }
  summation(1L, 1, m)
  s
}

alpha <- 2
x <- c(3,2)
Hypergeo(3, alpha, c(4,2), 5, x)
CharFunToolR::HypergeompFqMat(t(c(4,2)), 5, x, MAX=3)

x <- c(2,2,2)
Hypergeo(3, alpha, c(4,2), 5, x)
CharFunToolR::HypergeompFqMat(t(c(4,2)), 5, x, MAX=3)