R/tweedie.R

In tweedie: Evaluation of Tweedie Exponential Family Models

#############################################################################
ptweedie.inversion <- function(q, mu, phi,  power, exact = FALSE ){ 
  # Peter K Dunn 
  # 08 Feb 2000--2005
  
  y <- q
  cdf <- array( dim = length(y) )
  # Error checks
  if ( power < 1) stop("power must be greater than 1.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y < 0) ) stop("y must be a non-negative vector.")
  if ( any(mu <= 0) ) stop("mu must be positive.")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.")
  } else {
    mu <- array( dim = length(y), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  
  #its <- y
  for (i in (1:length(y))) {
    
    # This has been added to avoid an issue with *very* small values of y
    # causing the FORTRAN to die when p>2; 
    # reported by Johann Cuenin 09 July 2013 (earlier, but that was the easiest email I could find about it :->)
    if ( ( power > 2 ) & (y[i] < 1.0e-300) ) {
      ### THE  e-300  IS ARBITRARY!!!!  ###
      # Keep an eye on it; perhaps it needs changing
      
      # That is, y is very small: Use the limit as y->0 as the answer as FORTRAN has difficulty converging
      # I have kept the call to the FORTRAN (and for p>2, of course, the limiting value is 0).
      #  I could have done this differently
      # by redefining very small y as y=0.... but this is better methinks
      cdf[i] <- 0
    } else {
      
      #   cat(power,phi[i], y[i], mu[i], as.integer(exact),"\n")
      tmp <- .Fortran( "twcdf",
                       as.double(power),
                       as.double(phi[i]),
                       as.double(y[i]),
                       as.double(mu[i]),
                       as.integer( exact ),
                       as.double(0), # funvalue
                       as.integer(0), # exitstatus
                       as.double(0), # relerr
                       as.integer(0)) # its
#                       PACKAGE="tweedie")
      
      cdf[i] <- tmp[[6]]
    }
  }
  
  cdf
  
}



#############################################################################
dtweedie.dldphi.saddle <- function(phi, mu, power, y){
  # Calculates the derivative of log f wrt phi
  # where the density is the saddlepoint density
  
  # Peter Dunn
  # 13 August 2002
  
  dev <- tweedie.dev( power = power, 
                      y = y, 
                      mu = mu)
  l <-  (-1) / (2 * phi) + dev / (2 * phi ^ 2)
  
  -2 * sum(l)
}







#############################################################################
dtweedie.logl <- function(phi, y, mu, power) {
  # Computes the log-likelihood for
  # a Tweedie density.  
  
  # Peter Dunn
  # 26 April 2001
  
  sum( log( dtweedie( y = y, 
                      mu = mu, 
                      phi = phi, 
                      power = power) ) )
  
}





#############################################################################
logLiktweedie <- function(glm.obj, dispersion = NULL) {
  # Computes the log-likelihood for
  # a Tweedie glm.  
  
  # Peter Dunn
  # 19 October 2017
  
  p <- get("p", envir = environment(glm.obj$family$variance))
  if (p == 1) message("*** Tweedie index power = 1: Consider using  dispersion = 1  in call to  logLiktweedie().\n")
  
  AICtweedie(glm.obj, 
             dispersion = dispersion, 
             k = 0, 
             verbose = FALSE) / (-2)
}


#############################################################################
dtweedie.logl.saddle <- function( phi, power, y, mu, eps=0){
  # Calculates the log likelihood of Tweedie densities
  # where the density is the saddlepoint density
  
  # Peter Dunn
  # 01 May 2001
  sum( log( dtweedie.saddle(power = power, 
                            phi = phi, 
                            y = y, 
                            mu = mu,
                            eps = eps) ) )
  
}


#############################################################################
dtweedie.logv.bigp <- dtweedie.logv.bigp <- function( y, phi, power){ 
  # Peter K Dunn 
  # 02 Feb 2000 
  # 
  
  #
  # Error traps
  #
  
  if ( power < 2) stop("power must be greater than 2.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y <= 0) ) stop("y must be a strictly positive vector.")
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  
  #
  # Set up
  #
  
  p <- power
  a <- ( 2 - p ) / ( 1 - p )
  a1 <- 1 - a
  
  r <- -a1 * log(phi) - log(p - 2) - a * log(y) + a * log(p - 1)
  drop <- 37
  
  #
  # Now we find limits of summation, using Stirling's approximation
  # to approximate the gamma terms, and find max as k.max.
  #
  
  logz <- max(r)
  k.max <- max( y ^ (2 - p) / ( phi * (p - 2) ) )
  k <- max( 1, k.max )
  
  c <- logz + a1 + a * log(a)
  vmax <- k.max * a1
  estlogv <- vmax
  #
  # Now we search either side for when we can
  # ignore terms as they are negligible.
  # Remember we are using the envelope as the tests,
  # not individual terms in the series.
  
  # First:  the upper limit of k
  
  while ( estlogv > (vmax - drop) ) {
    k <- k + 2
    estlogv <- k*( c - a1 * log(k) )
  }
  
  hi.k <- ceiling(k)
  #
  # Now the lower limit of k
  #
  logz <- min(r)
  k.max <- min( y ^ (2 - p) / ( phi * (p - 2) ) )
  k <- max( 1, k.max )
  c <- logz + a1 + a * log(a)
  vmax <- k.max * a1
  estlogv <- vmax
  
  while ( (estlogv > (vmax - drop) ) && ( k >= 2) ) {
    k <- max(1, k - 2)
    estlogv <- k*( c - a1 * log(k) )
  }
  
  lo.k <- max(1, floor(k) )
  
  #
  # Now sum the series between established limits.
  # We ensure it works for vector y.
  #
  k <- seq(lo.k, hi.k) 
  o <- matrix( 1, nrow = length(y)) 
  
  g <- matrix( lgamma( 1 + a * k) - lgamma(1 + k),
               nrow = 1, 
               ncol = length(k) )
  
  og <- o %*% g
  A <- outer(r, k) + og
  C <- matrix( sin( -a * pi * k ) * (-1)^k,
               nrow = 1, 
               ncol = length(k) )
  
  C <- o %*% C
  m <- apply(A, 1, max)
  ve <- exp(A - m)
  sum.ve <- apply( ve * C, 1, sum )
  
  # Now be careful!  Because of the +/- nature of the sin term,
  # sum.ve can be very small but negative  due to subtractive
  # cancellation.  Treat those carefully and separately.
  
  neg.sum.ve <- (sum.ve <= 0)
  pos.sum.ve <- (sum.ve > 0)
  
  logv <- sum.ve
  sum.ve[neg.sum.ve] <- 0
  logv[neg.sum.ve] <- -Inf
  logv[pos.sum.ve] <- log( sum.ve[pos.sum.ve] ) + m[pos.sum.ve]
  
  list(lo = lo.k, 
       hi = hi.k, 
       logv = logv, 
       k.max = k.max )
  
}


#############################################################################
dtweedie.logw.smallp <- function(y, phi, power){ 
  #
  # Peter K Dunn
  # 02 Feb 2000
  #
  
  #
  # Error traps
  #
  
  if ( power < 1) stop("power must be between 1 and 2.")
  if ( power > 2) stop("power must be between 1 and 2.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y <= 0) ) stop("y must be a strictly positive vector.")
  
  #
  # Set up
  #
  p <- power
  a <- ( 2 - p ) / ( 1 - p )         # Note that a<0 for 1<p<2
  
  a1 <- 1 - a
  r <- -a * log(y) + a * log(p - 1) - a1 * log(phi) - log(2 - p)        # All terms to power j
  
  drop <- 37             # Accuracy of terms: exp(-37)
  #
  # Find limits of summation using Stirling's approximation
  # to approximate gamma terms, and find max as j.max
  #
  logz <- max(r)             # To find largest  j  needed
  j.max <- max( y ^ ( 2 - p ) / ( phi * (2 - p) ) )
  j <- max( 1, j.max )
  
  cc <- logz + a1 + a * log(-a)    #
  # This is all the terms to j-power, not needing any other j terms.
  # The other terms are introduced when we know the values of j
  
  wmax <- a1 * j.max
  estlogw <- wmax
  
  # First, the upper limit of j
  while(estlogw > (wmax - drop) ){
    j <- j + 2
    estlogw <- j * (cc - a1 * log(j))
  }
  
  hi.j <- ceiling(j)
  
  #
  #
  # Now the lower limit of j
  #
  #
  logz <- min(r) 
  j.max <- min( y ^ (2 - power) / ( phi * (2 - power) ) )
  
  j <- max( 1, j.max)
  wmax <- a1 * j.max 
  estlogw <- wmax 
  
  # First, optimize to find the location of the maximum
  while ( ( estlogw > (wmax - drop) ) && ( j >= 2) ) {
    j <- max(1, j - 2)
    oldestlogw <- estlogw
    estlogw <- j * (cc - a1 * log(j))
  }
  
  lo.j <- max(1, floor(j))
  # 
  # Now sum the series between established limits.
  # We ensure it works for vector y.
  
  j <- seq( lo.j, hi.j)      # sequence of j terms needed
  
  o <- matrix( 1, nrow = length(y))     # matrix of ones
  
  g <- matrix(lgamma( j + 1 ) + lgamma( -a * j ),
              nrow = 1, 
              ncol = hi.j - lo.j + 1)
  
  og <- o %*% g              # matrix of gamma terms
  A <- outer(r, j) - og      # the series, almost ready to sum
  m <- apply(A, 1, max)      # avoid overflow; find maximum values
  we <- exp( A - m )         # evaluate terms, less max.
  sum.we <- apply( we, 1, sum) # sum terms
  logw <- log( sum.we ) + m  # now restore max.
  
  list(lo = lo.j, 
       hi = hi.j, 
       logw = logw, 
       j.max = j.max )
  
}


#############################################################################
dtweedie <- function(y, xi = NULL, mu, phi, power = NULL)
{
  #
  # This is a function for determining Tweedie densities.
  # Two methods are employed:  cgf inversion (type=1)
  # and series evaluation (type=2 if 1<p<2; type=3 if p>2)).
  # This function uses bivariate interpolation to accurately
  # approximate the inversion in conjunction with the series.
  #
  #
  # FIRST, establish whether the series or the cgf
  # # inversion is the better method.
  #
  
  #
  # Here is a summary of what happens:
  #
  #   p                          |  Whats happpens
  # -----------------------------+-------------------
  #  1 -> p.lowest (=1.3)        |  Use series A
  #  p.lowest (=1.3) -> p=2      |  Use series A if
  #                              |        xix > xix.smallp (=0.8)
  #                              |      inversion with rho if
  #                              |        xix < xix.smallp (=0.8)
  #  p=2 -> p=3                  |  Use series B if
  #                              |        xix > xix.smallp (=0.8)
  #                              |      inversion with rho if
  #                              |        xix < xix.smallp (=0.8)
  #  p=3 -> p.highest (=4)       |  Use series B if xix > xix.bigp (=0.9)
  #                              |      inversion with 1/rho if
  #                              |        xix < xix.bigp (=0.9)
  #  p > 4                       |  Use series B if xix > 0.03
  #                              | Saddlepoint approximation if xix<= 0.03 (NEEDS FIXING)
  #
  
  # Sort out the xi/power notation
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  
  # Error checks
  if ( any(power < 1) )  stop( paste(index.par.long, "must be greater than 1.\n") )
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y < 0) ) stop("y must be a non-negative vector.\n")
  if ( any(mu <= 0) ) stop("mu must be positive.\n")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.\n")
  } else {
    mu <- array( dim = length(y), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.\n")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  density <- y
  
  # Special Cases
  if ( power == 3 ){
    density <- statmod::dinvgauss(x = y, 
                                  mean = mu, 
                                  dispersion = phi)
    return(density)
  }
  if ( power == 2 ) {
    density <- dgamma( rate = 1 / (phi * mu), 
                       shape = 1 / phi, 
                       x = y )
    return(density)
  }
  if ( power == 0) {
    density <- dnorm( mean = mu, 
                      sd = sqrt(phi), 
                      x = y )
    return(density)
  }
  if ( (power == 1) & (all(phi == 1))) {
    # Poisson case
    density <- dpois(x = y / phi, 
                     lambda = mu / phi )
    return(density)
  }
  
  # Set up
  id.type0 <- array( dim = length(y) )
  id.series <- id.type0
  id.interp <- id.type0
  
  ###   Now consider the cases 1<p<2   ###
  id.type0 <- (y == 0)
  if (any(id.type0)) {
    if ( power > 2 ) {
      density[id.type0] <- 0
    } else {
      lambda <- mu[id.type0] ^ (2 - power) / (phi[id.type0] * (2 - power))
      density[id.type0] <- exp( -lambda )
    }
  }
  
  
  xi <- array( dim = length(y) )
  xi[id.type0] <- 0
  xi[!id.type0] <- phi[!id.type0] * y[!id.type0] ^ (power - 2)
  xix <- xi / ( 1 + xi )
  
  if ( (power > 1) && (power <= 1.1) ) {
    id.series <- (!id.type0)
    if (any(id.series)){
      density[id.series] <- dtweedie.series(y = y[id.series],
                                            mu = mu[id.series], 
                                            phi = phi[id.series],
                                            power = power)
    }
    return(density = density)
  }
  
  if ( power == 1 ) { # AND phi not equal to one here
    id.series <- rep(TRUE, length(id.series))
    id.interp <- rep(FALSE, length(id.series))
  }
  if ( (power > 1.1) && (power <= 1.2) ) {
    id.interp <- ( (xix > 0) & (xix < 0.1) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 1.1
      p.hi <- 1.2
      xix.lo <- 0
      xix.hi <- 0.1
      np <- 15
      nx <- 25
    }
  }
  
  if ( (power > 1.2) && (power <= 1.3) ) {
    id.interp <- ( (xix > 0) & (xix < 0.3) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 1.2
      p.hi <- 1.3
      xix.lo <- 0
      xix.hi <- 0.3
      np <- 15
      nx <- 25
    }
  }
  if ( (power > 1.3) && (power <= 1.4) ) {
    id.interp <- ( (xix > 0) & (xix < 0.5) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 1.3
      p.hi <- 1.4
      xix.lo <- 0
      xix.hi <- 0.5
      np <- 15
      nx <- 25
    }
  }
  if ( (power > 1.4) && (power <= 1.5) ) {
    id.interp <- ( (xix > 0) & (xix < 0.8) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 1.4
      p.hi <- 1.5
      xix.lo <- 0
      xix.hi <- 0.8
      np <- 15
      nx <- 25
    }
  }
  if ( (power > 1.5) && (power < 2) ) {
    id.interp <- ( (xix > 0) & (xix < 0.9) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 1.5
      p.hi <- 2
      xix.lo <- 0
      xix.hi <- 0.9
      np <- 15
      nx <- 25
    }
  }
  
  ### Cases p>2   ###
  if ( (power > 2) && (power < 3) ) {
    id.interp <- ( (xix > 0) & (xix < 0.9) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 2
      p.hi <- 3
      xix.lo <- 0
      xix.hi <- 0.9
      np <- 15
      nx <- 25
    }
  }
  if ( (power >= 3) && (power < 4) ) {
    id.interp <- ( (xix > 0) & (xix < 0.9) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 3
      p.hi <- 4
      xix.lo <- 0
      xix.hi <- 0.9
      np <- 15
      nx <- 25
    }
  }
  if ( (power >= 4) && (power < 5) ) {
    id.interp <- ( (xix > 0) & (xix < 0.9) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 4
      p.hi <- 5
      xix.lo <- 0
      xix.hi <- 0.9
      np <- 15
      nx <- 25
    }
  }
  if ( (power >= 5) && (power < 7) ) {
    id.interp <- ( (xix > 0) & (xix < 0.5) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 5
      p.hi <- 7
      xix.lo <- 0
      xix.hi <- 0.5
      np <- 15
      nx <- 25
    }
  }
  if ( (power >= 7) && (power <= 10) ) {
    id.interp <- ( (xix > 0) & (xix < 0.3) )
    id.series <- (!(id.interp | id.type0))
    if ( any(id.interp)) {
      grid <- stored.grids(power)
      p.lo <- 7
      p.hi <- 10
      xix.lo <- 0
      xix.hi <- 0.3
      np <- 15
      nx <- 25
    }
  }
  
  if ( power > 10) {
    id.series <- (y != 0)
    id.interp <- (!(id.series | id.type0))
  }
  
  if (any(id.series)) {
    density[id.series] <- dtweedie.series(y = y[id.series],
                                          mu = mu[id.series], 
                                          phi = phi[id.series],
                                          power = power)
  }
  
  if (any(id.interp)) {
    dim( grid ) <- c( nx + 1, np + 1 )
    rho <- dtweedie.interp(grid, 
                           np = np, 
                           nx = nx,
                           xix.lo = xix.lo, 
                           xix.hi = xix.hi,
                           p.lo = p.lo, 
                           p.hi = p.hi,
                           power = power, 
                           xix = xix[id.interp] )
    dev <- tweedie.dev(power = power, 
                       mu = mu[id.interp], 
                       y = y[id.interp])
    front <- rho / (y[id.interp] * sqrt(2 * pi * xi[id.interp]))
    density[id.interp] <- front * exp(-1/(2 * phi[id.interp]) * dev)
  }
  
  
  density
  
}


#############################################################################
dtweedie.saddle <- function(y, xi=NULL, mu, phi, eps=1/6, power=NULL) {
  #
  # Peter K Dunn
  # 09 Jan 2002
  #
  #
  
  # Sort out the xi/power notation
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  
  # Error traps
  #
  if( any(phi <= 0) )
    stop("phi must be positive.")
  if( (power >= 1) & any(y < 0))
    stop("y must be a non-negative vector.")
  if( any(mu <= 0) )
    stop("mu must be positive.")
  if ( length(mu) == 1 )
    mu <- array( dim=length(y), mu )
  if ( length(phi) == 1 )
    phi <- array( dim=length(y), phi )
  #
  # END error checks
  #
  #   if ( any(y==0) ) y <- y + (1/6)
  y.eps <- y 
  
  if (power<2) y.eps <- y + eps
  # Nelder and Pregibon's idea for problems at y=0
  
  y0 <- (y == 0)
  density <- y
  #    if(any(y == 0)) {
  #        if(power >= 2) {
  #         # This is a problem!
  #            density[y0] <- 0 * y[y0]
  #        }
  #      else{
  
  dev <- tweedie.dev(y = y, 
                     mu = mu, 
                     power = power)
  density <- (2 * pi * phi * y.eps ^ power) ^ (-1/2) * exp( -dev / (2 * phi) )
  
  #      }
  #   }
  if ( any(y == 0) ){
    if((power >= 1) && (power < 2)) {
      lambda <- mu[y0] ^ (2 - power) / (phi[y0] * (2 - power))
      density[y0] <- exp( -lambda)
    } else {
      density[y0] <- 0
    }
  }
  
  #    if(any(y != 0)) {
  #     dev <- tweedie.dev(y=y[yp], mu=mu[yp], power=power)
  #      density[yp] <- (2*pi*phi[yp]*y[yp]^power)^(-1/2) * exp( -dev/(2*phi[yp]) )
  #    }
  density
}


#############################################################################
dtweedie.series.bigp <- function(power, y, mu, phi){ 
  
  # 
  # Peter K Dunn 
  # 02 Feb 2000 
  # 
  
  #
  # Error traps
  #
  
  if ( power < 2) stop("power must be greater than 2.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y <= 0) ) stop("y must be a strictly positive vector.")
  if ( any(mu <= 0) ) stop("mu must be positive.")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.")
  } else {
    mu <- array( dim = length(y), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  
  result <- dtweedie.logv.bigp(power = power, 
                               y = y, 
                               phi = phi)
  logv <- result$logv
  
  theta <- mu ^ (1 - power) / ( 1 - power )
  kappa <- mu ^ (2 - power) / ( 2 - power )
  
  logfnew <- (y * theta - kappa) / phi - log( pi * y) + logv
  f <- exp( logfnew )
  
  list(density = f, 
       logv = logv, 
       lo = result$lo, 
       hi = result$hi )
  
}


#############################################################################
dtweedie.series <- function(y, power, mu, phi){ 
  #
  # Peter K Dunn
  # 09 Jan 2002
  #
  
  #
  # Error traps
  #
  
  if ( power < 1) stop("power must be between 1 and 2.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y < 0) ) stop("y must be a non-negative vector.")
  if ( any(mu <= 0) ) stop("mu must be positive.")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.")
  } else {
    mu <- array( dim = length(y), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  
  y0 <- (y == 0 )
  yp <- ( y!=0 )
  density <- array( dim = length(y))
  
  if ( (power == 2) | (power == 1) ) { # Special cases
    if ( power == 2 ){
      density <- dgamma( y, 
                         shape = 1 / phi, 
                         rate = 1 / (phi * mu ) )
    }
    
    if ( (power == 1) ){
      density <- dpois(x = y / phi, 
                       lambda = mu / phi ) / phi  # Using identity: f(y; mu, phi) = c f(cy; c mu, c phi) for p=1.  Now set c = 1/phi
      if ( !all(phi == 1)){
        warnings("The density computations when phi=1 may be subject to errors using floating-point arithmetic\n")
      }
    }
    # CHANGED: 31 October, thanks to Dina Farkas' email 18 October 2012.  WAS:
    #   if ( (power == 1) & (all(phi==1)) ){
    #    density <- dpois(x=y/phi, lambda=mu/phi )
    #   }
  } else{
    
    if ( any(y == 0) ) {
      if ( power > 2 ) {
        density[y0] <- 0 * y[y0]
      }
      if ( ( power > 1) && (power < 2) ) {
        lambda <- mu[y0] ^ (2 - power) / ( phi[y0] * (2 - power) )
        density[y0] <- exp( -lambda )
      }
    }
    
    if ( any( y!= 0 ) ) { 
      if ( power > 2 ) {
        density[yp] <- dtweedie.series.bigp(power = power,
                                            mu = mu[yp], 
                                            y = y[yp],
                                            phi = phi[yp])$density
      }
      
      if ( ( power > 1 ) && ( power < 2 ) ) {
        density[yp] <- dtweedie.series.smallp(power = power,
                                              mu = mu[yp], 
                                              y = y[yp],
                                              phi = phi[yp])$density
      }
    }
  }
  
  density
  
}




#############################################################################
dtweedie.series.smallp <- function(power, y, mu, phi){ 
  
  #
  # Peter K Dunn
  # 02 Feb 2000
  #
  
  #
  # Error traps
  #
  
  if ( power < 1) stop("power must be between 1 and 2.")
  if ( power > 2) stop("power must be between 1 and 2.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y <= 0) & (power >= 2) ) stop("y must be a strictly positive vector.")
  if ( any(y < 0) & (power > 1) & (power < 2) ) stop("y must be a non-negative vector.")
  if ( any(mu <= 0) ) stop("mu must be positive.")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.")
  } else {
    mu <- array( dim = length(y), mu )
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
  }
  
  result <- dtweedie.logw.smallp( y = y, 
                                  power = power, 
                                  phi = phi)
  logw <- result$logw
  
  
  tau <- phi * (power - 1) * mu ^ ( power - 1 )
  lambda <- mu ^ (2 - power) / ( phi * (2 - power) )
  logf <- -y / tau - lambda - log(y) + logw
  f <- exp( logf )
  
  list(density = f, 
       logw = logw, 
       hi = result$hi, 
       lo = result$lo)
  
}




#############################################################################
ptweedie <- function(q, xi = NULL, mu, phi, power = NULL) {
  # Evaluates the cdf for
  # Tweedie distributions.
  
  # Peter Dunn
  # 01 May 2001
  
  # Sort out the xi/power notation
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  
  y <- q
  
  y.negative <- (y < 0)
  y.original <- y
  y <- y[ !y.negative ]
  
  # Error checks
  if ( any(power < 1) )  stop( paste(index.par.long, "must be greater than 1.\n") )
  if ( any(phi <= 0) ) stop("phi must be positive.")
  # if ( any(y < 0) ) stop("y must be a non-negative vector.\n")
  if ( any(mu <= 0) ) stop("mu must be positive.\n")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y.original) ) stop("mu must be scalar, or the same length as y.\n")
  } else {
    mu <- array( dim = length(y.original), mu )
    # A vector of all mu's
  }
  mu.original <- mu
  mu <- mu[ !y.negative ]
  
  if ( length(phi) > 1) {
    if ( length(phi) != length(y.original) ) stop("phi must be scalar, or the same length as y.\n")
  } else {
    phi <- array( dim = length(y.original), phi )
    # A vector of all phi's
  }
  phi.original <- phi
  phi <- phi[ !y.negative ]
  
  
  
  cdf.positives <- array( dim = length(y) )
  cdf <- y.original
  
  # Special Cases
  if ( power == 2 ) {
    f <- pgamma( rate = 1 / (phi * mu), 
                 shape = 1 / phi, 
                 q = y )
  }
  if ( power == 0) {
    f <- pnorm( mean = mu, 
                sd = sqrt(phi),
                q = y )
  }
  if ( power == 1) {
    f <- ppois(q = y, 
               lambda = mu / phi )
  }
  
  # Now, for p>2 the only option is the inversion
  if ( power> 2 ) {
    
    # However, when y/q is very small, 
    # the function should return 0;
    # sometimes it fails (when *very* small...).
    # This adjustment is made in ptweedie.inversion()
    f <- ptweedie.inversion(power = power, 
                            mu = mu, 
                            q = y, 
                            phi = phi)
  }
  
  # Now for 1<p<2, the two options are the series or inversion.
  # We avoid the series when p is near one, otherwise it is fine.
  
  # A few caveats.  Gustaov noted this case:
  # ptweedie(q=7.709933e-308, mu=1.017691e+01, phi=4.550000e+00, power=1.980000e+00)
  # which fails for the inversion, but seems to go fine for the series.
  # So the criterion should be a bit more detailed that just p<1.7...
  # But what?
  # Shallow second derivative of integrand is OK in principle... but hard to ascertain.
  
  # In a bug report by Gustavo Lacerda (April 2017), 
  # it seems that the original code here was a bit too 
  # harsh on the series, and a bit forgiving on the inversion.
  # Changed April 2017 to use the series more often.
  
  # ### OLD CODE:
  # if ( (power>1) & (power<2) ) {
  #     if ( power <1.7 ) {
  #        f <- ptweedie.series(power=power, q=y, mu=mu, phi=phi )
  #     } else{
  #        f <- ptweedie.inversion( power=power, q=y, mu=mu, phi=phi)
  #     }
  # }
  
  ### REVISED CODE:
  ### The choice of  1.999 is arbitrary.  Probably needs serious attention to decide properly
  ### Changed early 2017; thanks to Gustavo Lacerda
  if ( (power > 1) & (power < 2) ) {
    
    if ( power < 1.999) {
      f <- ptweedie.series(power = power, 
                           q = y, 
                           mu = mu, 
                           phi = phi )
    } else{
      f <- ptweedie.inversion( power = power, 
                               q = y, 
                               mu = mu, 
                               phi = phi)
    }
  }
  cdf[ !y.negative ] <- f
  cdf[  y.negative ] <- 0
  
  cdf[ is.infinite( cdf ) ] <- 1
  
  cdf[ cdf > 1 ] <- 1
  
  return(cdf)
  
}



#############################################################################
ptweedie.series <- function(q, power, mu, phi) {
  # 
  # Peter K Dunn 
  # 14 March 2000 
  #
  y <- q
  # Error checks
  if ( power < 1) stop("power must be between 1 and 2.\n")
  if ( power > 2) stop("power must be between 1 and 2.\n")
  if ( any(phi <= 0)) stop("phi must be positive.\n")
  if ( any(y < 0) ) stop("y must be a non-negative vector.\n")
  if ( any(mu <= 0) ) stop("mu must be positive.\n")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.\n")
  } else {
    mu <- array( dim = length(y), mu )
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.\n")
  } else {
    phi <- array( dim = length(y), phi )
  }
  
  
  # Set up
  p <- power
  lambda <- mu ^ (2 - p) / ( phi * (2 - p) )
  tau <- phi * (p - 1) * mu ^ ( p - 1 )
  alpha <- (2 - p) / ( 1 - p )
  
  
  # Now find the limits on N:
  # First the lower limit on N
  
  drop <- 39
  lambda <- max(lambda )
  logfmax <-  -log(lambda)/2
  estlogf <- logfmax
  N <- max( lambda )
  
  while ( ( estlogf > (logfmax - drop) ) & ( N > 1 ) ) {
    N <- max(1, N - 2)
    estlogf <- -lambda + N * ( log(lambda) - log(N) + 1 ) - log(N)/2
  }
  
  lo.N <- max(1, floor(N) )
  
  
  # Now for the upper limit on N
  
  lambda <- mu ^ (2 - p) / ( phi * (2 - p) )
  lambda <- min( lambda )
  logfmax <-  -log(lambda) / 2
  estlogf <- logfmax
  
  N <- max( lambda )
  
  while ( estlogf > (logfmax - drop) ) {
    N <- N + 1
    estlogf <- -lambda + N * ( log(lambda) - log(N) + 1 ) - log(N)/2
  }
  
  hi.N <- max( ceiling(N) )
  
  # Now evaluate between limits of N
  
  cdf <- array( dim = length(y), 0 )
  
  lambda <- mu ^ (2 - p) / ( phi * (2 - p) )
  tau <- phi * (p - 1) * mu ^ ( p - 1 )
  alpha <- (2 - p) / ( 1 - p )
  
  
  for (N in (lo.N : hi.N)) {
    
    # Poisson density
    pois.den <- dpois( N, lambda)
    
    # Incomplete gamma
    incgamma.den <- pchisq(2 * y / tau, 
                           -2 * alpha * N )
    
    # What we want
    cdf <- cdf + pois.den * incgamma.den
    
  }
  
  cdf <- cdf + exp( -lambda )
  its <- hi.N - lo.N + 1
  
  
  cdf
  
}



#############################################################################
stored.grids <- function(power){
  # This S-Plus function contains all the stored interpolation
  # grids for interpolating the Tweedie densities.
  
  # Peter Dunn
  # 18 April 2001
  
  # Grids stored
  #
  ##### PART 1:  1 < p < 2 ####
  #
  # 1.5 < p < 2 (A)
  if ( (power > 1.5) && ( power < 2) ) {
    grid <- c(0.619552027046892,
              -1.99653020187375,
              -8.41322974185725,
              -36.5410153849965,
              -155.674714429765,
              -616.623539139617,
              -2200.89580607229,
              -6980.89597754589,
              -19605.6031284706,
              -48857.4648924699,
              -108604.369946264,
              -216832.50696005,
              -391832.053251182,
              -646019.900543419,
              -979606.689484813,
              -1377088.40156071,
              -1808605.4504829,
              -2236021.71690435,
              -2621471.84422418,
              -2935291.72910167,
              -3160933.2895598,
              -3296019.98597692,
              -3350191.70672175,
              -3341188.71720722,
              -3290375.55904795,
              -3256251.04779701,
              0.55083671442337,
              2.79545515380255,
              6.43661784251532,
              0.0593276289419893,
              -88.7776267014812,
              -592.599945216246,
              -2715.1932676716,
              -9995.06110621652,
              -31004.2502694815,
              -82973.5646780374,
              -194644.329882664,
              -405395.126655972,
              -758085.608783886,
              -1285893.38788744,
              -1997446.83841278,
              -2866954.62503587,
              -3834679.24827227,
              -4818512.83125468,
              -5732403.6885235,
              -6504847.15089854,
              -7091578.60837028,
              -7479819.10447519,
              -7684747.8312771,
              -7739400.44941064,
              -7651256.22196861,
              -10780871.374931,
              -0.296443816521139,
              -0.740400436722078,
              13.5437223703132,
              132.570369303337,
              739.235236905108,
              3222.27476321793,
              11883.2578358337,
              37985.4069102717,
              106385.512302156,
              263186.759184411,
              579709.939191383,
              1146196.24715222,
              2051166.63217944,
              3349844.60761646,
              5033457.4421439,
              7014240.88296282,
              9135328.47624901,
              11203377.8720252,
              13031732.9028836,
              14478751.8662425,
              15470353.6878961,
              16004512.3987063,
              16147291.8371304,
              16087058.36173,
              17256553.6861257,
              -90196970.3516856,
              0.212508829604868,
              -1.25345600957563,
              -33.6344002424504,
              -229.413615894842,
              -961.2927083918,
              -3048.22512983433,
              -7865.94887016698,
              -16619.1472772416,
              -27426.4597246056,
              -29170.0602809743,
              6910.62618199269,
              133529.230915424,
              425031.289018802,
              962079.672230075,
              1805724.1265787,
              2970428.54627254,
              4407924.65252939,
              6009208.15039934,
              7623903.3639327,
              9090617.14760774,
              10273366.8474699,
              11112760.9034168,
              11783847.5998635,
              13967546.0605232,
              40580571.9968577,
              -1800585369.9502,
              -0.117096128441307,
              3.86814002139823,
              52.1297859521411,
              231.959187763769,
              297.55942409085,
              -2028.44588042558,
              -16035.7900988595,
              -69727.4207742442,
              -231001.098781009,
              -634711.96272568,
              -1500867.48368527,
              -3120952.565006,
              -5795192.65274283,
              -9728268.26803975,
              -14923011.3657198,
              -21127733.2735346,
              -27877621.6155107,
              -34629531.4516162,
              -40942366.5228353,
              -46618585.271072,
              -51693745.2168584,
              -56076194.3321918,
              -57846478.5394169,
              -39555239.6466308,
              244583291.9647,
              -19387464175.4177,
              -0.0251908119117217,
              -6.97888536964718,
              -57.2053622635109,
              -7.59149018702844,
              1940.74883722613,
              13776.3994355061,
              59688.9506688118,
              198967.267063524,
              554664.181209208,
              1342679.76760535,
              2880718.85057692,
              5556233.09746906,
              9746453.11803942,
              15705233.9839578,
              23444873.7106245,
              32630213.4195304,
              42473186.150293,
              51604301.4143779,
              57950838.5609574,
              58814629.2970318,
              51670532.360853,
              37073253.8293895,
              31112689.7825927,
              166639899.499491,
              2285908471.65836,
              -142300999941.473,
              0.239462332088407,
              9.7006247266445,
              24.6796376740648,
              -621.977746808668,
              -6061.14667092794,
              -28846.5231563469,
              -92974.9043842097,
              -229420.907671922,
              -455750.587752974,
              -720917.683643727,
              -798268.858897086,
              -142553.959799974,
              2227938.90487962,
              7722331.28668364,
              17902944.5216494,
              33799508.0633365,
              54727373.5621622,
              76699837.6804652,
              90928852.1062459,
              83791368.1450994,
              41369883.2049391,
              -33644574.3006822,
              -51466205.6969009,
              728822186.855893,
              12613725489.8081,
              -779054232341.425,
              -0.531163176360929,
              -9.62133388925249,
              83.5858879354598,
              1749.74435647472,
              10388.878276731,
              29739.3251394829,
              27228.8352302813,
              -137459.132798184,
              -797275.579291472,
              -2522656.85139472,
              -6002370.79311501,
              -11564879.4570243,
              -18293067.3260108,
              -23137574.7865943,
              -20926214.0959888,
              -6386212.74814566,
              20947334.9587757,
              49084025.5298454,
              45674483.7892477,
              -43357709.595878,
              -274284042.960864,
              -630264236.584573,
              -707988499.827136,
              2768071791.46789,
              55323323631.9224,
              -3346981844559.16,
              0.843274592042774,
              1.85634928776026,
              -303.233463115877,
              -2964.77206695662,
              -8035.18601516459,
              26623.9135533587,
              296203.748735318,
              1300209.48867761,
              4002931.81502221,
              10154091.5309971,
              23011763.2727463,
              48484290.5041288,
              96405481.1482104,
              181157702.817844,
              319621265.640707,
              523128884.608782,
              780656497.356576,
              1034030709.78307,
              1153235261.31943,
              934435310.723166,
              172736683.699636,
              -1059296923.94426,
              -1259513524.21048,
              11540211151.7843,
              202982012848.938,
              -11765811129127.3,
              -0.971763767126055,
              20.7278527481909,
              599.216068887268,
              2508.62108313741,
              -15799.3976275259,
              -190836.337882638,
              -873333.216664437,
              -2412364.04114457,
              -4372947.5863571,
              -4280053.82275206,
              4581124.92368874,
              37692021.4517723,
              125690565.165227,
              317577060.569135,
              672918414.800248,
              1235540175.83018,
              1981789079.14886,
              2741889556.9801,
              3118003228.13883,
              2478792089.71153,
              201910662.627241,
              -3440015764.71965,
              -3701941985.52277,
              36624484882.9033,
              631340208742.919,
              -35016783929182,
              0.465205915130579,
              -62.4878848976034,
              -713.029716677298,
              3463.89852940917,
              76620.1536606951,
              425050.797698706,
              1098749.49892239,
              830042.192958006,
              -3609253.07165929,
              -14019647.104919,
              -22008231.8565183,
              2879363.57573307,
              127763068.732814,
              471678257.276824,
              1199156009.13846,
              2449442490.02104,
              4177372240.464,
              5937970631.74612,
              6701001745.26596,
              4850964750.68393,
              -1252857624.6872,
              -11008101341.7195,
              -11613916592.1703,
              99208893729.6641,
              1723996766015.73,
              -90686535774445.6,
              1.35055309427014,
              109.632652005803,
              -15.551388984273,
              -19110.1792293505,
              -151248.805892165,
              -326520.136900253,
              1323004.57105885,
              11263794.9628682,
              40182893.7043995,
              99473459.937714,
              213007606.45782,
              462071289.943091,
              1036322096.46886,
              2233228139.73472,
              4392262214.22542,
              7778537815.0104,
              12326926829.6142,
              17068944624.2195,
              19342475640.096,
              14633223559.5287,
              -1238504576.84042,
              -24780929364.4666,
              -19774890981.5309,
              267278604099.047,
              4262049717379.08,
              -209109831966511,
              -4.90254965773785,
              -106.16559587471,
              2473.32253262609,
              39971.9971896279,
              105932.466006557,
              -940508.183476386,
              -7972770.02131952,
              -26249511.2854385,
              -38016600.0226778,
              28952783.4114192,
              296603175.897931,
              921440155.664983,
              2211825458.1676,
              4873052470.30472,
              9959226055.61093,
              18095018536.2409,
              28333607390.4192,
              37736361000.0061,
              42033152625.5253,
              35930405200.4038,
              11674689966.7752,
              -38568497842.5504,
              -74929319913.8016,
              440819854424.984,
              9367147142753.47,
              -436516234054193,
              9.21654006976202,
              -53.4970979572424,
              -6659.50032530198,
              -38474.8203617517,
              310488.103734294,
              3845463.04045484,
              14525744.8285833,
              14191442.4974447,
              -71954549.6112557,
              -299246254.252719,
              -416148031.634689,
              348026958.617362,
              2939406987.4056,
              7971636742.33878,
              16399104338.5733,
              31135394894.2163,
              55589206993.0977,
              84773664582.7313,
              93124980530.7634,
              37204300253.4301,
              -100809540619.219,
              -211523088540.365,
              124648961402.01,
              2454987217464.61,
              23136880987249.1,
              -833737986871973,
              -10.6044400547848,
              457.68066225078,
              9572.011521514,
              -37331.5997372033,
              -1178378.58119297,
              -5650741.83842416,
              1857995.74642602,
              109904372.796869,
              458568718.1367,
              915187897.2047,
              811019172.09732,
              92558199.5987161,
              1998081069.83546,
              13128498122.3811,
              36243363277.9486,
              61424177454.5951,
              73667085976.0338,
              86031421495.4807,
              166947042118.509,
              384532819727.471,
              591661667761.69,
              78943846328.2111,
              -2527580012466.73,
              -6916515802327.55,
              21469457257461.8,
              -1.4991232719638e+015,
              3.17664940253603,
              -1029.70926107305,
              -5007.37894552437,
              206353.644207492,
              1762192.60826409,
              -794281.059186644,
              -59885497.7612809,
              -267968880.144224,
              -316522642.887633,
              1411438830.94685,
              7247583561.26104,
              15913831315.8407,
              19845848486.0143,
              17971388397.4505,
              42579150295.9815,
              155992918698.532,
              356278537036.23,
              414110925668.561,
              -161408422221.427,
              -1685255831372.5,
              -3094542673159.69,
              -53087024607.4896,
              16845038275179.6,
              63212141053606.3,
              202776604073488,
              -2.43199634134677e+015 )
  }
  
  # 1.4 < p < 1.5 (B)
  if ( (power > 1.4) && (power <= 1.5) ) {
    grid <- c(0.52272502143588,
              -2.37970276267028,
              -4.50105382434541,
              2.01559335438892,
              26.7080425402173,
              53.1211365756307,
              56.5864772151065,
              83.6551659699119,
              218.128756577377,
              414.067827117006,
              586.302276594556,
              925.084114296478,
              1629.94052297329,
              2248.04335981282,
              2325.84860904575,
              3056.70827313074,
              5936.34902374365,
              7512.68000065931,
              -494.046409795432,
              -14632.889451627,
              2555.32119471521,
              103793.057384242,
              239732.321730146,
              57045.261742077,
              -1151133.66663602,
              -3657091.74594141,
              0.775913134802564,
              4.5833536253197,
              -6.51027140852265,
              -110.820592889565,
              -302.27161508231,
              -170.893952600532,
              597.013908647114,
              1036.46479866453,
              362.328670850611,
              1143.67350470972,
              5427.58842419152,
              7481.1938127866,
              3916.80746756334,
              10919.4206700089,
              37992.348893648,
              35578.2283061684,
              -48344.4290328334,
              -70796.935934355,
              307989.153472163,
              965852.855353596,
              280146.00264422,
              -4462973.64355295,
              -12247413.3846698,
              -7748551.19364459,
              47938914.6705982,
              196377667.248434,
              -1.15101738661085,
              -3.82183996501946,
              79.2652077448879,
              447.712841790259,
              183.895349671323,
              -4007.2875514913,
              -11554.3202406694,
              -11119.3558179111,
              -3321.49044452746,
              -24694.4788687522,
              -79648.9320167266,
              -73209.1251877769,
              -2021.74771903288,
              -154892.733325762,
              -580785.160917161,
              -291546.650641339,
              1473778.38160541,
              1519601.02525224,
              -7163122.30868759,
              -21848349.4729583,
              -6440961.41109336,
              106830967.379257,
              319671808.797567,
              306212673.299803,
              -944079002.967604,
              -4580578568.46053,
              2.02189579962441,
              -10.0772337050698,
              -303.024532759887,
              -704.351532760201,
              4973.76150820855,
              26142.4723427423,
              31405.9126438149,
              -53806.275525627,
              -149359.118457136,
              16092.9360420975,
              206115.379917112,
              -509055.061292639,
              -1505746.99603438,
              440398.507094803,
              4047882.74286314,
              -3759763.21945456,
              -28362810.7587706,
              -20497981.7777362,
              115049844.871681,
              326041126.501133,
              57527117.4877412,
              -1804448247.28013,
              -5607533850.24561,
              -6742827026.7323,
              11967660233.3111,
              137394861109.439,
              -2.64873193579361,
              68.3499218097747,
              718.463817413271,
              -1571.25771413107,
              -28554.9807730767,
              -66684.6392549349,
              112682.126175196,
              683976.718235879,
              693804.971402528,
              -977072.064073228,
              -991329.886033511,
              5636906.19208592,
              7342984.59011112,
              -20840829.3577062,
              -49193222.8209866,
              42836025.7323296,
              238149829.126778,
              -55474045.455026,
              -1800027569.31187,
              -4174218574.22397,
              -583545105.082449,
              22566852342.6238,
              72170080841.219,
              99488552940.42,
              -95607027136.1,
              2344789183218.79,
              0.265173429404147,
              -224.175784222991,
              -799.783571636822,
              14465.1992899084,
              83422.4433648738,
              -33747.4954822623,
              -1182034.20503755,
              -2562724.00673901,
              1327150.31866679,
              9096932.76169983,
              -8147841.7628163,
              -67997474.423206,
              -84381101.1698542,
              1167786.45920554,
              -290262963.143079,
              -1966559598.21217,
              -4495910335.48392,
              -2778049317.5256,
              9570853857.51408,
              20045770280.4768,
              -39053969152.1498,
              -304481715916.404,
              -861160659351.957,
              -1248601137158.62,
              1149200212746.68,
              152670892302636,
              12.5290983345227,
              499.005760177699,
              -1884.50478980683,
              -52292.5606105439,
              -102007.716028595,
              974234.695542423,
              4265255.61960463,
              -518423.706368479,
              -37847638.3377804,
              -94994579.0086004,
              -141430225.033227,
              -554077045.707268,
              -2621090414.90387,
              -7585723619.93884,
              -14489129113.1742,
              -23451598610.8142,
              -52316888100.9028,
              -154472187133.545,
              -386170628207.441,
              -654732126120.445,
              -458996711714.355,
              1207352849976.42,
              5003368112913.01,
              8162669646404.62,
              10014475965674.9,
              4.30099151054367e+015,
              -49.3184821796504,
              -659.190055418151,
              13760.6854640598,
              117281.576980704,
              -286688.694056081,
              -4514692.65697881,
              -11833557.5599173,
              -13306480.8841683,
              -177569989.290577,
              -1501492416.77122,
              -6595461136.24562,
              -20356400466.8674,
              -54254546277.6263,
              -140855615592.576,
              -350477425656.988,
              -772978728940.34,
              -1445553693725.57,
              -2349378847823.61,
              -3825995343158.07,
              -7818306940077.94,
              -19731827478910.9,
              -48657002077699.6,
              -98504762092397.2,
              -121989778390269,
              454260903660862,
              1.09114573020388e+017,
              123.410260081819,
              -197.403953092516,
              -40101.3176876088,
              -52951.2659905612,
              2144698.8673286,
              1739654.2351533,
              -141996380.772718,
              -1289342757.59316,
              -7156566175.76133,
              -32325427550.7015,
              -126169094996.878,
              -419281859332.1,
              -1174950784066.64,
              -2822039940568,
              -6037668177832.46,
              -12096555834995.5,
              -23612374044250.6,
              -45167080717259.5,
              -82416904456952.1,
              -138051430229770,
              -204193982413479,
              -251800695384389,
              -208815124278779,
              288829117097310,
              9.54395313223954e+015,
              2.24604455857544e+018,
              -210.739560523616,
              5859.20505890453,
              140328.667144469,
              278664.828817554,
              -10534350.9218635,
              -206287796.083656,
              -2616524539.48725,
              -22630598386.6276,
              -138968736141.657,
              -644711350455.876,
              -2387358710732.05,
              -7357036386363.44,
              -19420572655537.7,
              -44732883747785.5,
              -91030332840422.6,
              -165757883253373,
              -276320429460859,
              -441737780262524,
              -728934745124743,
              -1.31361122221926e+015,
              -2.49168560525055e+015,
              -4.31067214964099e+015,
              -4.70964682712407e+015,
              9.73624002713054e+015,
              2.04847011281059e+017,
              4.09466705815479e+019,
              278.504340481087,
              1227.95202954615,
              609240.949141676,
              8398669.53224352,
              -122964349.27965,
              -4120134389.77069,
              -51789735957.0207,
              -413475908382.113,
              -2400754425461.51,
              -10787324023927.9,
              -38888467331936.8,
              -115207923139742,
              -285196579464986,
              -596235892662438,
              -1.05638430488185e+015,
              -1.57370592252881e+015,
              -1.90815996658209e+015,
              -1.68347182504384e+015,
              -470951786697564,
              2.27515760918722e+015,
              8.65004127903261e+015,
              2.98774925717651e+016,
              1.1577555109681e+017,
              4.88806957320561e+017,
              3.5377171188424e+018,
              6.72365027054536e+020,
              1302.48019924681,
              259560.540734449,
              8128828.4793897,
              25490162.9304377,
              -3094217693.64666,
              -72611077275.6044,
              -860362188931.341,
              -6697568428379.51,
              -37943648370441.8,
              -165469263635562,
              -574639366352108,
              -1.61975362474037e+015,
              -3.72208462865178e+015,
              -6.84085667350701e+015,
              -9.24907873336246e+015,
              -5.63230537269553e+015,
              1.49215241150529e+016,
              7.00536184155271e+016,
              1.84186220493898e+017,
              3.92386129591378e+017,
              7.72461936211387e+017,
              1.59366458882732e+018,
              3.85548874292868e+018,
              1.14456878085735e+019,
              5.75359458832251e+019,
              1.02900591814312e+022,
              13286.5118009426,
              2546171.88500447,
              80098454.5079186,
              -575521939.430233,
              -60032299026.9678,
              -1202606432540.75,
              -13452077960416.8,
              -101362184753825,
              -559725106993583,
              -2.37495806056525e+015,
              -7.95202306660676e+015,
              -2.11888175159205e+016,
              -4.40669841594096e+016,
              -6.47957726367261e+016,
              -3.24083222588877e+016,
              1.83291719462516e+017,
              8.40260894845939e+017,
              2.37093537889737e+018,
              5.4317336077623e+018,
              1.10375221140593e+019,
              2.11436475962051e+019,
              4.09291799720987e+019,
              8.65527998606027e+019,
              2.14451316905219e+020,
              8.54313682198285e+020,
              1.49594792908599e+023,
              203942.274027855,
              26468105.5883892,
              462332761.175384,
              -25284736506.3552,
              -1094062620100.06,
              -19258710333321.7,
              -203938699172942,
              -1.48348319911623e+015,
              -7.94195851466556e+015,
              -3.25428940850727e+016,
              -1.03863076023785e+017,
              -2.5561068913892e+017,
              -4.48713805903153e+017,
              -3.40534549726198e+017,
              1.18853947957601e+018,
              6.7824313161356e+018,
              2.1555870304696e+019,
              5.41093701368525e+019,
              1.17698643723088e+020,
              2.32697280763342e+020,
              4.35104831467264e+020,
              8.06863397336643e+020,
              1.57397090527128e+021,
              3.47588746977727e+021,
              1.24037541169807e+022,
              2.082483901283e+024,
              2425722.75874371,
              168829967.462896,
              -6124967796.72205,
              -684794159351.057,
              -20099159890195.9,
              -312088468164664,
              -3.08510744950999e+015,
              -2.12874367646947e+016,
              -1.08016089267394e+017,
              -4.13278775037366e+017,
              -1.18371276048634e+018,
              -2.32798990878473e+018,
              -1.62965158901462e+018,
              9.87177193402913e+018,
              5.68548659646942e+019,
              1.93636952907764e+020,
              5.21964303390126e+020,
              1.21047675802231e+021,
              2.51556116094745e+021,
              4.81371375246211e+021,
              8.68998842506269e+021,
              1.52222181674142e+022,
              2.68233759254245e+022,
              5.09349061051762e+022,
              1.81405407714602e+023,
              2.73961792646156e+025,
              18727182.1067274,
              -3856365344.26257,
              -418983873891.394,
              -17333732195464.3,
              -378965106591800,
              -5.04766648385378e+015,
              -4.45572615875926e+016,
              -2.73473979347584e+017,
              -1.17835431873962e+018,
              -3.3155488196652e+018,
              -3.27497940087677e+018,
              2.43638130489573e+019,
              1.80581815169857e+020,
              7.62727068811908e+020,
              2.48326825663032e+021,
              6.79744241369944e+021,
              1.63174880753305e+022,
              3.52624193788357e+022,
              6.98893205641466e+022,
              1.28897681504313e+023,
              2.23908254454823e+023,
              3.69904872503628e+023,
              5.82169304070961e+023,
              8.85953055589028e+023,
              3.25024671656323e+024,
              3.52294080805527e+026 )
  }
  # 1.3 < p < 1.4 (C)
  
  if ( (power > 1.3) && (power <= 1.4) ) {
    grid <- c(0.897713960567798,
              -0.47040810564112,
              -1.43992915137813,
              -5.23984219701719,
              -15.0352102570745,
              -13.4481420654518,
              87.0680569141198,
              319.051289202574,
              -18.4506632683117,
              -2319.89441748794,
              -3118.4617838911,
              13713.8934655348,
              46257.8044960381,
              -39135.1988030987,
              -498172.85786169,
              -811816.300243115,
              2568821.88452726,
              15957456.1105785,
              30538856.3230172,
              -38265076.8403202,
              -445019432.131326,
              -1535324610.183,
              -2685571226.57687,
              1844994500.64329,
              33799853247.9255,
              148548951511.487,
              -0.0411924170235428,
              -0.280923312875219,
              1.85631524787589,
              35.4548526987637,
              182.842029466341,
              186.535484241998,
              -1885.62536085722,
              -7050.5577782269,
              3914.68802807428,
              72775.307234475,
              81966.2319220996,
              -546556.025881189,
              -1736496.33116182,
              2018406.08511184,
              22431395.8530421,
              38371009.900033,
              -121432141.656335,
              -823775133.379304,
              -1829474775.72627,
              1064685431.79,
              23026018217.409,
              92967101827.417,
              209244156825.792,
              108176144299.254,
              -1494179320052.26,
              -8820798547246.88,
              0.335385258677322,
              4.50345806518915,
              10.608217004451,
              -198.77147800601,
              -1523.07566142071,
              -1500.09456259484,
              23139.2326521245,
              83055.2779182394,
              -113519.71875397,
              -1189391.83133768,
              -880016.414378636,
              11594447.2485996,
              32944023.5477499,
              -57906962.3154716,
              -530418894.202135,
              -898586763.98809,
              3160587598.72435,
              22164670922.5832,
              54438271034.5952,
              -7405444239.37665,
              -615716148031.174,
              -2843896714598.56,
              -7637460249121.44,
              -9706633868213.88,
              27251052973973.8,
              220225910667289,
              -0.792263177969423,
              -21.2045420582639,
              -91.4401569487452,
              1101.7598394385,
              10118.5419251243,
              5726.65994681374,
              -215851.229167288,
              -669893.217282033,
              1905585.55556376,
              13591395.7207961,
              2340171.0149139,
              -169616173.653862,
              -405049039.820392,
              1178374485.90503,
              8651034672.11598,
              13041124217.199,
              -63654520224.3913,
              -427730047807.293,
              -1124157151846.8,
              -274204341784.276,
              11179066056046.7,
              58371933813335.8,
              179156520740047,
              326715205946158,
              -110966263637577,
              -7.30089254118998e+015,
              1.98088010280869,
              88.7205095644789,
              451.25451837431,
              -6238.38424401488,
              -56974.0505836968,
              23643.1529279025,
              1692687.00100628,
              4116879.50794305,
              -20945171.4628138,
              -107415978.096988,
              132557253.129546,
              2070091722.93035,
              3657706483.5209,
              -19429026052.7314,
              -118401385300.158,
              -176359135255.114,
              835686139242.973,
              5703797966574.08,
              14884329900045.1,
              -552681886101.364,
              -188818869070189,
              -1.00878535432539e+015,
              -3.36911029646898e+015,
              -7.50185785257662e+015,
              -5.12493009299165e+015,
              -1.75430627966896e+017,
              -6.01808947590792,
              -358.631844263425,
              -1646.82372824411,
              35116.9756113627,
              280756.645438638,
              -507546.635114422,
              -10020156.6805911,
              -4945705.75502759,
              272336260.80756,
              1104711765.97144,
              -568290834.123902,
              -14088675880.1234,
              -12950574108.0922,
              210192537963.341,
              859065438623.435,
              -747462235105.65,
              -20708041247987.7,
              -104755203973610,
              -298469477349205,
              -363165453972938,
              1.27630293527266e+015,
              1.00317977358032e+016,
              3.80411831558389e+016,
              9.62493626170414e+016,
              1.15786608399486e+017,
              -1.03118547220491e+019,
              22.952229421138,
              1405.69142020857,
              3967.15419654569,
              -181648.014611196,
              -997642.253786536,
              8248812.38538573,
              91248769.6906494,
              269940828.639906,
              457556637.301566,
              7822037651.92945,
              72574078872.7822,
              285799266923.8,
              186491437806.251,
              -3589451347123.35,
              -20096636648266.9,
              -55011152917522.8,
              -75141293052079.6,
              -68199466557531.5,
              -1.03828417550274e+015,
              -1.01420393701492e+016,
              -5.89973268530095e+016,
              -2.54452415456125e+017,
              -8.83842320399224e+017,
              -2.50090967000384e+018,
              -5.75293047799529e+018,
              -3.0288264126223e+020,
              -99.9878158346858,
              -5210.49423696155,
              4761.55706178008,
              1014983.38061518,
              6678222.58484158,
              18138060.5674288,
              523215144.725566,
              8388725586.42611,
              65308743699.6125,
              317410830648.742,
              1134002934158.94,
              3498250050884.25,
              9275186197008.01,
              5168202521004.28,
              -177786455579982,
              -1.541966007589e+015,
              -8.17352077554882e+015,
              -3.26417185677986e+016,
              -1.04838259404326e+017,
              -2.78718323122462e+017,
              -6.27055158391557e+017,
              -1.27711947908717e+018,
              -2.98493427948817e+018,
              -1.08753143842026e+019,
              -6.81489155059389e+019,
              -7.09103794715832e+021,
              445.7016632619,
              18118.8922546317,
              -59575.916393475,
              -1546683.20227499,
              66577416.4833025,
              1581334778.6952,
              18270427567.616,
              159620726230.676,
              1179819333919.81,
              7136701581206.33,
              33141064579502.4,
              108828204724357,
              187282794130913,
              -416027319900596,
              -5.25364928574569e+015,
              -2.85113559834757e+016,
              -1.22522620595716e+017,
              -4.79778248465207e+017,
              -1.80908798010203e+018,
              -6.6584955339892e+018,
              -2.38169794720652e+019,
              -8.23438092498963e+019,
              -2.75142145405642e+020,
              -8.98926942417278e+020,
              -3.30676128118007e+021,
              -1.39506678075182e+023,
              -1872.21347826633,
              -46283.567888389,
              1588637.59435705,
              58147592.3245715,
              1251197774.35208,
              23222717220.2576,
              321130227736.287,
              3181320366971.15,
              23454746025596.1,
              134242541498815,
              608809039493410,
              2.13579959983314e+015,
              4.94319739507647e+015,
              -642441041478663,
              -8.45372610438658e+016,
              -5.78902547364894e+017,
              -2.75526324472539e+018,
              -1.0840919422109e+019,
              -3.76905193553513e+019,
              -1.20991644215191e+020,
              -3.73703358932369e+020,
              -1.16027081029794e+021,
              -3.77617115231973e+021,
              -1.34025753565226e+022,
              -5.99737264585084e+022,
              -2.43233695351887e+024,
              7556.65517982232,
              215106.469973897,
              8930450.66728664,
              660119294.92026,
              21206807863.5722,
              398339343315.102,
              5190034087210.52,
              50568961091221.7,
              380791925478234,
              2.24722396598682e+015,
              1.03784283558222e+016,
              3.62656117204867e+016,
              8.1707753418827e+016,
              -1.42136634601931e+016,
              -1.3434681243398e+018,
              -8.98526482243229e+018,
              -4.27270839102888e+019,
              -1.72138580695355e+020,
              -6.27800609966133e+020,
              -2.15417480180286e+021,
              -7.16024766698911e+021,
              -2.36639068137229e+022,
              -7.98151888085422e+022,
              -2.8511432263888e+023,
              -1.2325305554076e+024,
              -4.09104391632781e+025,
              -22430.9879586296,
              1389052.54923454,
              207349675.587708,
              10195666599.2626,
              296796553336.983,
              5643928776333.53,
              74892895567209.3,
              732432041468099,
              5.49414382342695e+015,
              3.23706890283122e+016,
              1.50367841248151e+017,
              5.32587551175168e+017,
              1.22108658524644e+018,
              -2.40776349600474e+017,
              -2.10569831857027e+019,
              -1.42079098762757e+020,
              -6.7576612032549e+020,
              -2.69444367649145e+021,
              -9.63469651591202e+021,
              -3.22501888704947e+022,
              -1.04839963336164e+023,
              -3.43153835335537e+023,
              -1.17225158559623e+024,
              -4.36196591662948e+024,
              -1.98653397314499e+025,
              -6.88448478890627e+026,
              113277.793293176,
              18517911.460342,
              2450455304.46694,
              133732757996.437,
              3939467103342.68,
              73996025955050,
              973922420426797,
              9.51229943083948e+015,
              7.14628383239514e+016,
              4.211386770595e+017,
              1.94758654233723e+018,
              6.79901662283943e+018,
              1.47957074222858e+019,
              -9.95844999508158e+018,
              -3.10086406715516e+020,
              -2.02215387815306e+021,
              -9.57663684541226e+021,
              -3.84223024829774e+022,
              -1.39073983190476e+023,
              -4.72891356250685e+023,
              -1.5634530908607e+024,
              -5.19829965024042e+024,
              -1.79911666610008e+025,
              -6.76349804228038e+025,
              -3.07798865328979e+026,
              -1.20960393431303e+028,
              262418.415973756,
              250096978.583001,
              31172920464.8184,
              1642871008200.88,
              48188034311242.9,
              904074332252990,
              1.18607968360778e+016,
              1.15420861185206e+017,
              8.6517873709592e+017,
              5.09510325906077e+018,
              2.3530063632105e+019,
              8.13726112964088e+019,
              1.67367130237733e+020,
              -2.20436191917162e+020,
              -4.39432362856901e+021,
              -2.78056375326784e+022,
              -1.30625919423068e+023,
              -5.22563427007622e+023,
              -1.88958355853632e+024,
              -6.42962247956913e+024,
              -2.13314085393945e+025,
              -7.14689344739656e+025,
              -2.50319784231518e+026,
              -9.53595733424881e+026,
              -4.39966873221631e+027,
              -2.25020366975876e+029,
              5301456.42169776,
              2874720507.00564,
              367880928787.512,
              19643816240801.3,
              579437137354404,
              1.09310772186823e+016,
              1.44543936097009e+017,
              1.42253933395633e+018,
              1.08266855728833e+019,
              6.5042335654303e+019,
              3.08534291740846e+020,
              1.11258933788151e+021,
              2.54876590864897e+021,
              -1.16381734289056e+021,
              -5.15203133618314e+022,
              -3.49775856982996e+023,
              -1.70121626254362e+024,
              -6.97428046783589e+024,
              -2.57225393421212e+025,
              -8.89902171064717e+025,
              -2.99264095778419e+026,
              -1.01251337109722e+027,
              -3.56015066420899e+027,
              -1.34824053722105e+028,
              -6.3098335548906e+028,
              -4.46843594631534e+030,
              51969153.4318949,
              37970527243.0619,
              5016911719944.11,
              280591474635532,
              8.66499271311873e+015,
              1.71435134742076e+017,
              2.38968675963669e+018,
              2.4977218053678e+019,
              2.03817113129374e+020,
              1.33075091271939e+021,
              7.03087446839742e+021,
              2.99159383585797e+022,
              9.87178121143328e+022,
              2.15294869859065e+023,
              -1.70137308294605e+022,
              -3.20794310177486e+024,
              -2.12464634368172e+025,
              -9.96697641227305e+025,
              -3.96959854025141e+026,
              -1.44146335908902e+027,
              -5.00266337924097e+027,
              -1.73044335238055e+028,
              -6.20780217337989e+028,
              -2.41346445464175e+029,
              -1.2385425771224e+030,
              -9.58840383704782e+031 )
  }
  
  # 1.2 < p < 1.3 (D)
  if ( (power > 1.2) && (power <= 1.3) ) {
    grid <- c(0.959582735112296,
              -0.196662499936635,
              -0.301224766165899,
              -0.502046137683381,
              -1.5895374772965,
              -7.89601110410117,
              11.3307142810595,
              762.363953933895,
              7332.01631877668,
              29138.4582145098,
              -76339.1153989754,
              -1796970.2346991,
              -11726396.1851003,
              -33450515.2165733,
              98652180.9916008,
              1769595410.03626,
              11618980670.3263,
              49132193375.8717,
              121727250506.162,
              -108084821334.712,
              -3366026267162.29,
              -24251373433568,
              -125561039968382,
              -532052675965940,
              -1.87778088137625e+015,
              -5.20253340562286e+015,
              -0.0117920451301422,
              -0.0865135588789834,
              -0.305001740311124,
              0.805339553985313,
              45.1127229820372,
              543.393289708699,
              2052.1244600063,
              -23479.8112741772,
              -378633.827624817,
              -2160148.0589237,
              -71997.7261576385,
              94765697.6597504,
              790575791.875542,
              3071423100.48738,
              -1263803153.75136,
              -106897577346.909,
              -858303914821.876,
              -4302053589524.5,
              -14378154556128.4,
              -18792439949861.5,
              158250693996112,
              1.65848181037028e+015,
              1.01561192369251e+016,
              4.93870814001288e+016,
              2.03707862081918e+017,
              7.32155754220541e+017,
              0.0276487092578896,
              0.273727106195384,
              3.10493645078303,
              6.74883014953071,
              -791.517414710189,
              -14659.1601086036,
              -97513.5598522334,
              253826.352485024,
              9659750.39303433,
              72426555.4027922,
              108077398.518517,
              -2603530107.50454,
              -27142206424.4266,
              -129010093038.665,
              -115592228710.678,
              3373975177466.36,
              33069166217378.1,
              190970279120807,
              775088407219876,
              1.95531526988697e+015,
              -1.08844022433567e+015,
              -5.03041904941278e+016,
              -3.88737032793942e+017,
              -2.17996177307102e+018,
              -1.03733702308722e+019,
              -4.56687817243336e+019,
              0.00460342184603284,
              -1.4191835457503,
              -61.0896531126379,
              -722.851866770932,
              5719.49208464532,
              231804.146910697,
              2166106.69539304,
              160219.206657996,
              -174724415.264779,
              -1623587498.26411,
              -4008185392.86052,
              51909707843.9719,
              660864914267.222,
              3692989887945.52,
              7286367447296.07,
              -66279529836476.3,
              -829978185197398,
              -5.3925210853509e+015,
              -2.46096213563743e+016,
              -7.65535837790057e+016,
              -7.14762806859001e+016,
              1.1170095179102e+018,
              1.09817546139889e+019,
              6.89213311167972e+019,
              3.45878357931231e+020,
              9.50469739120499e+019,
              -0.135631224023216,
              9.82565018790751,
              752.94430372433,
              13626.2467481671,
              2250.45356087176,
              -2801646.63778468,
              -34431818.5075219,
              -60859775.5058763,
              2414520669.31045,
              26210017396.4224,
              71434279662.7282,
              -950461823661.22,
              -12870038897298.3,
              -75651275138190.2,
              -153087109622533,
              1.57838782399087e+015,
              2.08458508751872e+016,
              1.47352057035772e+017,
              7.65309425661461e+017,
              3.07056896036942e+018,
              8.78357943721663e+018,
              7.63366506587975e+018,
              -1.20231415827737e+020,
              -1.19381796423493e+021,
              -8.4234840567679e+021,
              -2.34505645250979e+023,
              2.86398306366288,
              4.76124711367818,
              -6457.28547283359,
              -164758.27150296,
              -589766.226217846,
              28067458.5191206,
              401139336.264519,
              416911989.44211,
              -40222720768.9054,
              -459374207097.421,
              -1682302721723.89,
              12458076906678.7,
              216651931418543,
              1.51961566257658e+015,
              5.4584519020853e+015,
              -5.19685396766197e+015,
              -2.27022332259378e+017,
              -1.85664800967564e+018,
              -9.86237368369173e+018,
              -3.61931853333836e+019,
              -6.07700089132706e+019,
              3.81186434434807e+020,
              4.88459347203027e+021,
              3.36478940057831e+022,
              1.56888688507016e+023,
              -1.17704898582546e+025,
              -30.8405895817803,
              -700.185732152437,
              47091.7507004153,
              1596708.4034919,
              7159288.38315759,
              -321689179.440809,
              -5431369518.38533,
              -21502235882.8679,
              277109738167.276,
              3653526301164.72,
              5678876796467.48,
              -241265821098598,
              -2.84217859105037e+015,
              -1.41758365590388e+016,
              1.445542591378e+016,
              9.32012909683014e+017,
              9.89863985847271e+018,
              7.14065483941308e+019,
              4.10045064062402e+020,
              1.97534191394671e+021,
              8.07769891709407e+021,
              2.70685067720764e+022,
              6.02927395033448e+022,
              -9.22337544984228e+022,
              -3.1648353300078e+024,
              -5.61603302303069e+026,
              250.323709544488,
              8545.64176293525,
              -338614.039444803,
              -14362564.9062655,
              -92763398.7027217,
              2174220601.06076,
              29854756860.9811,
              -212387816186.495,
              -8028244379539.38,
              -76503721502058.2,
              -140612055233030,
              4.75984421187865e+015,
              6.87179074752808e+016,
              5.46755810778416e+017,
              3.00178419969484e+018,
              1.15075325799871e+019,
              2.56621820643217e+019,
              -9.99967497065107e+018,
              -2.40812092499658e+020,
              1.000269417608e+021,
              2.80929218333346e+022,
              2.76237750806896e+023,
              1.92887184219875e+024,
              1.03524743837727e+025,
              2.54120088308415e+025,
              -1.67374835060249e+028,
              -1788.89737939434,
              -72709.9323496133,
              2484233.26866509,
              105824573.878487,
              46058548.2567914,
              -41008709091.9079,
              -661049991528.956,
              -2869882451384.02,
              34629600071672.9,
              584592196626057,
              3.85506398038373e+015,
              1.85280581626924e+016,
              2.40594389635909e+017,
              4.39015768159735e+018,
              5.68500042585068e+019,
              5.46101970218862e+020,
              4.19859227844407e+021,
              2.71529223079585e+022,
              1.52448364339692e+023,
              7.57284365298019e+023,
              3.350130760701e+024,
              1.30047897758996e+025,
              4.0904969542102e+025,
              5.75746965876781e+025,
              -8.76067257210248e+026,
              -4.18601579421029e+029,
              12278.7581956566,
              496755.477226686,
              -22174016.5110608,
              -982354213.175663,
              -7651083547.2257,
              74418384151.9456,
              -158221447812.226,
              -43635118674338.1,
              -493258537684304,
              3.84848286761028e+015,
              1.75887871994533e+017,
              2.68669123829776e+018,
              2.71502811657687e+019,
              2.09549316384593e+020,
              1.32297232667709e+021,
              7.18370559114638e+021,
              3.53542513917885e+022,
              1.68168650207952e+023,
              8.25241080182299e+023,
              4.30466100873181e+024,
              2.32738116938991e+025,
              1.23477038775956e+026,
              6.04371816144202e+026,
              2.416150272357e+027,
              6.92811011890208e+026,
              -8.30613806320748e+030,
              -79244.2377537492,
              -2985145.62680809,
              144884969.0038,
              3426964337.87191,
              -120514004452.734,
              -4568263757378.07,
              -41894084093912.4,
              490584024082356,
              1.9241818444463e+016,
              3.11732137423038e+017,
              3.73577725478084e+018,
              3.90759600611774e+019,
              3.77869191059759e+020,
              3.36829721948857e+021,
              2.71562326419553e+022,
              1.95921942858408e+023,
              1.26492353680313e+024,
              7.3517449416279e+024,
              3.87469051660304e+025,
              1.86123499753071e+026,
              8.13297579530344e+026,
              3.16817378492111e+027,
              1.01489843927767e+028,
              1.5807400240241e+028,
              -2.021313503304e+029,
              -1.45957971542625e+032,
              521788.722135468,
              9472176.27503808,
              -1834442451.42211,
              -67207743391.7279,
              -811715015878.492,
              -2843970947585.28,
              171797242776594,
              1.07249545636155e+016,
              3.34114929712701e+017,
              6.59909611038221e+018,
              9.29872040038554e+019,
              1.00743611185641e+021,
              8.83945953164831e+021,
              6.53568861391816e+022,
              4.21089667672736e+023,
              2.43823913494771e+024,
              1.30691447936152e+025,
              6.66339503265359e+025,
              3.30039394627076e+026,
              1.60447868150269e+027,
              7.63096240611393e+027,
              3.48620982415311e+028,
              1.46850577313574e+029,
              5.0140288773257e+029,
              -3.94001569880189e+029,
              -2.26766193173512e+033,
              -2825375.34945871,
              -61496352.0334524,
              2493173658.3079,
              -333237181754.782,
              -18996407047060.9,
              -185565155812803,
              9.2705029708914e+015,
              3.99039220138775e+017,
              8.61232487663729e+018,
              1.32051416047768e+020,
              1.60969350838808e+021,
              1.6490151337151e+022,
              1.46466548357974e+023,
              1.14951997029669e+024,
              8.07500862310061e+024,
              5.1281421080201e+025,
              2.96949619046387e+026,
              1.57946761642313e+027,
              7.76011998009655e+027,
              3.52921742493438e+028,
              1.47769012971021e+029,
              5.55807112615323e+029,
              1.72069255611094e+030,
              2.45634325498066e+030,
              -4.33879732862619e+031,
              -3.40784734034525e+034,
              17658315.6693612,
              -1121888367.87085,
              -181270807554.95,
              -7242038986871.72,
              -120082281227313,
              2.31291589803502e+015,
              2.20848277428698e+017,
              7.48543822457755e+018,
              1.61954010343843e+020,
              2.55218480608558e+021,
              3.13237500947909e+022,
              3.13138431695314e+023,
              2.6379745983487e+024,
              1.92552803357783e+025,
              1.24720050890185e+026,
              7.32106996964179e+026,
              3.9678324994807e+027,
              2.01741972525514e+028,
              9.74388679831022e+028,
              4.5059355340519e+029,
              1.99801024333608e+030,
              8.4241281541234e+030,
              3.28635082828553e+031,
              1.06272567840542e+032,
              -2.16825023052698e+032,
              -5.26834661824363e+035,
              -48488127.6322903,
              -9941368958.95411,
              -1701495949545.85,
              -115729090727230,
              -2.46271860575519e+015,
              5.70594335232465e+016,
              4.91770401202199e+018,
              1.54976753508251e+020,
              3.17533901633655e+021,
              4.84563283115304e+022,
              5.87562438473976e+023,
              5.8931457323464e+024,
              5.03145371691109e+025,
              3.73914198714188e+026,
              2.46316680755782e+027,
              1.46055114567583e+028,
              7.89759555011506e+028,
              3.9367543253666e+029,
              1.82414003560443e+030,
              7.89299615934093e+030,
              3.18195160616357e+031,
              1.17447804036354e+032,
              3.74155783420197e+032,
              7.41309337505642e+032,
              -9.98125051400098e+033,
              -9.74447097008623e+036,
              575794042.269378,
              -359696847349.006,
              -46582766175582.5,
              -2.5172826170398e+015,
              -4.85788166089857e+016,
              1.35929935720063e+018,
              1.12759088306295e+020,
              3.61755309466896e+021,
              7.52761720304016e+022,
              1.15671836837518e+024,
              1.39982433553636e+025,
              1.39034166332529e+026,
              1.16802699889137e+027,
              8.49867342754942e+027,
              5.46129812043187e+028,
              3.15144234524267e+029,
              1.65678656217278e+030,
              8.03580116507344e+030,
              3.63383187207197e+031,
              1.54337891859975e+032,
              6.16359372138549e+032,
              2.28426586046228e+033,
              7.47381246679803e+033,
              1.68250380067942e+034,
              -1.71346139855961e+035,
              -2.23719793136107e+038 )
  }
  
  # 1.1 < p < 1.2 (E)
  if ( (power >= 1.1) && (power <= 1.2) ) {
    grid <- c(0.989973425448979,
              -0.111802566436515,
              -0.128138198419178,
              -0.144020752096856,
              -0.156865568330961,
              -0.139657194482135,
              3.536695146237,
              303.956020562169,
              16971.561119878,
              711511.636211499,
              23703218.8754648,
              650925619.121534,
              15122827515.9695,
              302907419807.887,
              5300237420118.51,
              81627637427388.5,
              1.10626268382253e+015,
              1.30088329870981e+016,
              1.26350384665605e+017,
              8.34748022753782e+017,
              -1.48193266830101e+018,
              -1.83141490787536e+020,
              -4.36156639238075e+021,
              -8.28126332215202e+022,
              -1.08728181988638e+024,
              -1.54431621438905e+027,
              -0.00313462642593353,
              -0.0388702138075973,
              -0.0903770763700018,
              -0.178015741979943,
              -0.442864053335794,
              -13.4654296356606,
              -879.397544291817,
              -44326.6781688309,
              -1719945.84474206,
              -52165024.459267,
              -1245203281.74841,
              -23434690273.4911,
              -352258956554.21,
              -4722829656798.21,
              -83507513328313.6,
              -2.39585027904263e+015,
              -7.50927529719736e+016,
              -2.07270273612729e+018,
              -4.97162068932377e+019,
              -1.06766044761913e+021,
              -2.12429857352596e+022,
              -4.05944617085845e+023,
              -7.80448328338183e+024,
              -1.62547597498508e+026,
              -3.41127169961574e+027,
              5.72214427762114e+029,
              0.00507283741509471,
              0.0611283656924362,
              0.118275683282882,
              0.0638053314994744,
              -22.662058007833,
              -2973.06775411472,
              -243880.688361423,
              -13804096.0539272,
              -566905197.509389,
              -17310061654.4572,
              -391792916235.597,
              -6189711302728.5,
              -49123570530885.4,
              601194433055370,
              2.98825837615e+016,
              5.47657960626698e+017,
              3.5405751519937e+018,
              -1.26138966768788e+020,
              -5.9951048018011e+021,
              -1.65556052103626e+023,
              -3.74652939321509e+024,
              -7.7582256629277e+025,
              -1.58711907280473e+027,
              -3.48648576340569e+028,
              -7.77485808458045e+029,
              1.03410922602955e+031,
              0.000415849616352672,
              0.0106967682516999,
              0.158203224877345,
              4.49743806433715,
              -1145.2512510092,
              -257762.508706388,
              -25867475.247582,
              -1637247291.63787,
              -72907884189.7384,
              -2409044555520.72,
              -60655326006420.5,
              -1.17024127884226e+015,
              -1.69507534186056e+016,
              -1.73353934726338e+017,
              -1.19768287172967e+018,
              -1.6925948580369e+019,
              -8.40082861320748e+020,
              -3.19703011126093e+022,
              -9.01680307946117e+023,
              -2.10464678914191e+025,
              -4.37986494585875e+026,
              -8.62023394900792e+027,
              -1.70595448526874e+029,
              -3.65540979450112e+030,
              -7.98948970235937e+031,
              -8.72229109997141e+033,
              -0.000277296401191326,
              0.00785581341291909,
              7.70945326818073,
              1047.85637415002,
              12617.6081671164,
              -11318694.4184715,
              -1521632559.16444,
              -109651752055.073,
              -5365129636121.41,
              -195049584142469,
              -5.55257586362403e+015,
              -1.28875628635706e+017,
              -2.54416168432133e+018,
              -4.52586741722226e+019,
              -7.85882741468167e+020,
              -1.43675527082477e+022,
              -2.81191719043273e+023,
              -5.63847556272856e+024,
              -1.10380771755939e+026,
              -2.06538042850477e+027,
              -3.71079218491912e+028,
              -6.55949719267708e+029,
              -1.19216339970128e+031,
              -2.37651568510983e+032,
              -4.91119938960528e+033,
              -9.19298686871738e+035,
              -4.0453365352258e-005,
              0.808120615148473,
              462.33674259018,
              75430.0290598068,
              4331914.30073851,
              -205449452.816475,
              -53008737555.3822,
              -4488085605111.43,
              -241208396247575,
              -9.50346382864154e+015,
              -2.94900174959011e+017,
              -7.58627631848485e+018,
              -1.69435784683881e+020,
              -3.4411087434676e+021,
              -6.63519114738198e+022,
              -1.2514488828429e+024,
              -2.32881811283864e+025,
              -4.24297462063675e+026,
              -7.48476399731269e+027,
              -1.27258719106238e+029,
              -2.09917093736239e+030,
              -3.42920519941086e+031,
              -5.77122270489562e+032,
              -1.06549974420896e+034,
              -2.07816556556317e+035,
              -4.87536422889045e+037,
              0.00675158161047212,
              30.709935211898,
              17928.1999465246,
              3251282.07045102,
              259150408.193801,
              4065970794.89839,
              -1081487969698.15,
              -120396596572924,
              -7.24129443659267e+015,
              -3.08345729970917e+017,
              -1.02588642970565e+019,
              -2.83199105144724e+020,
              -6.79734935044741e+021,
              -1.47579853708826e+023,
              -2.99031990839464e+024,
              -5.7671579517374e+025,
              -1.06688395883241e+027,
              -1.89309297237781e+028,
              -3.21672741997238e+029,
              -5.24303649293225e+030,
              -8.26753906534395e+031,
              -1.28566458271899e+033,
              -2.0446711266049e+034,
              -3.53884600906098e+035,
              -6.61800844109463e+036,
              -1.75572778206776e+039,
              0.258679354840137,
              840.360258770242,
              507444.100015697,
              99442112.1901258,
              9368084176.4115,
              380375494385.308,
              -9258609923520.86,
              -2.23106090141392e+015,
              -1.58174211439346e+017,
              -7.3607104165464e+018,
              -2.61924232181302e+020,
              -7.67289085556648e+021,
              -1.94419837765708e+023,
              -4.42140613397016e+024,
              -9.26766241853947e+025,
              -1.82037390050501e+027,
              -3.37830777373629e+028,
              -5.94348917058151e+029,
              -9.9356636154327e+030,
              -1.58518493151072e+032,
              -2.43604478698234e+033,
              -3.67029806957723e+034,
              -5.60458388706875e+035,
              -9.22449768201003e+036,
              -1.6904092353277e+038,
              -4.89465564411877e+040,
              6.41800151101403,
              18182.0558161604,
              11380508.3514703,
              2361269590.73996,
              245103026658.95,
              12999718384259.7,
              144508294149080,
              -3.1075671302527e+016,
              -2.81566498473887e+018,
              -1.44690887425048e+020,
              -5.49152466538323e+021,
              -1.692387361606e+023,
              -4.47164792383022e+024,
              -1.05093110630748e+026,
              -2.25304900366832e+027,
              -4.47763141879538e+028,
              -8.32719512434137e+029,
              -1.45736744448634e+031,
              -2.41132989810688e+032,
              -3.79399183745648e+033,
              -5.72975750077454e+034,
              -8.44163520665995e+035,
              -1.2506421655368e+037,
              -1.98161951684487e+038,
              -3.63563587421314e+039,
              -1.14232297834824e+042,
              130.4805948356,
              319626.889155575,
              206944151.13871,
              44870875123.7641,
              4947555746666.1,
              294990039461152,
              6.86125098068619e+015,
              -3.63786637655628e+017,
              -4.49266130472303e+019,
              -2.52923657396585e+021,
              -1.01259704639731e+023,
              -3.24472934374372e+024,
              -8.83324503429516e+025,
              -2.1213437136394e+027,
              -4.60888186519505e+028,
              -9.21064839386007e+029,
              -1.71143638740235e+031,
              -2.97850907929798e+032,
              -4.88496458063687e+033,
              -7.6009042549272e+034,
              -1.13248736207305e+036,
              -1.64007876899463e+037,
              -2.37417054431763e+038,
              -3.65864990028412e+039,
              -6.84510916139291e+040,
              -2.35049544393324e+043,
              963.245715385098,
              4945125.74645911,
              3444934228.87347,
              766212865526.459,
              85513510174417.3,
              5.0967420057197e+015,
              1.1088205871047e+017,
              -7.45318428114027e+018,
              -8.70993014114392e+020,
              -4.89483272136854e+022,
              -1.96863338240187e+024,
              -6.33450769555172e+025,
              -1.72614090113824e+027,
              -4.13046184681991e+028,
              -8.89861663528251e+029,
              -1.7563102367483e+031,
              -3.21437889043872e+032,
              -5.50329625743619e+033,
              -8.87758069741474e+034,
              -1.35900361585279e+036,
              -1.99194876368027e+037,
              -2.83400436317181e+038,
              -4.01653931243524e+039,
              -6.05394644066378e+040,
              -1.16801065759799e+042,
              -4.42572447936007e+044,
              51233.824363294,
              71261780.1129107,
              41401713028.4753,
              8934258413489.44,
              988188857001118,
              5.68035904215963e+016,
              8.47175769846526e+017,
              -1.32231153172396e+020,
              -1.32772648545315e+022,
              -7.30083365623776e+023,
              -2.93724557054933e+025,
              -9.51930409068657e+026,
              -2.6176588474417e+028,
              -6.31768365981176e+029,
              -1.37023407131264e+031,
              -2.71586983800496e+032,
              -4.97904709617111e+033,
              -8.52052408526674e+034,
              -1.37148020128378e+036,
              -2.09209859073969e+037,
              -3.05157815956964e+038,
              -4.31198694357118e+039,
              -6.0516153334399e+040,
              -9.04746158745842e+041,
              -1.83089018411669e+043,
              -7.85792264499214e+045,
              -1497353.00381781,
              723041183.21408,
              798570121442.665,
              179799692875761,
              1.75400882707701e+016,
              5.90805036793554e+017,
              -4.49390336205661e+019,
              -7.12757601895464e+021,
              -4.99886165418825e+023,
              -2.40376807666057e+025,
              -8.92056369816489e+026,
              -2.71249252736027e+028,
              -7.02691152031126e+029,
              -1.59610221811831e+031,
              -3.25079147118034e+032,
              -6.04362056723916e+033,
              -1.04033863298669e+035,
              -1.67721342078464e+036,
              -2.55641247592603e+037,
              -3.71487059097827e+038,
              -5.19156012641781e+039,
              -7.05767849053118e+040,
              -9.54079997480186e+041,
              -1.37682173010483e+043,
              -2.81944608373897e+044,
              -1.35568431120779e+047,
              50890062.0717611,
              14367608874.8046,
              -2441428268667.74,
              -1.1038760824906e+015,
              -1.39748101423647e+017,
              -9.39850732312024e+018,
              -3.91099001869365e+020,
              -1.10065830730325e+022,
              -2.72302043182617e+023,
              -1.10817009473656e+025,
              -6.4123402563569e+026,
              -3.20104135366767e+028,
              -1.26609617693228e+030,
              -4.08925485679342e+031,
              -1.11710339860563e+033,
              -2.65351306408463e+034,
              -5.59921841894524e+035,
              -1.06786709570634e+037,
              -1.86801008380294e+038,
              -3.03739730383572e+039,
              -4.65152896853984e+040,
              -6.8118589248791e+041,
              -9.79544921447061e+042,
              -1.51514184913643e+044,
              -3.61763707621482e+045,
              -2.33635753945665e+048,
              -1632017496.10562,
              110317600316.58,
              330180395951750,
              6.7231021628206e+016,
              3.97828600685863e+018,
              -2.95784859710759e+020,
              -6.98592705045258e+022,
              -6.25659244735647e+024,
              -3.64861751928678e+026,
              -1.5833753360599e+028,
              -5.45880499049355e+029,
              -1.55850222438659e+031,
              -3.79783004095228e+032,
              -8.08967352013915e+033,
              -1.53636362179642e+035,
              -2.64622016202475e+036,
              -4.19608871034953e+037,
              -6.20828443819136e+038,
              -8.67464415143325e+039,
              -1.15730260817286e+041,
              -1.48916975227492e+042,
              -1.86658981818677e+043,
              -2.31349724473074e+044,
              -3.02888770555798e+045,
              -5.99835061959248e+046,
              -4.23733375004766e+049,
              38417214435.0839,
              -726039946748.135,
              -7.31431915372476e+015,
              -1.56745545878788e+018,
              -1.02395774227146e+020,
              4.67055145073008e+021,
              1.3399291038475e+024,
              1.20997806353346e+026,
              6.95795172708743e+027,
              2.95024225282919e+029,
              9.86887837596031e+030,
              2.71349946947451e+032,
              6.30942685965434e+033,
              1.26704450603725e+035,
              2.23282870445361e+036,
              3.49345248902422e+037,
              4.88842278953211e+038,
              6.12750579137622e+039,
              6.82428986341838e+040,
              6.55589778514714e+041,
              4.93603193443648e+042,
              1.70522518557545e+043,
              -3.25119471661806e+044,
              -1.19081432094687e+046,
              -5.38145156298819e+047,
              -9.5337915116242e+050 )
  }
  
  
  # 2 < p < 3 (F)
  
  if ( (power > 2) && (power < 3) ) {
    grid <- c(
      0.999320859444135,
      -0.00339863996520666,
      -0.017405834743935,
      -0.094415885413374,
      -0.474943122220438,
      -2.12217526114901,
      -8.28657088560929,
      -28.1983730239296,
      -83.9866444485264,
      -220.714016257597,
      -516.840853257854,
      -1089.80906886319,
      -2091.04653592782,
      -3687.84899560698,
      -6035.14688398691,
      -9245.16287671588,
      -13364.8056799522,
      -18368.6312623177,
      -24170.5284208701,
      -30652.0675548239,
      -37701.6909476229,
      -45257.9298753625,
      -53351.6680939871,
      -62148.0795745321,
      -71944.5580351856,
      -92584.6307265222,
      0.283341460167403,
      1.4198595014972,
      7.26323206149768,
      39.3470798317175,
      197.737872455163,
      882.775423162057,
      3442.34491552515,
      11680.2981292535,
      34586.2312661717,
      89941.2772399677,
      207071.437085722,
      425804.06636173,
      789124.09455765,
      1329886.46550453,
      2055949.50037776,
      2940437.73707706,
      3922432.53173063,
      4918772.14066394,
      5842643.29870822,
      6622148.2558007,
      7212994.5161919,
      7602717.43947442,
      7807280.48778773,
      7863250.68212219,
      7817746.83340438,
      7725050.76243782,
      -0.0600764380680103,
      -0.392050079124554,
      -1.58937633747626,
      -6.08664376399068,
      -21.2367358176583,
      -63.2532442004164,
      -151.966790509734,
      -262.308768568861,
      -171.730813184702,
      848.887604105088,
      4456.38515028148,
      13551.4810243792,
      32150.4572347722,
      64571.5085421313,
      114045.803910683,
      181267.784390237,
      263555.257000779,
      355057.226722376,
      447959.727920055,
      534210.69005017,
      607115.455356772,
      662364.32312859,
      698095.516091012,
      715453.61053991,
      709404.020526808,
      746938.05963745,
      0.0245245775014973,
      0.155539966613345,
      0.314103723867304,
      -0.974507538399082,
      -14.3303629251534,
      -93.9696323047437,
      -452.198879088332,
      -1754.67665143518,
      -5704.65115146005,
      -15896.7164855003,
      -38614.4509020089,
      -82918.5184841515,
      -159371.089358326,
      -277309.878496798,
      -441509.504984657,
      -649717.382851132,
      -892380.682977137,
      -1154918.2173388,
      -1421750.28004424,
      -1680648.33849184,
      -1926151.40276689,
      -2161369.37925071,
      -2399069.23604879,
      -2658615.89551391,
      -2996584.32616847,
      -3244309.2211475,
      -0.0114761862276732,
      -0.0592619852303128,
      0.149317524117004,
      3.1187422858429,
      23.1527784014401,
      124.86200951222,
      542.858009682394,
      1973.11003584607,
      6116.61892273702,
      16406.537857693,
      38528.1580303481,
      80007.829774633,
      148150.992289522,
      246211.713959957,
      368704.211125567,
      497587.798271825,
      601013.469645754,
      635274.165995275,
      549071.310137826,
      288010.646988658,
      -203166.790795375,
      -982564.60054751,
      -2119680.41798567,
      -3704741.92691329,
      -5931154.48660554,
      -8809844.1184798,
      0.00575921327240797,
      0.0165592415712064,
      -0.310297398420674,
      -3.50279461877254,
      -22.6809711897805,
      -113.076858198932,
      -465.191104437209,
      -1624.40007498999,
      -4906.66394460468,
      -13042.9156711335,
      -31056.3757863278,
      -67485.3530473821,
      -136357.774862024,
      -260555.840068835,
      -476923.168849248,
      -842508.117841961,
      -1440283.41961548,
      -2382568.94608435,
      -3811813.17259831,
      -5900950.25877304,
      -8858145.80320564,
      -12942579.2179964,
      -18500106.3718362,
      -26024539.6191355,
      -36343654.5689473,
      -51760152.8969368,
      -0.00300526348269945,
      0.00332582644494308,
      0.352539219455832,
      3.25105798234768,
      19.0591425935349,
      87.5918992144986,
      333.900689176084,
      1071.0953680497,
      2874.83712283637,
      6250.68546381489,
      9860.08999845749,
      5357.73872252884,
      -33652.9231630469,
      -170144.676973302,
      -528295.164157734,
      -1318787.55146601,
      -2857009.40475253,
      -5568572.63242854,
      -9983746.86426784,
      -16731892.6976634,
      -26555077.5194214,
      -40365482.5555169,
      -59377234.8688406,
      -85360882.1668992,
      -121173410.369109,
      -177831829.550814,
      0.00156864913600111,
      -0.0131184236007616,
      -0.348599078434079,
      -2.77780309250851,
      -14.4361149519266,
      -58.8320138727439,
      -203.569099773387,
      -645.03374745931,
      -2061.89411085477,
      -7067.36215375172,
      -25074.4673862501,
      -85193.394389473,
      -263487.697425328,
      -729232.04752068,
      -1807739.26844015,
      -4046984.43990115,
      -8268652.76656352,
      -15588975.4539116,
      -27412600.48323,
      -45427572.5997714,
      -71652140.9796545,
      -108601600.144071,
      -159667231.52998,
      -229874415.837102,
      -327427911.37296,
      -486267279.031203,
      -0.000751098308492496,
      0.0183607429610277,
      0.32560155770289,
      2.20414845756873,
      9.30362608159624,
      27.198466156795,
      41.1328701889614,
      -147.601617296233,
      -1778.48998961713,
      -10615.1337836503,
      -48121.4524064771,
      -180117.812148503,
      -576102.516202511,
      -1608834.11717014,
      -3989140.08790138,
      -8909582.40404006,
      -18156969.4268232,
      -34163633.7942921,
      -60005088.1200331,
      -99404892.9765117,
      -156858144.123099,
      -238028512.439587,
      -350642272.707905,
      -506310271.081071,
      -724411047.556449,
      -1085756135.90953,
      0.00023268753144256,
      -0.0214573124116083,
      -0.290521160221242,
      -1.52935064098677,
      -3.59013983096331,
      3.92703520664888,
      60.9297355136416,
      133.554643751598,
      -1021.17664549104,
      -11590.6683029842,
      -66685.879961975,
      -281254.304289336,
      -960925.358517155,
      -2788713.9792639,
      -7077864.08967289,
      -16041484.6298937,
      -33004923.3507173,
      -62506118.7968919,
      -110299815.326543,
      -183381109.539535,
      -290240024.52127,
      -441647425.450367,
      -652423811.084398,
      -945103927.138418,
      -1358058723.14185,
      -2050647607.71471,
      0.000141239642755454,
      0.0232583731015363,
      0.239064459228162,
      0.680275724339044,
      -3.14242038852734,
      -42.3056828515275,
      -274.070438394209,
      -1492.11348836561,
      -7656.84591539787,
      -36265.7134048327,
      -152004.565173214,
      -553517.300322874,
      -1753199.12822194,
      -4878532.09154446,
      -12079077.4876383,
      -26959794.0499863,
      -54922340.2612081,
      -103328174.078604,
      -181512150.271794,
      -300836965.808012,
      -475136824.235409,
      -722051548.703595,
      -1066029150.12698,
      -1544638613.28437,
      -2222964632.64711,
      -3372045968.09847,
      -0.000444860575100238,
      -0.0236847905695715,
      -0.161058411288511,
      0.407895924404698,
      10.7548689475349,
      73.0538463066281,
      281.046556675514,
      383.530300181306,
      -3421.68161502129,
      -33452.5094386606,
      -180847.896775903,
      -735888.039250597,
      -2456210.90539502,
      -7010955.43700536,
      -17571067.8749036,
      -39423967.3829222,
      -80439386.6066522,
      -151267929.917988,
      -265326130.629948,
      -438855452.428628,
      -691566344.59827,
      -1048598947.11475,
      -1544960914.68021,
      -2234995984.89575,
      -3214527134.98847,
      -4893933770.89838,
      0.000703732261226084,
      0.0219231227951077,
      0.0425149877267029,
      -1.78013554662686,
      -18.8024427670158,
      -106.878166315772,
      -491.309303563418,
      -2410.2059871961,
      -13085.0171426766,
      -67405.3817964102,
      -299027.541071522,
      -1121165.69781058,
      -3591096.47322529,
      -9998617.39241346,
      -24622492.5304239,
      -54479632.7305074,
      -109838493.819578,
      -204359729.041694,
      -354963747.711344,
      -581822254.35559,
      -909137795.541055,
      -1367637645.50103,
      -2000281379.29399,
      -2874669331.61444,
      -4113040594.89626,
      -6286557382.99022,
      -0.000959178477233191,
      -0.0172549225992488,
      0.119521088004283,
      3.25076886753573,
      24.469396043305,
      77.1613526002622,
      -413.359872941446,
      -8197.32770301892,
      -68568.1716842124,
      -394782.555626955,
      -1746273.13505017,
      -6261284.87873029,
      -18848727.5557636,
      -48909488.5602613,
      -111724871.937695,
      -228668100.175633,
      -425755681.178738,
      -730883280.355038,
      -1170935158.54954,
      -1770335830.80924,
      -2552404036.30653,
      -3544443902.10624,
      -4787694880.17175,
      -6355795903.3852,
      -8395778630.02711,
      -11697767237.6242,
      -0.000133829436888704,
      -0.00402825444711644,
      -0.39467858862294,
      -4.78902212824088,
      -23.5283465754233,
      -495.238264820242,
      -11644.3071017881,
      -143948.075333492,
      -1133949.74994485,
      -6423704.38453159,
      -28106118.540934,
      -99563713.8931762,
      -295413197.299012,
      -753554043.615336,
      -1687633184.19881,
      -3377063524.57688,
      -6129762082.98725,
      -10226991209.9182,
      -15871393620.5447,
      -23160885564.1831,
      -32101906559.4075,
      -42662039575.6347,
      -54852876075.286,
      -68836142772.9544,
      -85067835018.7585,
      -105221455769.378,
      -0.0159347723438438,
      -0.139677220958585,
      -0.501567031055687,
      -1.15800906623169,
      -0.28145406205746,
      -5860.29004848185,
      -152251.522259703,
      -1873731.64284049,
      -14679841.3179456,
      -82892307.4814804,
      -361969279.219145,
      -1280474002.74372,
      -3795015931.20803,
      -9670673302.28459,
      -21636448308.9111,
      -43251166651.4325,
      -78419385616.2198,
      -130679713013.754,
      -202537118746.243,
      -295128647709.38,
      -408391269373.725,
      -541723400932.74,
      -695004000950.63,
      -869822131123.467,
      -1070898637903.45,
      -1311173239757.92)
  }
  
  # 3<=p<4 (G)
  if ( (power >= 3) && (power < 4) ) {
    grid <- c(
      0.807045029092288,
      -0.742817271116966,
      -3.33628159647649,
      -17.1892563957248,
      -83.606924206003,
      -365.071192253248,
      -1400.71893678865,
      -4694.58855569616,
      -13766.4944192918,
      -35518.0066568666,
      -81236.2491747527,
      -166111.053093699,
      -306340.130694687,
      -514013.401873071,
      -791488.479479209,
      -1127820.45893735,
      -1499215.00485969,
      -1873687.95856514,
      -2218215.22065063,
      -2505773.30210073,
      -2720001.08609826,
      -2856885.35011282,
      -2921562.04807504,
      -2942860.2782209,
      -2724727.65113592,
      -6393599.83992389,
      -0.136539772699055,
      -0.3712773738108,
      -1.34060216897536,
      -6.31281195417687,
      -28.7743296376622,
      -120.000150401644,
      -444.643043043044,
      -1450.18256577969,
      -4160.32884893442,
      -10542.0837079396,
      -23750.546686822,
      -47945.825282172,
      -87450.1861981138,
      -145332.756352579,
      -221913.780561476,
      -313883.835136119,
      -414535.3380698,
      -515112.581103503,
      -606768.028120295,
      -682457.909013515,
      -738016.337582899,
      -773197.016326358,
      -785451.127656309,
      -832107.719453664,
      -65576.1860675511,
      -16681644.1359413,
      0.0390190571539757,
      0.217823357411596,
      1.05307617320898,
      5.46628587502473,
      26.7079674195229,
      116.777871575447,
      448.243317291106,
      1502.23344046383,
      4404.21226928853,
      11360.1463323601,
      25977.6589655352,
      53114.6692363956,
      97962.628202229,
      164425.342777476,
      253336.682488931,
      361326.799670235,
      480960.413361712,
      602200.055474944,
      714666.469016104,
      809833.902052876,
      882644.532433897,
      931195.500052511,
      963253.25246436,
      909728.615872677,
      2118955.81118029,
      -28223901.1772564,
      -0.0135344333535702,
      -0.105123408597516,
      -0.597134560351953,
      -3.30138792730478,
      -16.7678509371076,
      -75.2601777869174,
      -294.380322822538,
      -1000.62051658493,
      -2965.26791173967,
      -7710.19873496618,
      -17729.923018914,
      -36365.1109880616,
      -67103.1049888601,
      -112342.914312332,
      -172033.611258013,
      -242814.982348664,
      -318156.29116042,
      -389547.350581509,
      -448309.499169329,
      -487358.607042042,
      -502288.280654725,
      -491896.568009773,
      -454167.367329454,
      -436702.595676197,
      592554.019288957,
      -28175511.1369623,
      0.00516979878298409,
      0.049186369552098,
      0.313361035619778,
      1.82380477698879,
      9.55014857006874,
      43.7632364583752,
      173.85444871189,
      599.000056509861,
      1801.41110174687,
      4775.62836370548,
      11292.324476147,
      24116.159116646,
      47097.8558813758,
      85097.7219846563,
      143732.570485074,
      228950.891134358,
      346463.780768116,
      501132.836280051,
      696519.783935997,
      934868.31070775,
      1217737.35915428,
      1547476.14417311,
      1929941.85114461,
      2371249.26369004,
      3259938.37856758,
      -10481780.2845814,
      -0.00210569527115553,
      -0.0231405024519351,
      -0.161238923976278,
      -0.980926414519322,
      -5.27476113570811,
      -24.596921426023,
      -98.4679251076239,
      -335.875641752254,
      -966.009103975749,
      -2295.01927050472,
      -4285.53743080301,
      -5305.54520561608,
      298.828932844618,
      25527.1544820149,
      93698.0628522327,
      239219.913855273,
      504807.80240039,
      935434.395727312,
      1571191.70516859,
      2442224.38428192,
      3568464.03604572,
      4965686.14710397,
      6658844.98464334,
      8715913.57478651,
      11333531.2576485,
      12135648.0516836,
      0.000904992444600065,
      0.0111168052401889,
      0.0833664125811974,
      0.527134997928565,
      2.89808864704583,
      13.7770488849601,
      57.6536611318659,
      222.31978415136,
      830.973533709237,
      3052.29525500565,
      10683.2000738177,
      34231.1489571498,
      98080.965980282,
      249787.613056241,
      567853.696586215,
      1162730.51664402,
      2167379.38902465,
      3718634.07876582,
      5936186.98691647,
      8908850.23901353,
      12696377.3865465,
      17351472.7681057,
      22966767.528756,
      29776590.0504127,
      38478019.7505486,
      50867749.5850808,
      -0.000409107208596102,
      -0.00550111018040516,
      -0.0439376208280782,
      -0.288505996551621,
      -1.63212078318116,
      -7.75520903054617,
      -27.7954092957704,
      -45.8656633398937,
      259.464351955424,
      2854.74539011967,
      15714.5555070392,
      63230.6123114516,
      204776.271508092,
      558927.74579377,
      1323826.85797549,
      2780007.90439764,
      5266252.25523076,
      9131789.73343168,
      14682560.3322589,
      22146472.0731672,
      31679055.6650164,
      43421701.98937,
      57625070.1381395,
      74910568.9517981,
      97163770.0643605,
      127754334.17767,
      0.000194307002727792,
      0.00281756674493309,
      0.0237136332062326,
      0.159405381495387,
      0.898230867797415,
      4.60345984485858,
      27.6400397896721,
      205.661534317615,
      1448.56631969149,
      8291.67719559354,
      37903.115045638,
      141316.649850398,
      441087.624260472,
      1180252.50521153,
      2763322.5415243,
      5761704.22812121,
      10864619.9821838,
      18782495.180654,
      30138958.1231858,
      45402276.5830579,
      64900090.6565401,
      88943149.521878,
      118085299.880759,
      153669984.92336,
      199638879.192562,
      263053537.257288,
      -9.6813970562699e-005,
      -0.00149818930643216,
      -0.0133055656398246,
      -0.0945115298691158,
      -0.579908914141087,
      -2.47603853062905,
      4.13987484493922,
      197.011666853647,
      1993.36630427103,
      12842.763799835,
      61521.7609859012,
      234369.520358223,
      739791.084258406,
      1992499.07342665,
      4684730.93725107,
      9797208.48563922,
      18517157.6071782,
      32074859.2918227,
      51560378.2551523,
      77807929.4600213,
      111424466.116615,
      153009280.153805,
      203615359.468966,
      265712174.566699,
      346298256.693629,
      458843429.853769,
      5.0510737242787e-005,
      0.000825454278033001,
      0.00754236535471132,
      0.0518083089767653,
      0.274621994126378,
      1.98114382300995,
      29.5576814444073,
      370.652270203056,
      3159.13901640118,
      19407.1216982271,
      91464.1099352158,
      346424.449467708,
      1091742.85694851,
      2941228.74414017,
      6923668.63289443,
      14504109.9212339,
      27468203.8817623,
      47684788.4333904,
      76836997.6005713,
      116252578.587634,
      166950731.275237,
      229981255.740226,
      307139843.974733,
      402454023.216751,
      526912234.233937,
      702818878.857313,
      -2.7418070102492e-005,
      -0.000471538370173029,
      -0.0046165376568393,
      -0.036587699748026,
      -0.259306297425945,
      -0.513545608507554,
      24.5167655002064,
      422.079080450406,
      3864.47708330576,
      24364.9020805635,
      116380.447641915,
      444742.438017592,
      1411324.94709018,
      3824627.71731913,
      9050677.77134831,
      19052334.7104253,
      36248377.3563104,
      63208294.4047604,
      102300026.230609,
      155465995.050065,
      224290351.769315,
      310469434.35741,
      416810360.956061,
      549298562.77482,
      723618453.467978,
      973334765.794742,
      1.5473116821586e-005,
      0.000276604798913003,
      0.00264969418261052,
      0.0174900279365758,
      0.0766314928362542,
      1.35761435509042,
      35.91128539334,
      507.476924676718,
      4502.64745318921,
      28323.9262461144,
      135970.357806147,
      523202.778300938,
      1672239.7082978,
      4563033.2399971,
      10868160.2019504,
      23016985.4182579,
      44040524.526026,
      77210463.32594,
      125613056.798018,
      191876531.064486,
      278262231.602734,
      387270853.524674,
      522928872.391834,
      693428189.404835,
      919453049.166152,
      1247727868.46718,
      -8.8935422652274e-006,
      -0.000166030604502527,
      -0.00180320520193969,
      -0.0158064250790231,
      -0.112690635415925,
      0.605792100613844,
      35.8920609726225,
      549.294618680751,
      5007.87803105007,
      32039.3397103949,
      155789.404879026,
      605600.459856101,
      1951459.22132761,
      5359784.8952449,
      12832094.7037116,
      27287162.7391152,
      52377387.4480437,
      92054947.2665312,
      150060697.345423,
      229605705.993322,
      333504156.464628,
      464936857.308676,
      629016062.372982,
      835873005.440438,
      1110473643.66999,
      1511045682.99279,
      5.2836771853153e-006,
      0.000101503163727336,
      0.000993429510955962,
      0.00768536258973886,
      0.0762004564807366,
      2.5289608455215,
      61.1254255902454,
      833.768195052616,
      7354.03786958471,
      46348.1408509856,
      223057.664768771,
      859057.342688984,
      2741866.57069041,
      7454083.6797878,
      17650308.942397,
      37089910.5528786,
      70295592.6389546,
      121895207.707387,
      195915436.6333,
      295398279.000812,
      422647796.67442,
      580244039.797793,
      772834327.8546,
      1010119774.0892,
      1316484399.07765,
      1741035939.30733,
      -3.064377771159e-006,
      -5.980922030261e-005,
      -0.000636343983276073,
      -0.00078615973385692,
      0.104702137517939,
      4.49267742988801,
      99.9926447089914,
      1292.20378864764,
      10963.7639617,
      66756.6908405897,
      310381.696918354,
      1152600.74257852,
      3536665.37191601,
      9209232.17803997,
      20794869.4119386,
      41458235.608168,
      74105805.4180075,
      120360968.097685,
      179752435.234823,
      249520137.896054,
      325147964.097433,
      401275805.233735,
      472066391.64127,
      529481285.471072,
      561635505.992595,
      419706664.414131)
  }
  
  # 4 <= p < 5 (H)
  if ( (power >= 4) && (power < 5) ) {
    grid <- c(
      0.699925346022531,
      -0.967077543569014,
      -4.01132613691835,
      -20.1399473550733,
      -96.2422051602767,
      -415.365980035549,
      -1580.21518068211,
      -5262.33523148009,
      -15354.0312963513,
      -39454.0262740431,
      -89938.4982463926,
      -183391.698932919,
      -337402.996179853,
      -564968.193375439,
      -868383.335714499,
      -1235426.26432415,
      -1639939.91489061,
      -2046986.48980844,
      -2420647.49346844,
      -2731722.04945314,
      -2962469.78638038,
      -3110166.85555429,
      -3167407.14273507,
      -3347866.5224723,
      220066.29635677,
      -101623416.405294,
      -0.0856827146015816,
      -0.129194389632916,
      -0.275741081827042,
      -0.995185808644477,
      -3.44304072656882,
      -11.1738124760052,
      -32.3161811060365,
      -82.0820212947264,
      -181.186657587385,
      -345.1190334895,
      -560.487576581505,
      -755.308564356227,
      -780.02789704195,
      -420.782646552896,
      549.203902376783,
      2292.75919352824,
      4844.65304486468,
      8078.19818210918,
      11752.1076384033,
      15489.0074612448,
      19295.1080639715,
      20576.2965127781,
      43458.6905027808,
      -272065.299693591,
      7025098.66546232,
      -233107902.179985,
      0.0168969338431571,
      0.0658060047775105,
      0.252941889440879,
      1.1930373881614,
      5.46017018978347,
      22.7950867197324,
      84.5102486368491,
      275.66150811094,
      790.803379055727,
      2003.61163553412,
      4513.25595683332,
      9109.342888956,
      16611.6941580461,
      27601.4176540827,
      42137.2784788037,
      59588.6192203015,
      78680.9563361091,
      97750.8199011234,
      115126.506416694,
      129453.386054272,
      140101.469157844,
      145829.366024648,
      160224.556969174,
      -43787.4381679082,
      5529231.86708544,
      -215530714.274696,
      -0.0040335782993408,
      -0.0234817794259036,
      -0.110587472703672,
      -0.561343113249565,
      -2.70232831818259,
      -11.6848904632598,
      -44.472343214546,
      -148.047891896767,
      -431.69906130895,
      -1108.54940295543,
      -2525.46917566759,
      -5147.26241339181,
      -9467.85961656983,
      -15855.3598573824,
      -24383.2818096163,
      -34724.30268989,
      -46163.6433390437,
      -57734.3395387668,
      -68421.6708128944,
      -77354.6697154979,
      -83913.2169220258,
      -87816.2036212983,
      -87020.6389051429,
      -128645.623351391,
      1605434.61145123,
      -93593764.193113,
      0.00106368144379509,
      0.0078728133185123,
      0.0423052301671125,
      0.226374767743507,
      1.12595759435215,
      4.97840665743307,
      19.2538811192632,
      64.8572506879996,
      190.743003851787,
      492.585067098186,
      1125.23986845133,
      2291.87964977926,
      4195.71676979692,
      6958.1794888724,
      10533.1890357286,
      14663.2366154889,
      18914.8130296304,
      22801.0859452267,
      25968.0905843622,
      28413.4033257206,
      30721.5445983317,
      34434.2563191554,
      41999.1412856565,
      63751.8453401499,
      142903.993721121,
      -14773309.479337,
      -0.0002988416478197,
      -0.0026128755213848,
      -0.0155064767814191,
      -0.0866816879592185,
      -0.442392050536327,
      -1.99055510244982,
      -7.80955850461462,
      -26.7322098215555,
      -80.4281213186107,
      -215.042835786294,
      -517.452034952092,
      -1134.82818218957,
      -2293.86917654107,
      -4310.82691914404,
      -7575.43567664562,
      -12486.4405956532,
      -19320.884708294,
      -28032.1935354864,
      -37976.4690140615,
      -47537.7814603878,
      -53537.2289629646,
      -50176.3090563924,
      -27111.3184862328,
      35023.1543546029,
      97772.386910872,
      1945675.4370755,
      8.8111852916931e-005,
      0.000874298083941498,
      0.00562299485434772,
      0.0327030242649852,
      0.171726872656,
      0.790684116809639,
      3.1130199548177,
      10.1710449772459,
      26.346202329904,
      48.3548936424336,
      31.8771254839422,
      -186.395662514681,
      -1016.84516919046,
      -3258.27728769355,
      -8162.68516490687,
      -17312.2380768225,
      -32239.1706818026,
      -53802.4297241544,
      -81382.1575234866,
      -111877.895287325,
      -138254.603738189,
      -146999.311848418,
      -113518.232821567,
      3726.73086187838,
      273817.624118751,
      757199.553621758,
      -2.7085879774208e-005,
      -0.000297498761276123,
      -0.00204245827241299,
      -0.0121381151934212,
      -0.061733063855247,
      -0.26699858234589,
      -1.1122126687644,
      -5.31569083447484,
      -27.8731682359897,
      -134.874202708192,
      -553.327324835573,
      -1903.59891170908,
      -5566.04749679175,
      -14086.3410351293,
      -31399.2007566238,
      -62614.6499198601,
      -113234.321754671,
      -187884.675320979,
      -288785.120075029,
      -414032.58145669,
      -555329.364059289,
      -694159.031925558,
      -795349.360834586,
      -800767.324268166,
      -634936.677499812,
      -1390026.09670781,
      8.6547130615233e-006,
      0.000103571257343756,
      0.000763282155189074,
      0.00503495689305218,
      0.0339503192962971,
      0.201182583372745,
      0.679200600696824,
      -1.40781922364333,
      -34.7607658472872,
      -251.167735291424,
      -1214.87163835565,
      -4539.92881673297,
      -13946.9496782576,
      -36547.4380325321,
      -83834.1134982485,
      -171726.961316446,
      -319298.861981019,
      -546376.82435059,
      -870690.377723771,
      -1305228.69338149,
      -1855917.38613467,
      -2519159.50584102,
      -3280158.49857154,
      -4124134.63935406,
      -5101301.89955304,
      -8040658.58499969,
      -2.863420992543e-006,
      -3.6705711047603e-005,
      -0.000267736627702962,
      -0.00118987631537025,
      0.00500374982239699,
      0.0906349021665064,
      -0.03964765353084,
      -8.86995948576699,
      -96.4179103903694,
      -622.102159599423,
      -2925.97709670939,
      -10897.5978947285,
      -33690.7023380288,
      -89258.4433448129,
      -207545.201776322,
      -431710.061827054,
      -816252.874623963,
      -1422300.99639821,
      -2311804.80050625,
      -3543091.25767127,
      -5169984.23754927,
      -7246834.8839518,
      -9846927.88218808,
      -13128306.876303,
      -17545564.4897858,
      -25760072.0632545,
      9.862455038939e-007,
      1.3687405837123e-005,
      0.000130756634353942,
      0.00161205842445782,
      0.0260978720622342,
      0.205692870824142,
      -0.217742919739228,
      -19.0670830309275,
      -201.334161636633,
      -1297.1205613785,
      -6144.28137107098,
      -23122.4216546073,
      -72314.544985341,
      -193853.872583445,
      -455974.877979017,
      -958991.306803016,
      -1832262.21329924,
      -3224343.55883332,
      -5290401.34780242,
      -8183538.60025336,
      -12057410.3744032,
      -17089089.6467557,
      -23542978.7525315,
      -31949315.2967005,
      -43603277.3073714,
      -61461357.0483208,
      -3.3898916523125e-007,
      -4.553631043329e-006,
      -6.3146792385749e-006,
      0.000996685778798012,
      0.0283638961892905,
      0.205028928564685,
      -1.5209878776305,
      -40.9193475021547,
      -399.642312169516,
      -2535.0167059544,
      -12007.7742204405,
      -45413.275846636,
      -142947.623599072,
      -385712.615552065,
      -912728.368992133,
      -1929683.59832148,
      -3702976.31241841,
      -6539242.01547221,
      -10759257.0350438,
      -16681448.9101165,
      -24632750.0692673,
      -35008236.3642549,
      -48422810.8521226,
      -66089231.9656645,
      -90790259.3017384,
      -124782634.831774,
      1.3670317399056e-007,
      2.8922795451404e-006,
      8.0311153292009e-005,
      0.00209218748382464,
      0.042288144829039,
      0.224071380381849,
      -3.88595477558196,
      -80.3828018669491,
      -756.261098584473,
      -4760.28166814639,
      -22556.6461988553,
      -85542.2228609317,
      -270094.51932018,
      -730737.212456536,
      -1732553.94919498,
      -3667198.87154186,
      -7040018.87161964,
      -12429346.6875643,
      -20436837.3840823,
      -31660964.6132857,
      -46731285.933596,
      -66447851.8147418,
      -92106261.376575,
      -126229910.642334,
      -174276152.679859,
      -236992787.643919,
      -2.9515402119269e-008,
      1.4230937296803e-006,
      0.000121964247166344,
      0.00353296955768749,
      0.065444187405576,
      0.228124670269768,
      -8.75709483474552,
      -158.883665823477,
      -1458.29696267044,
      -9116.76392174209,
      -43122.5389213505,
      -163475.6881624,
      -516064.788911388,
      -1395522.20150944,
      -3305660.01902088,
      -6987365.36000612,
      -13391031.6773161,
      -23598208.7876355,
      -38732820.725386,
      -59928676.182718,
      -88431093.3806788,
      -125922659.671926,
      -175220715.536157,
      -241690986.941136,
      -336090539.887824,
      -458661737.622016,
      4.2195953578419e-008,
      4.7262596571073e-006,
      0.000268409542694463,
      0.00711771507088045,
      0.117709590133548,
      0.185527243324104,
      -20.5126719620795,
      -341.724384735611,
      -3067.30039938492,
      -18999.4745727416,
      -89434.205638909,
      -338004.516140821,
      -1064613.61087423,
      -2873576.28161018,
      -6796130.72493253,
      -14346593.3785094,
      -27467963.4091078,
      -48381210.7295654,
      -79426953.5719077,
      -123043357.109127,
      -182047880.504625,
      -260412755.530795,
      -364810110.446681,
      -507417982.521783,
      -710389534.723263,
      -981948479.261186,
      3.2858256896933e-008,
      8.6262465133508e-006,
      0.000619619047250753,
      0.0177352317421679,
      0.286522133653506,
      0.930762013518072,
      -34.6477998067267,
      -621.723655825581,
      -5696.25104328862,
      -35724.7343454217,
      -170023.838772441,
      -649865.670726276,
      -2071812.51834388,
      -5666452.90654193,
      -13596350.0025411,
      -29159538.9365487,
      -56804496.6608088,
      -101967267.133786,
      -170892032.312482,
      -270727256.793951,
      -410267304.386682,
      -601763640.206147,
      -864300376.261109,
      -1229059105.40151,
      -1743431919.71737,
      -2419331563.0187)
  }
  
  # 5 <= p < 7 (I)
  if ( (power >= 5) && (power < 7) ) {
    grid <- c(
      0.786434222068313,
      -0.346012079815851,
      -0.111403879165884,
      -0.37684834413327,
      -0.227942819345689,
      -0.748816199446288,
      -0.429657425382675,
      -1.74160208112774,
      -0.380029940723467,
      -4.75431964095465,
      3.04669074354556,
      -18.7254127758988,
      42.2595203965277,
      -127.075610501907,
      597.354383371277,
      -1783.48701130305,
      11845.8231537032,
      -56908.6196171327,
      443014.347945159,
      -3727125.91186479,
      41915265.5772462,
      -613145856.515179,
      12314442848.2845,
      -353565282768.735,
      15056060852441.6,
      -935512014119462,
      -0.0372367474176292,
      -0.0300330290598776,
      0.0194320567856015,
      -0.0286632698352744,
      0.0445003394042309,
      -0.060948400702475,
      0.132495375767898,
      -0.213544390269728,
      0.423844494115997,
      -1.61814802662431,
      -1.11118707905738,
      -26.5227940118242,
      -67.7293316745046,
      -425.458598662828,
      -1001.46373152249,
      -5317.32161731875,
      -6139.36538655163,
      -71044.3658618384,
      194667.687528656,
      -2788059.62676948,
      30424957.0834719,
      -488372280.977764,
      10689386756.2126,
      -340771077349.598,
      16472908240969.7,
      -1.18753639215006e+015,
      0.00272199211678462,
      0.00600457935696839,
      0.000771571810584139,
      0.00448332513732787,
      0.00189163733536341,
      0.00589886871984612,
      0.00239247793636055,
      -0.0286338123481378,
      -0.315494675308072,
      -2.47958187250564,
      -14.5308793302688,
      -72.977265238482,
      -303.738213023504,
      -1105.5766522905,
      -3387.71681429456,
      -9110.18620747135,
      -18243.3848321049,
      -29784.0538288672,
      75227.02534185,
      -48264.5575731943,
      8222854.19107454,
      -103843212.592982,
      2899046709.27306,
      -108457074214.473,
      6391125630551.74,
      -572642997373694,
      -0.000235725585525655,
      -0.000939712997116092,
      -0.000628940476056143,
      -0.000732868263909586,
      -0.00123615938011299,
      0.000322710464360014,
      0.0114320282205539,
      0.137507245882712,
      1.1142446434786,
      8.13970036772955,
      52.4668442088235,
      302.39429825616,
      1551.75362597238,
      7136.34444891635,
      29595.0973156719,
      111934.511057314,
      390212.830424631,
      1272109.66953064,
      3921732.69419177,
      11664373.480828,
      33562473.6192614,
      96275902.7551133,
      386971758.390603,
      -7779194558.35237,
      869116980003.18,
      -120470494442105,
      2.2131289395219e-005,
      0.000143553411002132,
      0.000170793457061453,
      0.000168846275141174,
      0.000879434613056662,
      0.00980484450010389,
      0.12801774812558,
      1.33913284198253,
      11.5110566068305,
      82.8906614784825,
      510.178293855477,
      2728.20581284727,
      12861.2337022356,
      54151.1413496475,
      206133.384228595,
      717759.486483854,
      2312738.60846531,
      6977878.20433018,
      19960158.1406384,
      54872807.7633121,
      147027968.000616,
      392467947.201468,
      1007010930.32944,
      4277401271.27259,
      -32175521793.3006,
      -5944774644892.04,
      -2.1158714077244e-006,
      -2.2245962793571e-005,
      -3.7048515147495e-005,
      2.6733273608329e-005,
      0.00159538838855198,
      0.0313650954321681,
      0.434832189395928,
      4.6653031612727,
      40.2540676573929,
      287.502036146911,
      1739.44171649193,
      9089.28709664162,
      41709.0038241202,
      170553.998156786,
      629707.080488873,
      2124898.97405362,
      6629823.09908051,
      19347946.9909511,
      53449119.9055938,
      141603511.259636,
      365139375.94298,
      931639294.908808,
      2399554030.90057,
      6365379770.8051,
      6629077555.34454,
      1325257347799.35,
      1.8706218898433e-007,
      3.5595867805602e-006,
      1.0461142077657e-005,
      0.000125822395739863,
      0.00343944847191744,
      0.0663161997681272,
      0.932504916856121,
      10.0628718510742,
      86.680128866696,
      614.451310218334,
      3674.06045960846,
      18916.695687826,
      85354.4120703696,
      342720.990999314,
      1241358.93947498,
      4106484.2634202,
      12551665.0365457,
      35852478.1688061,
      96827356.785738,
      250430599.497147,
      629495711.021131,
      1564844953.08431,
      3922964834.69878,
      10060076758.7741,
      28625631416.721,
      59232587864.0314,
      -1.0736419206779e-008,
      -5.5682318067157e-007,
      1.521132088195e-006,
      0.000169816881450556,
      0.00513471228094518,
      0.101321240932149,
      1.43787147451134,
      15.5501270999697,
      133.573873042078,
      940.8573571264,
      5575.65420221155,
      28399.6411897312,
      126608.618970022,
      501858.833477652,
      1793488.15346414,
      5851297.04271422,
      17631543.2949301,
      49624288.1787248,
      131965198.165583,
      335793721.545266,
      829834387.493952,
      2027957061.72163,
      5001954036.02077,
      12663935117.7047,
      32689069667.4124,
      18114334235.1189,
      -1.2346951657342e-009,
      1.2206227284162e-007,
      3.9297895016142e-006,
      0.000197419028183901,
      0.00602851603167602,
      0.120332460066855,
      1.71628779074478,
      18.5634802061232,
      158.929000041028,
      1113.04191204279,
      6547.27067245127,
      33064.4777866304,
      146043.954634495,
      573317.997343912,
      2028724.86709794,
      6553279.09218201,
      19550737.044832,
      54474555.29211,
      143383996.095597,
      361033832.339612,
      882831234.550636,
      2136130396.56957,
      5224797761.40674,
      13104715569.8387,
      32786630095.1341,
      37594565903.4031,
      7.6808927594202e-010,
      3.8884985498993e-009,
      3.1693128416737e-006,
      0.000185283673344334,
      0.00577829446251914,
      0.116358132374561,
      1.6638196322789,
      17.9690172702516,
      153.197902136379,
      1066.5313667415,
      6229.28866525285,
      31215.0838976393,
      136765.982020265,
      532563.383330351,
      1869646.06197352,
      5993641.65402298,
      17752031.0930607,
      49122405.022279,
      128443983.218799,
      321382731.452622,
      781396537.897085,
      1882476272.22521,
      4594300256.1123,
      11490195036.7978,
      27999004292.5526,
      22293146654.9453,
      -2.3198818654552e-010,
      1.7768432128195e-008,
      2.5170398599975e-006,
      0.000147600250383574,
      0.0046726297592754,
      0.0946608920777015,
      1.35378607562751,
      14.5687692163386,
      123.471744448172,
      853.19624444289,
      4941.88962460754,
      24548.7134403604,
      106623.582662544,
      411715.309842279,
      1434104.29957124,
      4564806.01096974,
      13434953.6203005,
      36971805.6020912,
      96215692.5085437,
      239814563.118944,
      581599219.830851,
      1400939948.6891,
      3430310154.31951,
      8605849617.10048,
      20412486855.5244,
      8789001626.70916,
      7.1680465084153e-011,
      8.2425411536541e-009,
      1.643130414875e-006,
      0.000101010537904274,
      0.00324768919113633,
      0.0660426805148569,
      0.941790746358288,
      10.06409355469,
      84.4700029565598,
      577.065123627681,
      3301.27606364533,
      16190.5323305796,
      69436.106953296,
      264898.991604471,
      912482.279496566,
      2875705.99281888,
      8390894.88493373,
      22923038.4921442,
      59301250.3489504,
      147168407.593844,
      356301731.33071,
      860860680.309354,
      2129050981.40544,
      5406333032.3096,
      12477283740.0382,
      -626655436.574466,
      -1.6592791250243e-011,
      4.0243341048155e-009,
      9.144573335573e-007,
      5.9220160549e-005,
      0.00193782694015436,
      0.039445147224364,
      0.55777780092105,
      5.874590816796,
      48.3891511020367,
      323.416138439675,
      1806.1076174669,
      8634.35391245784,
      36073.4833351782,
      134086.227392197,
      450406.616132974,
      1386158.30589757,
      3956483.9824662,
      10592449.579616,
      26908399.3711585,
      65779957.3387473,
      157925265.23842,
      383776594.467057,
      976109761.141651,
      2581703121.52396,
      5773630017.85034,
      -4997652480.24311,
      3.0254155247433e-012,
      6.480434669686e-010,
      3.8358388175966e-007,
      2.805883798994e-005,
      0.000946081848661616,
      0.0191810731660941,
      0.265048170865662,
      2.69054536709911,
      21.1078472872048,
      132.823226564838,
      690.014628342345,
      3029.32274273717,
      11460.2837105301,
      37977.6328286557,
      111736.598497042,
      294820.966835864,
      701306.659651741,
      1501566.94543448,
      2865072.81069675,
      4832518.62805497,
      7836886.63058089,
      19150010.6873593,
      89185710.2082667,
      421423170.310072,
      903955723.061903,
      -4583125093.25688,
      -4.6611488044825e-012,
      -1.238013643329e-009,
      6.0868733097722e-008,
      8.742211953457e-006,
      0.000328731139452505,
      0.00659863314758232,
      0.0844309956857507,
      0.742711442123046,
      4.57966556822816,
      18.4867246619517,
      26.504351580449,
      -278.749024654208,
      -3001.65994744203,
      -18452.141557694,
      -87884.1942798782,
      -354357.650582793,
      -1264077.22326902,
      -4102494.78750277,
      -12352555.6345389,
      -34956221.4649815,
      -93464824.2667665,
      -234738633.43955,
      -541385094.991335,
      -1108582688.25269,
      -2320229645.40791,
      -923285695.803153,
      -3.6747431215211e-012,
      -1.4003818727403e-009,
      3.2110392985838e-008,
      8.0692472963221e-006,
      0.000344137055565625,
      0.0076331695212749,
      0.108809443680536,
      1.09803763922457,
      8.2754406180277,
      47.8674513539372,
      212.505496482315,
      684.095100622171,
      1139.77502759503,
      -3574.08105308545,
      -43796.4476172793,
      -252629.656994052,
      -1124806.21665939,
      -4304377.68043652,
      -14822703.8520193,
      -47098602.6891216,
      -140045686.50728,
      -391792106.570437,
      -1027557566.96914,
      -2504246941.83045,
      -5877038428.29484,
      -2114725574.7047)
  }
  
  # 7 <= p <= 10 (J)
  if ( (power >= 7) && (power <= 10) ) {
    grid <- c(
      0.77536534618895,
      -0.430563489144848,
      0.248933090844648,
      -0.776691071939754,
      1.42435320848802,
      -3.96378215199569,
      10.6082763732752,
      -32.6646389815812,
      107.594392018696,
      -385.875807619758,
      1554.08106263416,
      -6562.13654153062,
      33311.2180271377,
      -175799.39356045,
      1100307.32618603,
      -8076936.20134624,
      64818006.8415247,
      -672272627.675492,
      7951846880.99446,
      -121489926465.507,
      2346719230615.01,
      -60597733349803.8,
      2.17252324837994e+015,
      -1.13969296410995e+017,
      9.19134050303596e+018,
      -1.13643544485344e+021,
      -0.0228287393493532,
      -0.0194681344424241,
      0.0389473545390899,
      -0.084197161744077,
      0.22422473111649,
      -0.619022478530716,
      1.92391840210546,
      -5.54851462125942,
      31.3124540235083,
      32.9263463638783,
      1526.62316779162,
      8457.00783887994,
      82602.3748106012,
      435383.429491403,
      3109338.73107016,
      12375814.1074866,
      91668768.6318141,
      102755151.774348,
      4154874222.20171,
      -39332097778.8331,
      956298556688.97,
      -25984409388026.1,
      1.02539880527099e+015,
      -6.00621364640884e+016,
      5.55452727721829e+018,
      -8.10234314733144e+020,
      0.000990449241892617,
      0.00245658725606492,
      -0.00243739393253416,
      0.00484724481040902,
      -0.00926456300347046,
      0.0281306359896192,
      0.156239298430668,
      4.17858979355244,
      62.7520656235757,
      780.015382056436,
      7975.4456209287,
      68984.4386566011,
      515762.532144664,
      3384991.53977911,
      19858782.3438055,
      105265045.358707,
      514394379.656016,
      2317182902.4062,
      10175938926.1775,
      37656980976.6547,
      252620644930.558,
      -2174477649470.25,
      136934481128099,
      -9.59053217215135e+015,
      1.11722695207039e+018,
      -2.10701957364907e+020,
      -4.739402348122e-005,
      -0.000231852992426791,
      7.5710640928626e-005,
      -0.000134980675522712,
      0.000354522998021474,
      0.0192668478637093,
      0.520510022906938,
      10.7916294399501,
      169.819567472824,
      2141.74949227175,
      22251.8592082094,
      195421.11468649,
      1479886.25999519,
      9836509.12146506,
      58265320.1095556,
      312023644.738493,
      1530894452.40315,
      6980944196.86386,
      29986172284.6235,
      123833899892.611,
      492927274255.622,
      2078899780231.5,
      6850572520142.16,
      -127562576321897,
      5.83782301712251e+016,
      -2.12174087566134e+019,
      2.1161776328178e-006,
      2.0906300316289e-005,
      7.2258740877989e-006,
      -5.5308991127345e-006,
      0.000518091195594951,
      0.0211849717747649,
      0.651125893562394,
      14.0232486435453,
      229.888540942629,
      2996.41818368473,
      32119.513737,
      290589.577507762,
      2266622.81206172,
      15520120.3234281,
      94786645.0963208,
      523940347.205171,
      2658194740.25114,
      12552187877.3901,
      55962360223.2959,
      239035772367.959,
      991951859205.356,
      4061174998315.7,
      15708445458792.5,
      118144223622463,
      -5.07332493452949e+015,
      6.4479142073582e+016,
      -5.2809354056961e-008,
      -1.8733028739588e-006,
      -2.173112075266e-006,
      -1.4999243577998e-006,
      -1.6130806107812e-005,
      0.0062857185977922,
      0.291729587926022,
      7.72340499458176,
      144.463644174106,
      2082.77967708294,
      24311.1198183042,
      237511.034262544,
      1991173.97650715,
      14616262.9410656,
      95580256.036951,
      565523735.27937,
      3071700431.36731,
      15529390455.3399,
      74043372155.9458,
      336948400501.14,
      1477696012185.8,
      6275079442094.89,
      25790455939029,
      96328050938412.8,
      190447249714597,
      1.00992445696896e+017,
      -6.18972021646e-009,
      1.6716212953349e-007,
      1.819435975603e-007,
      -1.4425153664549e-005,
      -0.000610612063278709,
      -0.0167634347715261,
      -0.311332599080969,
      -3.94138325688247,
      -29.7238495417058,
      6.2555652775239,
      4042.6061664409,
      72350.5669137632,
      851581.13806843,
      7901256.5282858,
      61774104.4395478,
      422241976.017542,
      2587726006.34361,
      14498413734.2669,
      75443525203.5096,
      369254749083.404,
      1715435173133.93,
      7593311590379.09,
      31824142191324.9,
      121720773218198,
      396847554120315,
      -4.20744529251827e+015,
      1.6263773828611e-009,
      -1.5358949677752e-008,
      -2.4571305966315e-007,
      -1.8428794010876e-005,
      -0.000937415911513019,
      -0.0306800211195246,
      -0.708070117922537,
      -12.1299638962403,
      -159.03668317182,
      -1618.62979619491,
      -12693.3036944218,
      -72284.2761884275,
      -217166.845678661,
      1044589.95460791,
      23105099.1246563,
      228679026.770271,
      1721827459.30268,
      11024395805.3666,
      62955703260.573,
      329358475506.641,
      1604477935025.79,
      7336606897800.17,
      31379245002770.4,
      121522052593940,
      347347465989858,
      370932937610640,
      -2.6808728840779e-010,
      7.9020135767862e-010,
      -1.5514883197396e-007,
      -1.6276332336879e-005,
      -0.000876562698279347,
      -0.0303398354210813,
      -0.738225701613061,
      -13.3408979363142,
      -185.677983898719,
      -2038.39425331329,
      -17876.6499694802,
      -124905.284696409,
      -670434.749401926,
      -2347535.07095901,
      618376.316323485,
      94377618.194901,
      987830337.872776,
      7299222860.15746,
      45155167569.2116,
      248286662532.449,
      1249349399406.29,
      5836370519195.45,
      25333886840723.4,
      99430753191700.9,
      299121530322441,
      625102006387724,
      4.0255310414341e-011,
      -4.0075628102411e-010,
      -1.0747626938796e-007,
      -1.1032436842021e-005,
      -0.000614076074677421,
      -0.0218656078360402,
      -0.54523287925633,
      -10.0746614750899,
      -143.319283331831,
      -1611.54363682135,
      -14561.2179931157,
      -106231.861503575,
      -615725.868813294,
      -2629632.88870542,
      -5405277.88267556,
      35284319.2277274,
      545097935.998791,
      4473496821.19931,
      29049471850.145,
      164111610940.043,
      839519077635.141,
      3964583399582.53,
      17362504452757.4,
      69068653702018.9,
      218680025952334,
      491443538772321,
      -5.5715190324839e-012,
      -1.6407435736789e-010,
      -5.7888379466468e-008,
      -6.1468686512493e-006,
      -0.000349575517522456,
      -0.0126681684563627,
      -0.320310340444328,
      -5.984181303379,
      -85.8933670684766,
      -973.163367650677,
      -8853.94643946026,
      -65051.8256679697,
      -380371.836571965,
      -1650092.93269853,
      -3630335.63566798,
      19513096.3150577,
      321805563.17688,
      2676104393.70507,
      17469092428.2355,
      98920174459.9976,
      506594361070.731,
      2395294913189.68,
      10527983468557,
      42415349208286.2,
      142117239446790,
      424957076182122,
      9.0943122447293e-013,
      -7.6362405646957e-011,
      -2.713727877671e-008,
      -2.928550895988e-006,
      -0.000169272201914165,
      -0.00620868754605586,
      -0.158241844155348,
      -2.96853021751442,
      -42.6197249239258,
      -480.920925303812,
      -4332.65386909201,
      -31222.9264278609,
      -175486.089577641,
      -686688.578405898,
      -700984.246915129,
      17606789.7093272,
      210224332.113958,
      1618919804.83086,
      10205343393.9837,
      56638977400.7164,
      286412959136.382,
      1344441166878.75,
      5905112905298.04,
      24102001115608.5,
      86044212091524.3,
      330694952809613,
      -7.2206827308685e-014,
      -2.9539633014077e-011,
      -1.1047488592785e-008,
      -1.2213639725909e-006,
      -7.167213125224e-005,
      -0.00265286298445461,
      -0.0678318219234301,
      -1.2685871768583,
      -18.0238255068467,
      -199.279184183516,
      -1731.66032795093,
      -11672.9085946704,
      -56607.411359395,
      -125551.570376765,
      1014540.31262727,
      16556288.6301798,
      145428947.733874,
      1003290094.11754,
      5972683400.78656,
      32007315379.0508,
      158215899080.749,
      732613701543.087,
      3206486334222.12,
      13279647471500.9,
      50798414623688.9,
      237038150364079,
      4.3609819166476e-014,
      -8.2019208583295e-012,
      -3.8681869833068e-009,
      -4.4736540762345e-007,
      -2.681092767123e-005,
      -0.00100096482939058,
      -0.02552849391404,
      -0.470191072698517,
      -6.46906795194434,
      -67.4507414346289,
      -524.848782738057,
      -2753.53939768566,
      -4110.14385245369,
      104460.727972183,
      1538900.68620935,
      14024098.498795,
      101919977.04296,
      638471700.403706,
      3583252863.16384,
      18458609364.3726,
      88808087125.7936,
      404428978510.69,
      1762294434820.44,
      7420587100023.73,
      30447350490820.6,
      161257587722473,
      1.102927814593e-014,
      -2.9833499013762e-013,
      -1.0237201386522e-009,
      -1.357500505407e-007,
      -8.5194681661709e-006,
      -0.000321263444569207,
      -0.00803242863440226,
      -0.139760568118627,
      -1.70845294566464,
      -13.7392013834384,
      -43.3312232463909,
      672.210515304308,
      14566.0279110099,
      170927.78078521,
      1525821.93755474,
      11366312.5854366,
      73868184.8847336,
      430610301.867199,
      2297726058.33628,
      11406906364.7078,
      53429633348.726,
      239204245061.543,
      1037569229691.63,
      4440670923491.02,
      19385786140661,
      108390210481645,
      2.4429367439691e-015,
      1.4097646088969e-012,
      -8.1610167897946e-011,
      -2.2790490204496e-008,
      -1.5951868755541e-006,
      -5.4515698129209e-005,
      -0.000902748062402921,
      0.000476641583177264,
      0.403718060688101,
      11.3250139182509,
      195.75060974824,
      2526.10978407101,
      26230.959211931,
      228725.678544486,
      1724448.556825,
      11494565.1027429,
      68980892.2243655,
      378540210.929366,
      1926018147.60076,
      9203415540.67664,
      41820439362.7481,
      183069601908.131,
      783979860688.641,
      3359830680474.02,
      15066614162007.1,
      80531744425420.1)
  }
  
  grid
}






#############################################################################
tweedie.dev <- function(y, mu, power)
{
  # 
  # Peter K Dunn 
  # 29 Oct 2000 
  # 
  
  p <- power
  if(p == 1) {
    dev <- array( dim = length(y) )
    mu <- array( dim = length(y), mu )
    dev[y != 0] <- y[y != 0] * log(y[y != 0] / mu[y != 0]) - (y[y != 0] - mu[y != 0])
    dev[y == 0] <- mu[y == 0]
  } else{
    if(p == 2) {
      dev <- log(mu / y) + (y / mu) - 1
    } else{
      if (p == 0) {
        dev <- (y - mu) ^ 2
        dev <- dev / 2
      } else{
        dev <- (y ^ (2 - p)) / ((1 - p) * (2 - p)) - 
          (y * (mu ^ (1 - p)))/(1 - p) + 
          (mu ^ (2 - p)) / (2 - p)
      }
    }
  }
  dev * 2
}


#############################################################################
dtweedie.dldphi <- function(phi, mu, power, y ){
  # Calculates the log-likelihood
  # function, wrt phi, for p>2.  In particular, it returns
  #    sum{ d(log f)/d(phi) } = d( log-likelihood )/d(phi).
  # The mle of phi can be found, therefore, by setting this to zero.
  # y  is generally a vector of observed values.
  
  #
  # Peter Dunn
  # 31 Jan 2001
  #
  
  if ( (power != 2 ) & ( power != 1 ) ) {
    
    k <- phi ^ (1 / (power - 2))
    # cat("k=",k,"\n")
    # cat("phi=",phi,"\n")
    # cat("power=",power,"\n")
    if ( k < 1  & k > 0 ) {
      # Use the transform f(y; mu, phi) = c f(c*y; c*mu, c^(2-p)*phi)
      # and differentiate with c=phi^(1/(p-2)):
      #    d log f / d phi = c^(2-p) * {df(cy; c*mu, 1)/dphi} / f(cy; c*mu, 1)
      f <- dtweedie( y = k * y, 
                     power = power, 
                     mu = k * mu, 
                     phi = 1 )
      d <- dtweedie.dlogfdphi( y = k * y, 
                               power = power, 
                               mu = k * mu, 
                               phi = 1 )
      # Note:  We need dlogf/dphi = dlogf.dphi * f
      top <- d * f
      d <- -2* sum( top / f * k ^ (2 - power) )
      
    } else{
      # Compute directly
      d <- -2 * sum( dtweedie.dlogfdphi(y = y, 
                                        power = power, 
                                        mu = mu, 
                                        phi = phi) )
    }
  } else{
    # Cases p == 1 and  p == 2 
    d <- -2 * sum( dtweedie.dlogfdphi(y = y, 
                                      power = power, 
                                      mu = mu, 
                                      phi = phi) )
  }
  d
}



#############################################################################
dtweedie.dlogfdphi <- function(y, mu, phi, power)
{
  #
  # Calculates d(log f)/d(phi) for the Tweedie
  # densities.
  # It is used, for example, in mle fitting of phi.  We would then
  # sum over  y  and set this function to 0.
  #
  #
  # Peter Dunn
  # 31 Jan 2001
  #
  
  p <- power
  a <- (2 - p) / (1 - p)
  if(length(phi) == 1) {
    phi <- array(dim = length(y), phi)
  }
  if(length(mu) == 1) {
    mu <- array(dim = length(y), mu)
  }
  A <- (y * mu ^ (1 - p)) / (phi ^ 2 * (p - 1))
  B <- mu^(2 - p)/(phi ^ 2 * (2 - p))
  
  if(power > 2) {
    f <- array(dim = c(length(y)))
    # Here, we evaluate logv and everything as normal.
    # If logv has infinite values then we resort to other tactics.
    
    kv <- dtweedie.kv.bigp(power = power, 
                           phi = phi, 
                           y = y)$kv
    dv.dphi <- (kv * (a - 1)) / phi
    out.logv <- dtweedie.logv.bigp(power = power, 
                                   phi = phi, 
                                   y = y)
    
    # Now see if this causes problems.
    logv <- out.logv$logv
    
    # Now detect problem computing  logv  and remedy them
    probs <- (is.infinite(logv)) | (is.nan(logv)) | (y < 1)
    
    if(any(probs)) {
      
      # OK then:  Troubles computing log(V).
      # best we can do is use definition I think.
      delta <- 1.0e-5
      a1 <- dtweedie(power = power, 
                     phi = phi[probs], 
                     mu = mu[probs], 
                     y = y[probs])
      a2 <- dtweedie(power = power, 
                     phi = phi[probs] + delta, 
                     mu = mu[probs], 
                     y = y[probs])
      f[probs] <- (log(a2) - log(a1) ) / delta
      
    }
    
    f[!probs] <- A[!probs] + B[!probs] + dv.dphi[ !probs] / exp(logv[!probs])
  }
  #  END p>2
  
  if(power == 2) {
    f <-   -log(y)  + ( y / mu ) + digamma(1 / phi) - 1 + log( mu * phi )
    f <- f / (phi ^ 2)
  }
  
  if(power == 1) {
    
    f <- mu - y - y * log(mu / phi) + y * digamma(1 + (y / phi))
    f <- f / (phi ^ 2)
  }
  
  if((power > 1) && (power < 2)) {
    # We need to treat y==0 separately.
    # We fill  f  with the Y=0 case
    f <- array(  dim = length(y), 
                 mu ^ (2 - power) / ( phi ^ 2 * (2 - power) )  )
    jw <- dtweedie.jw.smallp(power = power, 
                             phi = phi[y > 0], 
                             y = y[y > 0])$jw
    dw.dphi <- (jw * (a - 1)) / phi[y > 0]
    logw <- dtweedie.logw.smallp(power = power, 
                                 phi = phi[y > 0], 
                                 y = y[y > 0])$logw
    f[y>0] <- A[y > 0] + B[y > 0] + dw.dphi / exp(logw)
  }
  f
}


#############################################################################
dtweedie.interp <- function(grid, nx, np, xix.lo, xix.hi,
                            p.lo, p.hi, power, xix) {
  # Does the interpolation calculation
  
  # Peter Dunn
  # 17 April 2001
  
  # Set-up
  # xi-dimension
  jt <- seq(0, nx, 1)
  ts <- cos((2 * jt + 1)/(2 * nx + 2) * pi)
  ts <- ((xix.hi + xix.lo) + ts * (xix.hi - xix.lo)) / 2    #
  
  # p-dimension
  jp <- seq(0, np, 1)
  ps <- cos((2 * jp + 1) / (2 * np + 2) * pi)
  ps <- ((p.hi + p.lo) + ps * (p.hi - p.lo))/2    #
  
  rho <- array(dim = nx + 1)
  dd1 <- array(dim = nx + 1)    #
  
  # interpolate to find divided differences
  for(i in 1:(nx + 1)) {
    dd1[i] <- grid[i, np + 1]
    for(k in seq(np, 1, -1)) {
      dd1[i] <- dd1[i] * (power - ps[k]) + grid[i, k]
    }
  }
  
  # Now use divided difference to interpolate rho or 1/rho
  rho <- dd1[nx + 1]
  for(k in seq(nx, 1, -1)) {
    rho <- rho * (xix - ts[k]) + dd1[k]
  }
  
  if ( power >= 3) {
    rho <- 1 / rho
  }
  rho
  
}



#############################################################################
dtweedie.inversion <- function(y, power, mu, phi, exact=TRUE, method=3){ 
  # 
  # Peter K Dunn 
  # 06 Aug 2002
  #
  
  # Error checks
  if ( power < 1) stop("power must be greater than 1.")
  if ( any(phi <= 0)) stop("phi must be positive.")
  #if ( power>2) if ( any(y <= 0) ) stop("y must be a positive vector.")
  #if ( (power>1 & (power<2) ) if ( any(y < 0) ) stop("y must be a non-negative vector.")
  if ( any(mu <= 0) ) stop("mu must be positive.")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) {
      stop("mu must be scalar, or the same length as y.")
    }
  } else {
    mu <- array( dim = length(y), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  save.method <- method
  if ( !is.null(method)){
    if ( length(method) > 1 ) {
      method <- save.method <- method[1]
    }
    if ( !(method %in% c(1, 2, 3)) ) stop("method must be 1, 2 or 3 (or left empty).")
  }
  
  y.len <- length(y)
  density <- y
  its <- y
  verbose <- FALSE
  
  if ( is.null(method)){
    method <- array( dim = length(y))
  } else {
    method <- array( method, dim = length(y))
  }
  
  for (i in (1:y.len)) {
    
    # There are three approaches, each a product of a simple bit
    # and a complicated bit computed in FORTRAN
    # We choose Method 3 if no other is requested.
    #
    # The methods are documented in Dunn and Smyth (2008):
    # Method 1: Evaluate a: compute a(y, phi) = f(y; 1, phi)
    # Method 2: Rescale the mean to 1
    # Method 3: Rescale y to 1 and evaluate b.
    
    if ( y[i] <= 0 ) {
      if ( (power > 1) & (power < 2) ) {
        if ( y[i] == 0 ) {
          density[i] <- exp( -mu[i] ^ (2 - power) / ( phi[i] * (2 - power) ) )
        } else {
          density[i] <- 0
        }
      } else {
        density[i] <- 0
      }
    } else {
      # Here, y > 0
      m2 <- 1 / mu[i]
      
      theta <- ( mu[i] ^ (1 - power) - 1 ) / ( 1 - power )
      if ( ( abs(power - 2 ) ) < 1.0e-07 ){
        kappa <- log(mu[i]) + (2 - power) * ( log(mu[i]) ^ 2 ) / 2
      } else {
        kappa <- ( mu[i] ^ (2 - power) - 1 ) / (2 - power)
      }
      m1 <- exp( (y[i]*theta - kappa ) / phi[i] )
      
      dev <- tweedie.dev(y = y[i], 
                         mu = mu[i],
                         power = power )
      m3 <- exp( -dev/(2 * phi[i]) ) / y[i]
      
      min.method <- min( m1, m2, m3 )
      
      # Now if no method requested, find the notional "optimal"
      if ( is.null(method[i]) ) {
        if ( min.method == m1 ){
          use.method <- 1
        } else {
          if ( min.method == m2 ) {
            use.method <- 2
          } else {
            use.method <- 3
          }
        }
      } else {
        use.method <- method[i]
      }
      
      # Now use the method
      # NOTE: FOR ALL  METHODS, WE HAVE mu=1
      
      if ( use.method == 2 ) {
        tmp <- .Fortran( "twpdf",
                         as.double(power),
                         as.double(phi[i] / (mu[i] ^ (2 - power)) ), # phi
                         as.double(y[i] / mu[i]), # y
                         as.double(1), # mu
                         as.integer( exact ), #exact as an integer
                         as.integer( verbose ), #verbose as an integer
                           # NOW WHAT FOLLOWS ARE THE OUTPUTS:
                         as.double(0), # funvalue
                         as.integer(0), # exitstatus
                         as.double(0), # relerr
                         as.integer(0)) # its
#                         PACKAGE="tweedie")
        
        den <- tmp[[7]]
        density[i] <- den * m2
        
      } else {
        if ( use.method == 1 ) {
          tmp <- .Fortran( "twpdf",
                           as.double(power),
                           as.double(phi[i]), # phi
                           as.double(y[i]), # y
                           as.double(1), # mu
                           as.integer( exact ), #exact as an integer
                           as.integer( verbose ), #verbose as an integer
                             # NOW WHAT FOLLOWS ARE THE OUTPUTS:
                           as.double(0), # funvalue
                           as.integer(0), # exitstatus
                           as.double(0), # relerr
                           as.integer(0)) # its
#                           PACKAGE="tweedie")
          
          den <- tmp[[7]]
          density[i] <- den * m1
          
        } else { # use.method == 3
          tmp <- .Fortran( "twpdf",
                           as.double(power),
                           as.double(phi[i] / (y[i] ^ (2 - power))), # phi
                           as.double(1), # y
                           as.double(1), # mu
                           as.integer( exact ), #exact as an integer
                           as.integer( verbose ), #verbose as an integer
                              # NOW WHAT FOLLOWS ARE THE OUTPUTS:
                           as.double(0), # funvalue
                           as.integer(0), # exitstatus
                           as.double(0), # relerr
                           as.integer(0)) # its
#                           PACKAGE="tweedie")
          
          den <- tmp[[7]]
          density[i] <- den * m3
        }
      }
    }
  }
  
  density
  
}


#############################################################################
dtweedie.jw.smallp <- function(y, phi, power ){ 
  #
  # Peter K Dunn
  # 18 Jun 2002
  #
  
  #
  # Error traps
  #
  
  if ( power < 1) stop("power must be between 1 and 2.")
  if ( power > 2) stop("power must be between 1 and 2.")
  if ( any(phi <= 0) ) stop("phi must be strictly positive.")
  if ( any(y <= 0) ) stop("y must be a strictly positive vector.")
  
  #
  # Set up
  #
  p <- power
  a <- ( 2 - p ) / ( 1 - p )         # Note that a<0 for 1<p<2
  
  a1 <- 1 - a
  r <- -a * log(y) + a * log(p - 1) - a1 * log(phi) - log(2 - p)        # All terms to power j
  
  drop <- 37             # Accuracy of terms: exp(-37)
  #
  # Find limits of summation using Stirling's approximation
  # to approximate gamma terms, and find max as j.max
  #
  logz <- max(r)             # To find largest  j  needed
  j.max <- max( y^( 2 - p ) / ( phi * (2 - p) ) )
  j <- max( 1, j.max )
  
  c <- logz + a1 + a * log(-a)
  wmax <- a1 * j.max
  estlogw <- wmax
  
  # First, the upper limit of j
  
  while(estlogw > (wmax - drop) ){
    j <- j + 2
    estlogw <- j * (c - a1 * log(j))
  }
  
  hi.j <- ceiling(j)
  
  # Now the lower limit of j
  logz <- min(r) 
  j.max <- min( y ^ ( 2 - power ) / ( phi * (2 - power) ) )
  
  j <- max( 1, j.max)
  wmax <- a1 * j.max 
  estlogw <- wmax 
  
  while ( ( estlogw > (wmax - drop) ) && ( j >= 2) ) {
    j <- max(1, j - 2)
    estlogw <- j*(c - a1 * log(j))
  }
  
  lo.j <- max(1, floor(j))
  # 
  # Now sum the series between established limits.
  # We ensure it works for vector y.
  
  j <- seq( lo.j, hi.j)                 # sequence of j terms needed
  o <- matrix( 1, nrow = length(y))     # matrix of ones
  
  g <- matrix(lgamma( j + 1 ) + lgamma( -a * j ),
              nrow = 1, 
              ncol = hi.j - lo.j + 1) 
  logj <- matrix(log(j),
                 nrow = 1, 
                 ncol = hi.j - lo.j + 1) 
  og <- o %*% g                             # matrix of gamma terms
  ologj <- o %*% logj 
  A <- outer(r, j) - og + ologj        # the series, almost ready to sum
  m <- apply(A, 1, max)                    # avoid overflow; find maximum values
  we <- exp( A - m )                    # evaluate terms, less max.
  sum.we <- apply( we, 1, sum)            # sum terms
  jw <- sum.we * exp( m )                # now restore max.
  # Since derivs may be negative, can't use log-scale
  
  list(lo = lo.j, 
       hi = hi.j, 
       jw = jw, 
       j.max = j.max )
  
}



#############################################################################
dtweedie.kv.bigp <- function(y, phi, power){ 
  # 
  # Peter K Dunn 
  # 18 Jun 2002
  # 
  
  #
  # Error traps
  #
  
  if ( power < 2) stop("power must be greater than 2.")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y <= 0) ) stop("y must be a strictly positive vector.")
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  
  
  #
  # Set up
  #
  
  p <- power
  a <- ( 2 - p ) / ( 1 - p )
  a1 <- 1 - a
  
  r <- -a1 * log(phi) - log(p - 2) - a * log(y) + a * log(p - 1)
  drop <- 37
  
  #
  # Now we find limits of summation, using Stirling's approximation
  # to approximate the gamma terms, and find max as k.max.
  #
  
  logz <- max(r)
  k.max <- max( y ^ (2 - p) / ( phi * (p - 2) ) )
  k <- max( 1, k.max )
  
  c <- logz + a1 + a * log(a)
  vmax <- k.max * a1
  estlogv <- vmax
  #
  # Now we search either side for when we can
  # ignore terms as they are negligible.
  # Remember we are using the envelope as the tests,
  # not individual terms in the series.
  
  # First:  the upper limit of k
  
  while ( estlogv > (vmax - drop) ) {
    k <- k + 2
    estlogv <- k * ( c - a1 * log(k) )
  }
  
  hi.k <- ceiling(k)
  #
  # Now the lower limit of k
  #
  logz <- min(r)
  k.max <- min( y ^ (2 - p) / ( phi * (p - 2) ) ) 
  k <- max( 1, k.max )
  c <- logz + a1 + a * log(a)
  vmax <- k.max * a1
  estlogv <- vmax
  
  while ( (estlogv > (vmax - drop) ) && ( k >= 2) ) {
    k <- max(1, k - 2)
    estlogv <- k * ( c - a1 * log(k) )
  }
  
  lo.k <- max(1, floor(k) )
  
  #
  # Now sum the series  between established limits.
  # We ensure it works for vector y.
  #
  k <- seq(lo.k, hi.k) 
  o <- matrix( 1, nrow = length(y)) 
  
  g <- matrix( lgamma( 1 + a * k) - lgamma(1 + k),
               nrow = 1, 
               ncol  =length(k) )
  logk <- matrix( log(k),
                  nrow = 1, 
                  ncol = length(k) )
  
  og <- o %*% g
  ologk <- o %*% logk
  
  A <- outer(r, k) + og + ologk
  
  C <- matrix( sin( -a * pi * k ) * (-1) ^ k,
               nrow = 1, 
               ncol = length(k) )
  C <- o %*% C
  
  m <- apply(A, 1, max)
  ve <- exp(A - m)
  sum.ve <- apply( ve*C, 1, sum )
  kv <- sum.ve * exp( m )
  # Since derivs may be negative, can't use log-scale
  
  list(lo = lo.k, 
       hi = hi.k, 
       kv = kv, 
       k.max = k.max )
  
}


#############################################################################
qtweedie <- function(p, xi = NULL, mu, phi, power = NULL){
  
  # Sort out the xi/power notation
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  # Error checks
  if ( any(power < 1) )  stop( paste(index.par.long, "must be greater than 1.\n") )
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(p < 0) ) stop("p must be between zero and one.\n")
  if ( any(p > 1) ) stop("p must be between zero and one.\n")
  if ( any(mu <= 0) ) stop("mu must be positive.\n")
  if ( length(mu) > 1) {
    if ( length(mu) != length(p) ) stop("mu must be scalar, or the same length as p.\n")
  } else {
    mu <- array( dim = length(p), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(p) ) stop("phi must be scalar, or the same length as p.\n")
  } else {
    phi <- array( dim = length(p), phi )
    # A vector of all phi's
  }
  # Now, if mu is a vector, make power correct length
  if ( length(power) == 1 ) {
    if ( length(mu) > 1 ) {
      power <- rep( power, length(mu) ) 
    }
  }
  
  #
  # End error checks
  #
  
  len <- length(p) 
  
  # Some monkeying around to explicitly account for the cases p=1 and p=0
  ans <- ans2 <- rep( NA, length = len )
  if ( any(p == 1) ) ans2[p == 1] <- Inf
  if ( any(p == 0) ) ans2[p == 0] <- 0
  
  ans     <-  ans[ ( (p > 0) & (p < 1) ) ]
  mu.vec  <-  mu[ ( (p > 0) & (p < 1) ) ]
  phi.vec <-  phi[ ( (p > 0) & (p < 1) ) ]
  p.vec   <- p[ ( (p > 0) & (p < 1) ) ]
  
  for (i in (1 : length(ans)) ) {
    
    mu.1 <- mu.vec[i]
    phi.1 <- phi.vec[i]
    p.1 <- p.vec[i]  # This is the  qtweedie()  input p (a probability)
    pwr <- power[i]  # This is the Tweedie power, xi
    
    prob <- p.1 # Rename p to avoid confusion with  pwr: This is the  qtweedie()  input p (a probability)
    
    if ( pwr < 2 ) {
      qp <- qpois(prob, 
                  lambda = mu.1 / phi.1)
      if ( pwr == 1 ) ans[i] <- qp   
    }
    
    qg <- qgamma(prob,  
                 rate = 1 / (phi.1 * mu.1), 
                 shape = 1 / phi.1 )
    if ( pwr == 2 ) ans[i] <- qg
    
    # Starting values
    # - for 1<pwr<2, linearly interpolate between Poisson and gamma
    if ( (pwr > 1) & ( pwr < 2) ) {
      start <- (qg - qp) * pwr + (2 * qp - qg)
    }
    
    # - for pwr>2, start with gamma
    if ( pwr > 2 ) start <- qg
    
    # Solve!
    if ( ( pwr > 1) & (pwr < 2) ) { # This gives a *lower* bound on the value of the answer (if y>0)
      step <- dtweedie(y = 0, 
                       mu = mu.1, 
                       phi = phi.1, 
                       power = pwr)
      # This is P(Y = 0), the discrete "step"
      
      if ( prob <= step ) {
        ans[i] <- 0
      }
    }
    
    if ( is.na(ans[i]) ) { # All cases except Y=0 when 1 < pwr < 2
      
      
      pt2 <- function( q, 
                       mu, 
                       phi, 
                       pwr, 
                       p.given = prob ){ 
        
        ptweedie(q = q, 
                 mu = mu, 
                 phi = phi, 
                 power = pwr ) - p.given
      }
      
      pt <- pt2( q = start, 
                 mu = mu.1, 
                 phi = phi.1, 
                 pwr = pwr,
                 p.given = prob)
      
      #      cat("*** Start   =",start,"; pt =",pt,"\n")
      
      if ( pt == 0 ) ans2[i] <- start
      
      if ( pt > 0 ) { 
        loop <- TRUE
        start.2 <- start
        #start.2 <- 0   # Perhaps set this if too many attempts otherwise
        while ( loop ) {
          # Try harder
          start.2 <- 0.5 * start.2
          if (pt2( q = start.2, 
                   mu.1, 
                   phi.1, 
                   pwr, 
                   p.given = prob )<0 ) loop = FALSE
          #      cat(">>> Start.2 =",start.2,"; pt = ",pt2( q=start.2, mu.1, phi.1, pwr, p.given=prob ),"\n")
          # RECALL:  We are only is this part of the loop if  pt>0
        }
      }
      
      #      cat("*** Start.2 =",start.2,"\n")
      if ( pt < 0) {
        loop <- TRUE
        start.2 <- start
        
        while ( loop ) {
          # Try harder
          start.2 <- 1.5 * (start.2 + 2)
          if (pt2( q = start.2, 
                   mu.1, 
                   phi.1, 
                   pwr, 
                   p.given = prob ) > 0 ) loop = FALSE
          # RECALL:  We are only is this part of the loop if  pt<0
        }
      }
      
      #cat("start, start.2 =",start, start.2,"\n")
      #cat("pt2(start, start.2) =",
      #   pt2(start, mu=mu.1, phi=phi.1, pwr=pwr, p.given=prob), 
      #   pt2(start.2, mu=mu.1, phi=phi.1, pwr=pwr, p.given=prob),"\n")
      
      out <- uniroot(pt2, 
                     c(start, start.2), 
                     mu = mu.1, 
                     phi = phi.1, 
                     p = pwr, 
                     p.given = prob )
      #print(out)
      
      ans[i] <- uniroot(pt2, 
                        c(start, start.2), 
                        mu = mu.1, 
                        phi = phi.1, 
                        p = pwr, 
                        p.given = prob, 
                        tol = 0.000000000001 )$root
    }
    
  }
  
  ans2[ is.na(ans2) ] <-  ans
  ans2
}



