######################################################################
#Author: <blinded>
#Date:   15 November 2004
#
#Code segment implementing various tools in connection with
#hypergeometric functions.
#
#This includes:
#  * Rising factorials (aka. Pochhammer symbols) and their log version
#  * Barnes extended hypergeometric form
#  * Stirling numbers of the first kind (signed/unsigned version)
######################################################################

######################################################################
# The Pochhammer symbol (x)^(k) is the rising factorial
#  (x+k-1)!/(x-1)! = (x+k-1) \cdot x = gamma(x+n)/gamma(x)
# Negative values are handled by the following identity:
#  (x)^((k)) = (-x)_{(k)} * (-1)^k => (-x)^(k) = (x)_{(k)}*(-1)^k
#
# Params:
#  x - base term (has to be integer)
#  k - also needs to be integer
#
# Returns:
#  (x)^k
######################################################################

#pochhammer <- function(x,k) {
  #return(prod( x:(x+k-1)))
#   return(gamma(x+k)/gamma(x))
#}

##log version - works only if x>0
#lpochhammer <- function(x,k) {
#  return(lgamma(x+k)-lgamma(x))
#}


##Call it rising factorials instead. This version is ok with neg vals
risef <- function(x,k) {
  ifelse(x<0, fallf(-x,k)*(-1)^k, gamma(x+k)/gamma(x))
}

##log version, only works when the risef is positive!!
lrisef <- function(x,k) {
  ifelse(x<0, lfallf(-x,k)+log((-1)^k), lgamma(x+k)-lgamma(x))
}

######################################################################
#Falling factorial function x!/(x-k)!
#
#Params:
# x - Base
# k - modifying index
######################################################################
fallf <- function(x,k) {
  ifelse(x<0, risef(-x,k)*(-1)^k, exp(lfactorial(x)-lfactorial(x-k)))
}

lfallf <- function(x,k) {
  ifelse(x<0, lrisef(-x,k)+log((-1)^k), lfactorial(x)-lfactorial(x-k))
}


######################################################################
#Coarse approximation of the hypergeometric function
#by just reducing the infinite sum to a sum with maxk terms
#
#Parameters:
# maxk -- approximate sum by how many terms?
#
#Note:
# all elements in a,b have to be positive!!
######################################################################
hypergeometric <- function(a,b,z,maxk=2000) {
  sum <- 0
  negab <- (min(c(sign(a),sign(b))) == -1)
  ##Coarse (when z~1) sum approximation of hypergeometric function.
  for (k in 0:maxk) {
    ##Operate on log scale as long as possible.
    if (!is.null(a)) {
      t1 <- sum(sapply(a,function(ai) { lrisef(ai,k)}))
    } else {
      t1 <- 0
    }
    if (!is.null(b)) {
      t2 <- sum(sapply(b,function(bi) { lrisef(bi,k)}))
    } else {
      t2 <- 0
    }
    #cat("k = ",k,"\n")
    #print(t1)
    #print(t2)
    sum <- sum + exp(t1+k*log(z)-lfactorial(k)-t2)
  }
 
  return(sum)
}


######################################################################
#Stirling numbers of first kind (signed or version)
#Algorithm by Bill Venables
#http://www.biostat.wustl.edu/archives/html/s-news/2003-08/msg00022.html
#
# Basically this algorithm starts by computing:
#  [0,0],...,[0,n] and then uses the recursion formula
#  [n+1,k] = [n,k-1] - n [n,k] to calculate
#  [1,0],...,[1,n]  until [n,0],...[n,n] by single steps
#
# Note that [0,k]=0, [n,0]=0, but [0,0]=1
#
# Parameters
#   n - calculate the vector [n,1],...,[n,n]
#   signed - boolean indicating whether we want the signed or unsigned version
#
# Returns: 
#   Vector containing ([n,1],...,[n,n])
######################################################################

stirling1 <- function(n,signed=T) {
  v <- c(1., 0)
  if(n <= 1) return(rev(v)[-1])
  for(j in 1:(n - 1))
    v <- c(v, 0) - j * c(0, v)
  #Above gives [n,n],...,[n,0]
  stirlingnumbers <- rev(abs(v))[-1]

  #Compute signed version?
  if (signed) { stirlingnumbers <- (-1)^(n-1:n) *   stirlingnumbers }

  return(stirlingnumbers)
  
}