#############################################################################
rtweedie <- function(n, xi = NULL, mu, phi, power = NULL){
  
  # Sort out the xi/power notation
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  # Error checks
  if ( any(power < 1) )  stop( paste(index.par.long, "must be greater than 1.\n") )
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( n < 1 ) stop("n must be a positive integer.\n")
  if ( any(mu <= 0) ) stop("mu must be positive.\n")
  if ( length(mu) > 1) {
    if ( length(mu) != n ) stop("mu must be scalar, or of length n.\n")
  } else {
    mu <- array( dim = n, mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != n ) stop("phi must be scalar, or of length n.\n")
  } else {
    phi <- array( dim = n, phi )
    # A vector of all phi's
  }
  if (power==1) {
    rt <- phi * rpois(n, 
                      lambda = mu / ( phi  ) )
  }
  
  if (power == 2) {
    alpha <- (2 - power) / (1 - power)
    gam <- phi * (power - 1) * mu ^ (power - 1)
    rt <- rgamma( n, 
                  shape = 1 / phi, 
                  scale = gam )
  }
  
  if ( power > 2) {
    rt <- qtweedie( runif(n),
                    mu = mu,
                    phi = phi, 
                    power = power)
  }
  
  if ( (power > 1) & (power < 2) ) {
    # Two options:  As above or directly.
    # Directly is faster
    rt <- array( dim = n, NA)
    
    lambda <- mu ^ (2 - power) / ( phi * (2 - power) )
    alpha <- (2 - power) / (1 - power)
    gam <- phi * (power - 1) * mu ^ (power - 1)
    
    N <- rpois(n, 
               lambda = lambda)
    for (i in (1:n) ){
      rt[i] <- rgamma(1, 
                      shape = -N[i] * alpha, 
                      scale = gam[i])
    }
  }
  as.vector(rt)
}


#############################################################################
tweedie.profile <- function(formula, p.vec = NULL, xi.vec = NULL, link.power = 0, 
                            data, weights, offset, fit.glm = FALSE, 
                            do.smooth = TRUE, do.plot = FALSE, 
                            do.ci = do.smooth, eps = 1/6,
                            control = list( epsilon = 1e-09, maxit = glm.control()$maxit, trace = glm.control()$trace ),
                            do.points = do.plot, method = "inversion", conf.level = 0.95, 
                            phi.method = ifelse(method == "saddlepoint", "saddlepoint", "mle"), verbose = FALSE, add0 = FALSE) {
  # verbose gives feedback on screen:
  #    0 : minimal (FALSE)
  #    1 : small amount (TRUE)
  #    2 : lots
  # Determine the value of  p  to use to fit a Tweedie distribution
  # A profile likelihood is used (can can be plotted with do.plot=TRUE)
  
  # The plot can be (somewhat) recovered using (where  out  is the
  # returned object, with smoothing):
  #     plot( out$x, out$y)
  # The actual points can then be added:
  #     points( out$p, out$L)
  # (Note that out$L and out$y are the same if smoothing is  not requested.  
  # Likewise for out$p and out$x.)
  
  # Peter Dunn
  # 07 Dec 2004
  
  
  if ( is.logical( verbose ) ) {
    verbose <- as.numeric(verbose)
  }
  
  if (verbose >= 1 ) {
    cat("---\n This function may take some time to complete;\n")
    cat(" Please be patient.  If it fails, try using  method=\"series\"\n")
    cat(" rather than the default  method=\"inversion\"\n")
    cat(" Another possible reason for failure is the range of p:\n")
    cat(" Try a different input for  p.vec\n---\n")
  }
  
  
  cl <- match.call()
  mf <- match.call()
  m <- match(c("formula", "data", "weights","offset"), names(mf), 0L)
  mf <- mf[c(1, m)]
  mf$drop.unused.levels <- TRUE
  mf[[1]] <- as.name("model.frame")
  mf <- eval(mf, parent.frame())
  mt <- attr(mf, "terms")
  Y <- model.response(mf, "numeric")
  X <- if (!is.empty.model(mt))
    model.matrix(mt, mf, contrasts)
  else matrix(, NROW(Y), 0)
  weights <- as.vector(model.weights(mf))
  if (!is.null(weights) && !is.numeric(weights))
    stop("'weights' must be a numeric vector")
  offset <- as.vector(model.offset(mf))
  if (!is.null(weights) && any(weights < 0))
    stop("negative weights not allowed")
  if (!is.null(offset)) {
    if (length(offset) == 1)
      offset <- rep(offset, NROW(Y))
    else if (length(offset) != NROW(Y))
      stop(gettextf("number of offsets is %d should equal %d (number of observations)",
                    length(offset), NROW(Y)), domain = NA)
  }
  
  
  
  ### NOW some notation stuff
  xi.notation <- TRUE
  if ( is.null(xi.vec) & !is.null(p.vec) ){ # If p.vec given, but not xi.vec
    xi.vec <- p.vec
    xi.notation <- FALSE
  }
  
  if ( is.null(p.vec) & !is.null(xi.vec) ){ # If xi.vec given, but not p.vec
    p.vec <- xi.vec
    xi.notation <- TRUE
  }
  
  if ( is.null( p.vec ) & is.null(xi.vec)) { # Neither xi.vec or p.vec are given
    if ( any(Y == 0) ){
      p.vec <- seq(1.2, 1.8, 
                   by = 0.1)
    } else {
      p.vec <- seq(1.5, 5, 
                   by = 0.5)
    }
    xi.vec <- p.vec
    xi.notation <- TRUE
  }
  
  ### AT THIS POINT, we have both xi.vec and p.vec declared, and both are the same
  ### but we stick with using xi.vec hereafter
  
  
  # Determine notation to use in output (consistent with what was supplied by the user)
  index.par <- ifelse( xi.notation, "xi", "p")
  
  # Find if values of p/xi have been specified that are not appropriate
  xi.fix <- any( xi.vec <= 1 & xi.vec != 0, na.rm = TRUE)
  
  if ( xi.fix ) {
    xi.fix.these <- (xi.vec <= 1 & xi.vec != 0)
    xi.vec[ xi.fix.these ] <- NA 
    if ( length(xi.vec) == 0 ) {
      stop(paste("No values of", index.par, "between 0 and 1, or below 0, are possible.\n"))
    } else {
      cat("Values of", index.par, "between 0 and 1 and less than zero have been removed: such values are not possible.\n")
    }
  }
  
  # Warnings
  if ( any( Y == 0 ) & any( (xi.vec >= 2) | (xi.vec <= 0) ) ) {
    xi.fix.these <- (xi.vec >= 2 | xi.vec <= 0)
    xi.vec[ xi.fix.these ] <- NA 
    if ( length(xi.vec) == 0 ) {
      stop(paste("When the response variable contains exact zeros, all values of", index.par, "must be between 1 and 2.\n"))
    } else {
      cat("When the response variable contains exact zeros, all values of", index.par, "must be between 1 and 2; other values have been removed.\n")
    }
  }
  
  
  # Remove any NAs in the p/xi values
  xi.vec <- xi.vec[ !is.na(xi.vec) ] 
  
  if ( do.smooth & (length(xi.vec) < 5) ) {
    warning(paste("Smoothing needs at least five values of", index.par,".") )
    do.smooth <- FALSE
    do.ci <- FALSE
  }
  if ( (conf.level >= 1) | (conf.level <= 0)  ){
    stop("Confidence level must be between 0 and 1.")
  }
  
  if ( !do.smooth & do.ci ) {
    do.ci <- FALSE
    warning("Confidence intervals only computed if  do.smooth=TRUE\n")
  }
  
  if ( add0 ) {
    p.vec  <- c( 0, p.vec )
    xi.vec <- c( 0, xi.vec )
  }
  
  
  cat(xi.vec,"\n")
  
  
  # Some renaming
  ydata <- Y
  model.x <- X
  
  # Now, fit models!  Need the Tweedie class
  
  # First define the function to *minimize*;
  # since we want a *maximum* likelihood estimator,
  # define the negative.
  dtweedie.nlogl <- function(phi, y, mu, power) {
    ans <- -2 * sum( log( dtweedie( y = y, 
                                    mu = mu, 
                                    phi = phi, 
                                    power = power ) ) )
    if ( is.infinite( ans ) ) {
      # If infinite, replace with saddlepoint estimate?
      ans <- sum( tweedie.dev(y = y, 
                              mu = mu, 
                              power = power) ) / length( y )
    }
    attr(ans, "gradient") <- dtweedie.dldphi(y = y, 
                                             mu = mu, 
                                             phi = phi, 
                                             power = power)
    
    ans
  }
  
  
  # Set up some parameters
  xi.len <- length(xi.vec)
  phi <- NaN
  
  L <- array( dim = xi.len )
  phi.vec <- L
  
  # Now most of these are for debugging:
  b.vec  <- L
  c.vec  <- L
  mu.vec <- L
  b.mat  <- array( dim = c(xi.len, length(ydata) ) )
  
  for (i in (1:xi.len)) {
    
    if ( verbose > 0) {
      cat( paste(index.par," = ", xi.vec[i], "\n", sep="") )
    } else {
      cat(".")
    }

    # We find the mle of phi.
    # We try to bound it as well as we can, 
    # then use  uniroot to solve for it.
    
    # Set things up
    p <- xi.vec[i]
    
    phi.pos <- 1.0e2
    bnd.neg <- -Inf
    bnd.pos <- Inf
    
    if (verbose == 2) cat("* Fitting initial model:")
    
    # Fit the model with given p
    catch.possible.error <- try(
      fit.model <- glm.fit( x = model.x, 
                            y = ydata, 
                            weights = weights, 
                            offset = offset,
                            control = control,
                            family = statmod::tweedie(var.power = p, 
                                                      link.power = link.power)),
      silent = TRUE
    )
    
    
    skip.obs <- FALSE
    if ( is( catch.possible.error, "try-error" ) ) {
      skip.obs <- TRUE 
    }
    
    if( skip.obs ) {
      warning(paste("  Problem near ", index.par, " = ", p,"; this error reported:\n     ",
                    catch.possible.error,
                    " Examine the data and function inputs carefully.") )
      
      # NOTE:  We need epsilon to be very small to make a
      #        smooth likelihood plot
      mu <- rep(NA, length(ydata) )
      
    } else {
      mu <- fitted( fit.model )
    }
    if (verbose == 2) cat(" Done\n")
    
    ### -- Start:  Estimate phi
    if (verbose >= 1) cat("* Phi estimation, method: ", phi.method)
    
    if( skip.obs ) {
      if (verbose >= 1) cat("; but skipped for this obs\n")
      phi.est <- phi <- phi.vec[i] <- NA
    } else {
      if ( phi.method == "mle"){
        
        if (verbose >= 1) cat(" (using optimize): ")
        
        # Saddlepoint approx of phi:
        phi.saddle <- sum( tweedie.dev(y = ydata, 
                                       mu = mu, 
                                       power = p) ) / length( ydata )
        
        if ( is.nan(phi) ) {
          
          # NOTE:  This is only used for first p.
          #        The remainder use the previous estimate of phi as a starting point
          
          # Starting point for phi:  use saddlepoint approx/EQL.
          phi.est <- phi.saddle
          
        } else {
          phi.est <- phi
        }
        low.limit <- min( 0.001, phi.saddle / 2)
        
        #    ans <- nlm(p=phi.est, f=dtweedie.nlogl, 
        #                hessian=FALSE,
        #                power=p, mu=mu, y=data)
        if ( p != 0 ) {
          ans <- optimize(f = dtweedie.nlogl, 
                          maximum = FALSE, 
                          interval = c(low.limit, 
                                       10*phi.est),
                          power = p, 
                          mu = mu, 
                          y = ydata )
          phi <- phi.vec[i] <- ans$minimum
        } else {
          phi <- phi.vec[i] <- sum( (ydata - mu) ^ 2 ) / length(ydata)
        }
        if (verbose >= 1) cat(" Done (phi =", phi.vec[i], ")\n")
        
      } else{ # phi.method=="saddlepoint")
        
        if (verbose >= 1) cat(" (using mean deviance/saddlepoint): ")
        phi <- phi.est <- phi.vec[i] <- sum( tweedie.dev(y = ydata, 
                                                         mu = mu, 
                                                         power = p) )/ length( ydata )
        if (verbose >= 1) cat(" Done (phi =", phi, ")\n")
      }
    }
    
    ### -- Done:  estimate phi
    

    # Now determine log-likelihood at this p
    # Best to use the same density evaluation for all:  series or inversion
    
    if (verbose >= 1) cat("* Computing the log-likelihood ")
    
    # Now compute the log-likelihood
    if (verbose >= 1) cat("(method =", method, "):")
    
    if ( skip.obs ) {
      if (verbose >= 1) cat(" but skipped for this obs\n")
      L[i] <- NA
    } else {   
      if ( method == "saddlepoint") {
        L[i] <- dtweedie.logl.saddle(y = ydata, 
                                     mu = mu, 
                                     power = p, 
                                     phi = phi, 
                                     eps = eps)
      } else {
        if (p == 2) {
          L[i] <- sum( log( dgamma( rate = 1 / (phi * mu), 
                                    shape = 1 / phi, 
                                    x = ydata ) ) )
        } else {
          if ( p == 1 ) {
            if ( phi == 1 ){
              # The phi==1 part added 30 Sep 2010.
              L[i] <- sum( log( dpois(x = ydata / phi, 
                                      lambda = mu / phi ) ) )
            } else { # p=1, but phi is not equal to zero
              # As best as we can, we determine if the values
              # of  ydata  are multiples of  phi
              y.on.phi <- ydata/phi
              close.enough <- array( dim = length(y.on.phi))
              for (i in (1 : length(y.on.phi))){
                if (isTRUE(all.equal(y.on.phi, 
                                     as.integer(y.on.phi)))){
                  L[i] <- sum( log( dpois(x = round(y / phi), 
                                          lambda = mu / phi ) ) )
                } else {
                  L[i] <- 0
                }
              }
            }
          } else {
            if ( p == 0 ) {
              L[i] <- sum( dnorm(x = ydata, 
                                 mean = mu, 
                                 sd = sqrt(phi), 
                                 log = TRUE) )
            } else {
              L[i] <- switch(
                pmatch(method, 
                       c("interpolation", "series", "inversion"), 
                       nomatch = 2),
                "1" = dtweedie.logl( mu = mu, 
                                     power = p, 
                                     phi = phi, 
                                     y = ydata),
                "2" = sum( log( dtweedie.series( y = ydata, 
                                                 mu = mu, 
                                                 power =p, 
                                                 phi =phi) ) ),
                "3" = sum( log( dtweedie.inversion( y = ydata, 
                                                    mu = mu, 
                                                    power = p, 
                                                    phi = phi) ) ) )
            }
          }
        }
      }
    }
    
    if (verbose >= 1) {
      cat(" Done: L =", L[i], "\n")
    } 
  }
  

  if ( verbose == 0 ) cat("Done.\n")
  ### Smoothing
  # If there are infs etc in the log-likelihood,
  # the smooth can't happen.  Perhaps this can be worked around.
  # But at present, we note when it happens to ensure no future errors
  # (eg in computing confidence interval for p).
  
  # y and x are the smoothed plot produced; here we set them
  # as NA so they are defined when the function is returned to
  # the workspace
  y <- NA
  x <- NA
  
  if ( do.smooth ) {
    L.fix <- L
    xi.vec.fix <- xi.vec
    phi.vec.fix <- phi.vec
    if ( any( is.nan(L) ) | any( is.infinite(L) ) | any( is.na(L) ) ) {
      retain.these <- !( ( is.nan(L) ) | ( is.infinite(L) ) | ( is.na(L) ) )
      xi.vec.fix <- xi.vec.fix[ retain.these ]
      phi.vec.fix <- phi.vec.fix[ retain.these ]
      L.fix <- L.fix[ retain.these ]
      
      #         p.vec.fix <- p.vec.fix[ !is.infinite(L.fix) ]
      #         phi.vec.fix <- phi.vec.fix[ !is.infinite(L.fix) ]
      #         L.fix <- L.fix[ !is.infinite(L.fix) ]
      
      if (verbose >= 1) cat("Smooth perhaps inaccurate--log-likelihood contains  Inf  or  NA.\n")
    }
    #else {
    
    if ( length( L.fix ) > 0 ) {
      if (verbose >= 1) cat(".")
      if (verbose >= 1) cat(" --- \n")
      if (verbose >= 1) cat("* Smoothing: ")
      # Smooth the points
      # - get smoothing spline
      ss <- splinefun( xi.vec.fix, 
                       L.fix )
      
      # Plot smoothed data
      xi.smooth <- seq(min(xi.vec.fix), 
                       max(xi.vec.fix), 
                       length = 50 )
      L.smooth <- ss(xi.smooth )
      
      if ( do.plot) {
        keep.these <- is.finite(L.smooth) & !is.na(L.smooth)
        L.smooth <- L.smooth[ keep.these ] 
        xi.smooth <- xi.smooth[ keep.these ] 
        if ( verbose>=1 & any( !keep.these ) ) {
          cat(" (Some values of L are infinite or NA for the smooth; these are ignored)\n")
        }
        
        yrange <- range( L.smooth, na.rm = TRUE )
        
        plot( yrange ~ range(xi.vec),
              type = "n",
              las = 1,
              xlab = ifelse(xi.notation, 
                            expression(paste( xi, " index")), 
                            expression(paste( italic(p), " index")) ),
              ylab = expression(italic(L)))
        lines( xi.smooth, L.smooth,
               lwd = 2)
        rug( xi.vec )
        if (do.points) {
          points( L ~ xi.vec, 
                  pch = 19)
        }
        if (add0) lines(xi.smooth[xi.smooth < 1], 
                        L.smooth[xi.smooth < 1], 
                        col = "gray",
                        lwd = 2)
      }
      x <- xi.smooth
      y <- L.smooth
      
    } else {
      cat("  No valid values of the likelihood computed: smooth aborted\n",
          "  Consider trying another value for the input  method.\n")
    }
  } else {
    if ( do.plot) {
      
      keep.these <- is.finite(L) & !is.na(L)
      xi.vec <- xi.vec[ keep.these ] 
      L <- L[ keep.these ] 
      phi.vec <- phi.vec[ keep.these ]
      if ( verbose >= 1 & any( keep.these ) ) {
        cat(" Some values of L are infinite or NA, and hence ignored\n")
      }
      
      yrange <- range( L, na.rm = TRUE )
      # Plot the data we have
      plot( yrange ~ range(xi.vec),
            type = "n",     
            las = 1,
            xlab = ifelse( xi.notation, 
                           expression(paste(xi, " index")), 
                           expression(paste(italic(p), " index")) ),
            ylab=expression(italic(L)))
      lines( L ~ xi.vec, 
             lwd = 2)
      
      rug( xi.vec )
      if (do.points) {
        points( L ~ xi.vec, 
                pch = 19)
      }
      
    }
    x <- xi.vec
    y <- L
  }
  if (verbose >= 2) cat(" Done\n")
  
  ### Maximum likelihood estimates
  if ( do.smooth ){
    # WARNING:  This spline does not necessarily pass
    #           through the given points.  This should
    #           be seen as a deficiency in the method.
    
    if (verbose >= 2) cat(" Estimating phi:  ")
    
    # Find maximum from profile
    
    L.max <- max(y, na.rm = TRUE)
    xi.max <- x[ y == L.max ]
    
    # In some unusual cases, when add0=TRUE, the maximum of p/xi
    # may be between 0 and 1.
    # In this situation, we... set the mle to either 0 or 1,
    # whichever gives the larger values of the log-likelihood
    if ( (xi.max > 0)  & ( xi.max < 1 ) ) {
      L.max <- max( c( y[xi.vec == 0], 
                       y[xi.vec == 1]) )
      xi.max <- xi.vec[ L.max == y ]
      cat("MLE of", index.par, "is between 0 and 1, which is impossible.",
          "Instead, the MLE of", index.par, "has been set to", xi.max,
          ".  Please check your data and the call to tweedie.profile().")
    }
    
    # Now need to find mle of  phi  at this very value of  p
    # - Find bounds
    
    phi.1 <-   2 * max( phi.vec.fix, na.rm = TRUE )
    phi.2 <- 0.5 * min( phi.vec.fix, na.rm = TRUE )
    
    if ( phi.1 > phi.2 ) {
      phi.hi <- phi.1
      phi.lo <- phi.2 
    } else {
      phi.hi <- phi.2
      phi.lo <- phi.1
    }
    
    if (verbose >= 2) cat(" Done\n")
    
    # Now, if the maximum happens to be at an endpoint,
    # we have to do things differently:
    if ( (xi.max == xi.vec[1]) | (xi.max == xi.vec[length(xi.vec)]) ) {
      
      if ( xi.max == xi.vec[1]) phi.max <- phi.vec[1]
      if ( xi.max == xi.vec[length(xi.vec)]) phi.max <- phi.vec[length(xi.vec)]
      #      phi.max <- phi.lo # and is same as phi.max
      warning("True maximum possibly not detected.")
      
    } else {
      
      # - solve
      
      if ( phi.method == "saddlepoint"){
        
        mu <- fitted( glm.fit( y = ydata, 
                               x = model.x, 
                               weights = weights, 
                               offset = offset,
                               family = statmod::tweedie(xi.max, 
                                                         link.power = link.power)))
        phi.max <- sum( tweedie.dev(y = ydata, 
                                    mu = mu, 
                                    power = xi.max) ) / length( ydata )
        
      } else {
        #        phi.max <- nlm(p=phi.est, f=dtweedie.nlogl, 
        #                    hessian=FALSE,
        #                    power=p, mu=mu, y=data, 
        #                    gradient=dtweedie.dldphi)$estimate
        mu <- fitted( glm.fit( y = ydata, 
                               x = model.x, 
                               weights = weights, 
                               offset = offset,
                               family = statmod::tweedie(xi.max, 
                                                         link.power = link.power)))
        phi.max <- optimize( f = dtweedie.nlogl, 
                             maximum = FALSE, 
                             interval = c(phi.lo, phi.hi ), 
                             # set lower limit phi.lo???
                             power = xi.max, 
                             mu = mu,
                             y = ydata)$minimum
        
        #       Note that using the Hessian often produces computational problems
        #       for little (no?) advantage.
        
      }
    }
  } else {
    # Find maximum from computed data
    if (verbose >= 2) cat(" Finding maximum likelihood estimates:  ")
    
    L.max <- max(L)
    
    xi.max   <- xi.vec  [ L == L.max ]
    phi.max <- phi.vec[ L == L.max ]
    
    if (verbose >= 2) cat(" Done\n")
  }
  
  # Now report
  if ( verbose >= 2 ) {
    cat( "ML Estimates:  ", index.par, "=",xi.max, " with phi=", phi.max, " giving L=", L.max, "\n")
    cat(" ---\n")
  }
  
  # Now find approximate, nominal 95% CI info
  ht <- L.max - ( qchisq(conf.level, 1) / 2 )
  ci <- array( dim = 2, NA )
  
  
  if ( do.ci ) {
    if (verbose == 2) cat("* Finding confidence interval:")
    if ( !do.smooth ) {
      warning("Confidence interval may be very inaccurate without smoothing.\n")
      y <- L
      x <- xi.vec
    }
    
    if ( do.plot ) {
      abline(h = ht, 
             lty = 2)
      title( sub=paste("(", 100*conf.level, "% confidence interval)",
                       sep = "") )
    }
    
    # Now find limits on p:
    # - fit a smoothing spline the other way:  based on L predicting p
    #   so we can ensure that the limits on  p  are found OK
    
    # --- Left side ---
    cond.left <- (y < ht ) & (x < xi.max )
    if ( all(cond.left == FALSE) ) {
      warning("Confidence interval cannot be found: insufficient data to find left CI.\n")
    }else{
      # Now find it!
      approx.left <- max( x[cond.left] )
      index.left <- seq(1, length(cond.left) )
      index.left <- index.left[x == approx.left]
      left.left <- max( 1, index.left - 5)
      left.right <- min( length(cond.left), 
                         index.left + 5)
      
      # Left-side spline
      ss.left <- splinefun( y[left.left:left.right],
                            x[left.left: left.right] )
      ci.new.left <- ss.left( ht )
      
      ci[1] <- ci.new.left
    }
    
    # --- Right side ---
    cond.right <- (y < ht ) & (x > xi.max )
    if ( all( cond.right == FALSE ) ) {
      warning("Confidence interval cannot be found: insufficient data to find right CI.\n")
    }else{
      approx.right <- min( x[cond.right] )
      index.right <- seq(1, length(cond.right) )
      index.right <- index.right[x == approx.right]
      right.left <- max( 1, 
                         index.right - 5)
      right.right <- min( length(cond.left), 
                          index.right + 5)
      
      # Right-side spline
      ss.right <- splinefun( y[right.left : right.right], 
                             x[right.left : right.right] )
      ci.new.right <- ss.right(ht )
      
      ci[2] <- ci.new.right
    }
    if (verbose==2) cat(" Done\n")
  }
  
  if ( fit.glm ) {
    out.glm <- glm.fit( x = model.x, 
                        y = ydata, 
                        weights = weights, 
                        offset = offset,
                        family = statmod::tweedie(var.power = xi.max, 
                                                  link.power = link.power) )
    
    if ( xi.notation){
      out <- list( y = y,
                   x = x,
                   ht = ht,
                   L = L, 
                   xi = xi.vec, 
                   xi.max = xi.max,
                   L.max = L.max, 
                   phi = phi.vec, 
                   phi.max = phi.max, 
                   ci = ci,
                   method = method, 
                   phi.method = phi.method,
                   glm.obj = out.glm)
    } else {
      out <- list( y = y, 
                   x = x, 
                   ht = ht, 
                   L = L, 
                   p = p.vec, 
                   p.max = xi.max, 
                   L.max = L.max, 
                   phi = phi.vec, 
                   phi.max = phi.max, 
                   ci = ci,
                   method = method, 
                   phi.method = phi.method,
                   glm.obj = out.glm)
    }
  } else {
    if (xi.notation ){
      out <- list( y = y, 
                   x = x, 
                   ht = ht, 
                   L = L, 
                   xi = xi.vec,
                   xi.max = xi.max, 
                   L.max = L.max, 
                   phi = phi.vec, 
                   phi.max = phi.max, 
                   ci = ci,
                   method = method, 
                   phi.method = phi.method)
      
    } else {
      out <- list( y = y, 
                   x = x, 
                   ht = ht, 
                   L = L, 
                   p = p.vec, 
                   p.max = xi.max, 
                   L.max = L.max, 
                   phi = phi.vec, 
                   phi.max = phi.max, 
                   ci = ci,
                   method = method, 
                   phi.method = phi.method)
    }
    
  }
  
  if ( verbose ) cat("\n")
  
  invisible( out )
  
}



#############################################################################
dtweedie.stable <- function(y, power, mu, phi)
{
  # Error checks
  if ( power < 1) stop("power must be greater than 2.\n")
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(y < 0) ) stop("y must be a non-negative vector.\n")
  if ( any(mu <= 0) ) stop("mu must be positive.\n")
  if ( length(mu) > 1) {
    if ( length(mu) != length(y) ) stop("mu must be scalar, or the same length as y.\n")
  } else {
    mu <- array( dim = length(y), mu )
    # A vector of all mu's
  }
  if ( length(phi) > 1) {
    if ( length(phi) != length(y) ) stop("phi must be scalar, or the same length as y.\n")
  } else {
    phi <- array( dim = length(y), phi )
    # A vector of all phi's
  }
  
  density <- y
  alpha <- (2 - power) / (1 - power)
  beta <- 1
  k <- 1  # The parameterization used
  delta <- 0
  gamma <- phi * (power - 1) * 
    ( 1 / (phi * (power - 2)) * cos( alpha * pi / 2 ) ) ^ (1 / alpha)
  
  #        require(stabledist) # Needed for  dstable
  ds <- stabledist::dstable(y,
                            alpha = alpha, 
                            beta = beta, 
                            gamma = gamma, 
                            delta = delta, 
                            pm = k)
  density <- exp((y * mu ^ (1 - power) / (1 - power) - mu ^ (2 - power) / (2 - power)) / phi) * ds
  
  
  density
}


#############################################################################
tweedie.plot <- function(y, xi = NULL, mu, phi, type = "pdf", power = NULL, 
                         add =FALSE, ...) {
   
  # Sort out the xi/power notation
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  if ( ( power < 0 ) | ( ( power > 0 ) & ( power < 1 ) ) ) {
    stop( paste("Plots cannot be produced for", index.par.long, "=", power, "\n") )
  }   
  
  is.pg <- ( power > 1 ) & ( power < 2 )
  
  if ( type == "pdf") {
    fy <- dtweedie( y = y, 
                    power = power, 
                    mu = mu,
                    phi = phi)
  } else {
    fy <- ptweedie( q = y, 
                    power = power, 
                    mu = mu,
                    phi = phi)
  }
  
  if ( !add ) {
    if ( is.pg ) {
      plot( range(fy) ~ range( y ), 
            type = "n", ...)
      if ( any( y == 0 ) ) { # The exact zero
        points( fy[y == 0] ~ y[y == 0], 
                pch = 19, ... )
      }
      if ( any( y > 0 ) ) { # The exact zero
        lines( fy[y > 0]   ~ y[y > 0], 
               pch = 19, ... )
      }
    } else {  # Not a Poison-gamma dist
      plot( range(fy) ~ range( y ), 
            type = "n", ...)
      lines( fy ~ y, 
             pch = 19, ... )
    }
  } else {# Add; no new plot
    if ( is.pg ) {
      if ( any( y == 0 ) ) { # The exact zero
        points( fy[y == 0] ~ y[y == 0], 
                pch = 19, ... )
      }
      if ( any( y > 0 ) ) { # The exact zero
        lines( fy[y > 0]   ~ y[y > 0], 
               pch = 19, ... )
      }
    } else {  # Not a Poison-gamma dist
      lines( fy ~ y, 
             pch = 19, ... )
    }
    
  }
  return(invisible(list(y = fy, 
                        x = y) ))
  
}
#############################################################################
AICtweedie <- function( glm.obj, dispersion = NULL, k = 2, verbose = TRUE){ 
  # New  dispersion  input for (e.g.) Poisson case, added October 2017
  
  wt <- glm.obj$prior.weights
  n <- length(glm.obj$residuals)
  edf <- glm.obj$rank  # As used in logLik.glm()
  
  mu <- fitted( glm.obj )
  y  <- glm.obj$y
  p <- get("p", envir = environment(glm.obj$family$variance))
  
  if ( is.null(dispersion)) {  # New section
    if (p == 1 & verbose) message("*** Tweedie index power = 1: Consider using  dispersion=1  in call to  AICtweedie().\n")
    dev <- deviance(glm.obj)
    disp <- dev / sum(wt)  # In line with Gamma()$aic
    edf <- edf + 1  # ADD one as we are estimating phi too
  } else {
    disp <- dispersion
  }
  
  den <- dtweedie( y = y, 
                   mu = mu, 
                   phi = disp, 
                   power = p)
  AIC <- -2 * sum( log(den) * wt) 
  
  return( AIC + k * (edf) )
  
}

#############################################################################

tweedie.convert <- function(xi = NULL, mu, phi, power = NULL){
  ### ADDED 14 July 2017
  
  if ( is.null(power) & is.null(xi) ) stop("Either xi or power must be given\n")
  xi.notation <- TRUE
  if ( is.null(power) ) {   # Then  xi  is given
    if ( !is.numeric(xi)) stop("xi  must be numeric.\n")
    power <- xi
  } else {
    xi.notation <- FALSE
  }
  if ( is.null(xi) ) {   # Then   power  is given
    if ( !is.numeric(xi)) stop("power  must be numeric.\n")
    xi.notation <- FALSE
    xi <- power
  }
  if ( xi != power ) {
    cat("Different values for xi and power given; the value of xi used.\n")
    power <- xi
  }
  index.par       <- ifelse( xi.notation, "xi", "p")
  index.par.long  <- ifelse( xi.notation, "xi", "power")
  
  
  # Error checks
  if ( power < 1)      stop( paste(index.par.long, "must be greater than 1.\n") )
  if ( power >= 2)     stop( paste(index.par.long, "must be less than 2.\n") )
  if ( any(phi <= 0) ) stop("phi must be positive.")
  if ( any(mu <= 0) )  stop("mu must be positive.\n")
  
  if( length(mu) != length(phi) ){
    if ( length(mu)  == 1 ) mu  <- array(dim = length(phi), mu  ) 
    if ( length(phi) == 1 ) phi <- array(dim = length(mu),  phi ) 
  }
  # Now  mu  and  phi  will be the same length if one of them happened to be a scalar, so this works:
  if( length(mu) != length(phi) ) stop("phi and mu must be scalars, or the same length.\n")
  
  lambda <- ( mu ^ (2 - xi) ) / ( phi * (2 - xi) )  # Poisson distribution mean
  alpha  <- (2 - xi)  / (xi - 1)           # gamma distribution alpha (shape)
  gam    <- phi * (xi - 1) * mu ^ (xi - 1)   # gamma distribution beta  (scale)
  p0     <- exp( -lambda )
  phi.g  <- (2 - xi) * (xi - 1) * phi ^ 2 * mu ^ ( 2 * (xi-1) )
  mu     <- gam / phi
  
  list( poisson.lambda = lambda, 
        gamma.shape = alpha, 
        gamma.scale = gam, 
        p0 = p0,
        gamma.mean = mu,
        gamma.phi = phi.g)
  
}