##################################################################################### # Jim Gatheral, 2023 # Various rational approximations to the rough Heston solution # # Code simplified and speeded up; Padé66 added, February 2024 ##################################################################################### ######################################################################## # Padé approximations to h(a,tau) ######################################################################## ######################################################################## h.Pade22 <- function (params)function(a, tau) { H <- params$H rho <- params$rho nu <- params$nu al <- H + 1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp - (0 + (0+1i)) * rho * a aa <- sqrt(a * (a + (0 + (0+1i))) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa b1 <- -a * (a + (0+1i))/2 * 1/gamma(1 + al) b2 <- -b1 * lamTilde * nu * gamma(1 + al)/gamma(1 + 2 * al) g0 <- rm/nu g1 <- ifelse(al == 1, 0, -1/aa * 1/gamma(1 - al) * g0/nu) den <- g0^2 + b1 * g1 q1 <- (b1 * g0 - b2 * g1)/den q2 <- (b1^2 + b2 * g0)/den p1 <- b1 p2 <- b2 + b1*q1 y <- tau^al h.pade <- (p1 * y + p2 * y^2)/(1 + q1 * y + q2 * y^2) return(h.pade) } ######################################################################## h.Pade33 <- function(params)function(a,tau){ H <- params$H rho <- params$rho nu <- params$nu al <- H+1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp -(0+1i) * rho * a aa <- sqrt(a * (a + (0+1i)) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa b1 <- -a*(a+1i)/2 * 1/gamma(1+al) b2 <- -b1 * lamTilde * nu * gamma(1+al)/gamma(1+2*al) b3 <- (-b2 * lamTilde * nu + nu^2 * b1^2/2) * gamma(1+2*al)/gamma(1+3*al) b4 <- (-b3 * lamTilde *nu + nu^2 * b1*b2)* gamma(1+3*al)/gamma(1+4*al) g0 <- rm/nu g1 <- -1/(aa*nu)*1/gamma(1-al)*g0 g2 <- -1/(aa*nu)*(gamma(1-al)/gamma(1-2*al)*g1 -1/2*nu^2*g1*g1) g3 <- -1/(aa*nu)*(gamma(1-2*al)/gamma(1-3*al)*g2 -nu^2 * g1*g2) den <- g0^3 +2*b1*g0*g1-b2*g1^2+b1^2*g2+b2*g0*g2 p1 <- b1 p2 <- (b1^2*g0^2 + b2*g0^3 + b1^3*g1 + b1*b2*g0*g1 - b2^2*g1^2 +b1*b3*g1^2 +b2^2*g0*g2 - b1*b3*g0*g2)/den q1 <- (b1*g0^2 + b1^2*g1 - b2*g0*g1 + b3*g1^2 - b1*b2*g2 -b3*g0*g2)/den q2 <- (b1^2*g0 + b2*g0^2 - b1*b2*g1 - b3*g0*g1 + b2^2*g2 - b1*b3*g2)/den q3 <- (b1^3 + 2*b1*b2*g0 + b3*g0^2 -b2^2*g1 +b1*b3*g1 )/den p3 <- g0*q3 y <- tau^al h.pade <- (p1*y + p2*y^2 + p3*y^3)/(1 + q1*y + q2*y^2 + q3*y^3) return(h.pade) } ######################################################################## h.Pade44 <- function (params) function(a, x) { H <- params$H rho <- params$rho nu <- params$nu al <- H + 1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp - (0 + (0+1i)) * rho * a aa <- sqrt(a * (a + (0 + (0+1i))) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa b1 <- -a * (a + (0+1i))/2 * 1/gamma(1 + al) b2 <- -b1 * lamTilde * nu * gamma(1 + al)/gamma(1 + 2 * al) b3 <- (-b2 * lamTilde * nu + nu^2 * b1^2/2) * gamma(1 + 2 * al)/gamma(1 + 3 * al) b4 <- (-b3 * lamTilde * nu + nu^2 * b1 * b2) * gamma(1 + 3 * al)/gamma(1 + 4 * al) g0 <- rm/nu g1 <- ifelse(al == 1, 0, -1/(aa * nu) * 1/gamma(1 - al) * g0) g2 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - al)/gamma(1 - 2 * al) * g1 - 1/2 * nu^2 * g1 * g1)) g3 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - 2 * al)/gamma(1 - 3 * al) * g2 - nu^2 * g1 * g2)) den <- (g0^4 + 3 * b1 * g0^2 * g1 + b1^2 * g1^2 - 2 * b2 * g0 * g1^2 + b3 * g1^3 + 2 * b1^2 * g0 * g2 + 2 * b2 * g0^2 * g2 - 2 * b1 * b2 * g1 * g2 - 2 * b3 * g0 * g1 * g2 + b2^2 * g2^2 - b1 * b3 * g2^2 + b1^3 * g3 + 2 * b1 * b2 * g0 * g3 + b3 * g0^2 * g3 - b2^2 * g1 * g3 + b1 * b3 * g1 * g3) q1 <- (b1 * g0^3 + 2 * b1^2 * g0 * g1 - b2 * g0^2 * g1 - 2 * b1 * b2 * g1^2 + b3 * g0 * g1^2 - b4 * g1^3 + b1^3 * g2 - b3 * g0^2 * g2 + b2^2 * g1 * g2 + b1 * b3 * g1 * g2 + 2 * b4 * g0 * g1 * g2 - b2 * b3 * g2^2 + b1 * b4 * g2^2 - b1^2 * b2 * g3 - b2^2 * g0 * g3 - b1 * b3 * g0 * g3 - b4 * g0^2 * g3 + b2 * b3 * g1 * g3 - b1 * b4 * g1 * g3)/den q2 <- (b1^2 * g0^2 + b2 * g0^3 + b1^3 * g1 - b3 * g0^2 * g1 + b1 * b3 * g1^2 + b4 * g0 * g1^2 - b1^2 * b2 * g2 + b2^2 * g0 * g2 - 3 * b1 * b3 * g0 * g2 - b4 * g0^2 * g2 - b2 * b3 * g1 * g2 + b1 * b4 * g1 * g2 + b3^2 * g2^2 - b2 * b4 * g2^2 + b1 * b2^2 * g3 - b1^2 * b3 * g3 + b2 * b3 * g0 * g3 - b1 * b4 * g0 * g3 - b3^2 * g1 * g3 + b2 * b4 * g1 * g3)/den q3 <- (b1^3 * g0 + 2 * b1 * b2 * g0^2 + b3 * g0^3 - b1^2 * b2 * g1 - 2 * b2^2 * g0 * g1 - b4 * g0^2 * g1 + b2 * b3 * g1^2 - b1 * b4 * g1^2 + b1 * b2^2 * g2 - b1^2 * b3 * g2 + b2 * b3 * g0 * g2 - b1 * b4 * g0 * g2 - b3^2 * g1 * g2 + b2 * b4 * g1 * g2 - b2^3 * g3 + 2 * b1 * b2 * b3 * g3 - b1^2 * b4 * g3 + b3^2 * g0 * g3 - b2 * b4 * g0 * g3)/den q4 <- (b1^4 + 3 * b1^2 * b2 * g0 + b2^2 * g0^2 + 2 * b1 * b3 * g0^2 + b4 * g0^3 - 2 * b1 * b2^2 * g1 + 2 * b1^2 * b3 * g1 - 2 * b2 * b3 * g0 * g1 + 2 * b1 * b4 * g0 * g1 + b3^2 * g1^2 - b2 * b4 * g1^2 + b2^3 * g2 - 2 * b1 * b2 * b3 * g2 + b1^2 * b4 * g2 - b3^2 * g0 * g2 + b2 * b4 * g0 * g2)/den p1 <- b1 p2 <- b2 + b1*q1 p3 <- b3 + b1*q2 + b2*q1 p4 <- g0 * q4 y <- x^al h.pade <- (p1 * y + p2 * y^2 + p3 * y^3 + p4 * y^4)/(1 + q1 * y + q2 * y^2 + q3 * y^3 + q4 * y^4) return(h.pade) } ######################################################################## h.Pade55 <- function (params)function(a, tau) { H <- params$H rho <- params$rho nu <- params$nu al <- H + 1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp - (0 + (0+1i)) * rho * a aa <- sqrt(a * (a + (0 + (0+1i))) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa b1 <- -a * (a + (0+1i))/2 * 1/gamma(1 + al) b2 <- -b1 * lamTilde * nu * gamma(1 + al)/gamma(1 + 2 * al) b3 <- (-b2 * lamTilde * nu + nu^2 * b1^2/2) * gamma(1 + 2 * al)/gamma(1 + 3 * al) b4 <- (-b3 * lamTilde * nu + nu^2 * b1 * b2) * gamma(1 + 3 * al)/gamma(1 + 4 * al) b5 <- (-b4 * lamTilde * nu + nu^2 * (1/2 * b2 * b2 + b1 * b3)) * gamma(1 + 4 * al)/gamma(1 + 5 * al) g0 <- rm/nu g1 <- ifelse(al == 1, 0, -1/(aa * nu) * 1/gamma(1 - al) * g0) g2 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - al)/gamma(1 - 2 * al) * g1 - 1/2 * nu^2 * g1 * g1)) g3 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - 2 * al)/gamma(1 - 3 * al) * g2 - nu^2 * g1 * g2)) g4 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - 3 * al)/gamma(1 - 4 * al) * g3 - nu^2 * (1/2 * g2 * g2 + g1 * g3))) den <- (-g0**5 - 4*b1*g0**3*g1 - 3*b1**2*g0*g1**2 + 3*b2*g0**2*g1**2 + 2*b1*b2*g1**3 - 2*b3*g0*g1**3 + b4*g1**4 - 3*b1**2*g0**2*g2 - 3*b2*g0**3*g2 - 2*b1**3*g1*g2 + 2*b1*b2*g0*g1*g2 + 4*b3*g0**2*g1*g2 - b2**2*g1**2*g2 - 2*b1*b3*g1**2*g2 - 3*b4*g0*g1**2*g2 + b1**2*b2*g2**2 - 2*b2**2*g0*g2**2 + 4*b1*b3*g0*g2**2 + b4*g0**2*g2**2 + 2*b2*b3*g1*g2**2 - 2*b1*b4*g1*g2**2 - b3**2*g2**3 + b2*b4*g2**3 - 2*b1**3*g0*g3 - 4*b1*b2*g0**2*g3 - 2*b3*g0**3*g3 + 2*b1**2*b2*g1*g3 + 4*b2**2*g0*g1*g3 + 2*b4*g0**2*g1*g3 - 2*b2*b3*g1**2*g3 + 2*b1*b4*g1**2*g3 - 2*b1*b2**2*g2*g3 + 2*b1**2*b3*g2*g3 - 2*b2*b3*g0*g2*g3 + 2*b1*b4*g0*g2*g3 + 2*b3**2*g1*g2*g3 - 2*b2*b4*g1*g2*g3 + b2**3*g3**2 - 2*b1*b2*b3*g3**2 + b1**2*b4*g3**2 - b3**2*g0*g3**2 + b2*b4*g0*g3**2 - b1**4*g4 - 3*b1**2*b2*g0*g4 - b2**2*g0**2*g4 - 2*b1*b3*g0**2*g4 - b4*g0**3*g4 + 2*b1*b2**2*g1*g4 - 2*b1**2*b3*g1*g4 + 2*b2*b3*g0*g1*g4 - 2*b1*b4*g0*g1*g4 - b3**2*g1**2*g4 + b2*b4*g1**2*g4 - b2**3*g2*g4 + 2*b1*b2*b3*g2*g4 - b1**2*b4*g2*g4 + b3**2*g0*g2*g4 - b2*b4*g0*g2*g4) q1 <- (-(b1*g0**4) - 3*b1**2*g0**2*g1 + b2*g0**3*g1 - b1**3*g1**2 + 4*b1*b2*g0*g1**2 - b3*g0**2*g1**2 - b2**2*g1**3 - 2*b1*b3*g1**3 + b4*g0*g1**3 - b5*g1**4 - 2*b1**3*g0*g2 - b1*b2*g0**2*g2 + b3*g0**3*g2 + 4*b1**2*b2*g1*g2 + 2*b1*b3*g0*g1*g2 - 2*b4*g0**2*g1*g2 + 2*b2*b3*g1**2*g2 + b1*b4*g1**2*g2 + 3*b5*g0*g1**2*g2 - 2*b1*b2**2*g2**2 + b1**2*b3*g2**2 - 2*b1*b4*g0*g2**2 - b5*g0**2*g2**2 - b3**2*g1*g2**2 - b2*b4*g1*g2**2 + 2*b1*b5*g1*g2**2 + b3*b4*g2**3 - b2*b5*g2**3 - b1**4*g3 - b1**2*b2*g0*g3 + b2**2*g0**2*g3 + b4*g0**3*g3 - 2*b1**2*b3*g1*g3 - 4*b2*b3*g0*g1*g3 - 2*b5*g0**2*g1*g3 + b3**2*g1**2*g3 + b2*b4*g1**2*g3 - 2*b1*b5*g1**2*g3 + b2**3*g2*g3 - b1**2*b4*g2*g3 + b3**2*g0*g2*g3 + b2*b4*g0*g2*g3 - 2*b1*b5*g0*g2*g3 - 2*b3*b4*g1*g2*g3 + 2*b2*b5*g1*g2*g3 - b2**2*b3*g3**2 + b1*b3**2*g3**2 + b1*b2*b4*g3**2 - b1**2*b5*g3**2 + b3*b4*g0*g3**2 - b2*b5*g0*g3**2 + b1**3*b2*g4 + 2*b1*b2**2*g0*g4 + b1**2*b3*g0*g4 + 2*b2*b3*g0**2*g4 + b1*b4*g0**2*g4 + b5*g0**3*g4 - b2**3*g1*g4 + b1**2*b4*g1*g4 - b3**2*g0*g1*g4 - b2*b4*g0*g1*g4 + 2*b1*b5*g0*g1*g4 + b3*b4*g1**2*g4 - b2*b5*g1**2*g4 + b2**2*b3*g2*g4 - b1*b3**2*g2*g4 - b1*b2*b4*g2*g4 + b1**2*b5*g2*g4 - b3*b4*g0*g2*g4 + b2*b5*g0*g2*g4)/den q2 <- (-(b1**2*g0**3) - b2*g0**4 - 2*b1**3*g0*g1 - b1*b2*g0**2*g1 + b3*g0**3*g1 + 2*b1**2*b2*g1**2 + b2**2*g0*g1**2 - b4*g0**2*g1**2 + b1*b4*g1**3 + b5*g0*g1**3 - b1**4*g2 - b1**2*b2*g0*g2 - 2*b2**2*g0**2*g2 + 3*b1*b3*g0**2*g2 + b4*g0**3*g2 - 2*b1*b2**2*g1*g2 - 4*b1*b4*g0*g1*g2 - 2*b5*g0**2*g1*g2 - b2*b4*g1**2*g2 + b1*b5*g1**2*g2 + 2*b1*b2*b3*g2**2 - 2*b1**2*b4*g2**2 - b3**2*g0*g2**2 + 3*b2*b4*g0*g2**2 - 2*b1*b5*g0*g2**2 + b3*b4*g1*g2**2 - b2*b5*g1*g2**2 - b4**2*g2**3 + b3*b5*g2**3 + b1**3*b2*g3 + 3*b1**2*b3*g0*g3 + 3*b1*b4*g0**2*g3 + b5*g0**3*g3 + b2**3*g1*g3 - 2*b1*b2*b3*g1*g3 + b1**2*b4*g1*g3 + b3**2*g0*g1*g3 - b2*b4*g0*g1*g3 - b3*b4*g1**2*g3 + b2*b5*g1**2*g3 - b2**2*b3*g2*g3 - b1*b3**2*g2*g3 + 3*b1*b2*b4*g2*g3 - b1**2*b5*g2*g3 - b3*b4*g0*g2*g3 + b2*b5*g0*g2*g3 + 2*b4**2*g1*g2*g3 - 2*b3*b5*g1*g2*g3 + b2*b3**2*g3**2 - b2**2*b4*g3**2 - b1*b3*b4*g3**2 + b1*b2*b5*g3**2 - b4**2*g0*g3**2 + b3*b5*g0*g3**2 - b1**2*b2**2*g4 + b1**3*b3*g4 - b2**3*g0*g4 + b1**2*b4*g0*g4 - b2*b4*g0**2*g4 + b1*b5*g0**2*g4 + b2**2*b3*g1*g4 + b1*b3**2*g1*g4 - 3*b1*b2*b4*g1*g4 + b1**2*b5*g1*g4 + b3*b4*g0*g1*g4 - b2*b5*g0*g1*g4 - b4**2*g1**2*g4 + b3*b5*g1**2*g4 - b2*b3**2*g2*g4 + b2**2*b4*g2*g4 + b1*b3*b4*g2*g4 - b1*b2*b5*g2*g4 + b4**2*g0*g2*g4 - b3*b5*g0*g2*g4)/den q3 <- (-(b1**3*g0**2) - 2*b1*b2*g0**3 - b3*g0**4 - b1**4*g1 - b1**2*b2*g0*g1 + 2*b2**2*g0**2*g1 - b1*b3*g0**2*g1 + b4*g0**3*g1 + b1*b2**2*g1**2 - 2*b1**2*b3*g1**2 - 2*b2*b3*g0*g1**2 - b5*g0**2*g1**2 + b2*b4*g1**3 - b1*b5*g1**3 + b1**3*b2*g2 + 3*b1**2*b3*g0*g2 + 3*b1*b4*g0**2*g2 + b5*g0**3*g2 + 2*b3**2*g0*g1*g2 - 2*b2*b4*g0*g1*g2 - b3*b4*g1**2*g2 + b2*b5*g1**2*g2 - b1*b3**2*g2**2 + b1*b2*b4*g2**2 - b3*b4*g0*g2**2 + b2*b5*g0*g2**2 + b4**2*g1*g2**2 - b3*b5*g1*g2**2 - b1**2*b2**2*g3 + b1**3*b3*g3 + b2**3*g0*g3 - 4*b1*b2*b3*g0*g3 + 3*b1**2*b4*g0*g3 - 2*b3**2*g0**2*g3 + b2*b4*g0**2*g3 + b1*b5*g0**2*g3 - b2**2*b3*g1*g3 + 3*b1*b3**2*g1*g3 - b1*b2*b4*g1*g3 - b1**2*b5*g1*g3 + 3*b3*b4*g0*g1*g3 - 3*b2*b5*g0*g1*g3 - b4**2*g1**2*g3 + b3*b5*g1**2*g3 + b2*b3**2*g2*g3 - b2**2*b4*g2*g3 - b1*b3*b4*g2*g3 + b1*b2*b5*g2*g3 - b4**2*g0*g2*g3 + b3*b5*g0*g2*g3 - b3**3*g3**2 + 2*b2*b3*b4*g3**2 - b1*b4**2*g3**2 - b2**2*b5*g3**2 + b1*b3*b5*g3**2 + b1*b2**3*g4 - 2*b1**2*b2*b3*g4 + b1**3*b4*g4 + b2**2*b3*g0*g4 - 2*b1*b3**2*g0*g4 + b1**2*b5*g0*g4 - b3*b4*g0**2*g4 + b2*b5*g0**2*g4 - b2*b3**2*g1*g4 + b2**2*b4*g1*g4 + b1*b3*b4*g1*g4 - b1*b2*b5*g1*g4 + b4**2*g0*g1*g4 - b3*b5*g0*g1*g4 + b3**3*g2*g4 - 2*b2*b3*b4*g2*g4 + b1*b4**2*g2*g4 + b2**2*b5*g2*g4 - b1*b3*b5*g2*g4)/den q4 <- (-(b1**4*g0) - 3*b1**2*b2*g0**2 - b2**2*g0**3 - 2*b1*b3*g0**3 - b4*g0**4 + b1**3*b2*g1 + 4*b1*b2**2*g0*g1 - b1**2*b3*g0*g1 + 4*b2*b3*g0**2*g1 - b1*b4*g0**2*g1 + b5*g0**3*g1 - b2**3*g1**2 + b1**2*b4*g1**2 - 2*b3**2*g0*g1**2 + 2*b1*b5*g0*g1**2 + b3*b4*g1**3 - b2*b5*g1**3 - b1**2*b2**2*g2 + b1**3*b3*g2 - 2*b2**3*g0*g2 + 2*b1*b2*b3*g0*g2 + b3**2*g0**2*g2 - 2*b2*b4*g0**2*g2 + b1*b5*g0**2*g2 + 2*b2**2*b3*g1*g2 - 4*b1*b2*b4*g1*g2 + 2*b1**2*b5*g1*g2 - b4**2*g1**2*g2 + b3*b5*g1**2*g2 - b2*b3**2*g2**2 + b2**2*b4*g2**2 + b1*b3*b4*g2**2 - b1*b2*b5*g2**2 + b4**2*g0*g2**2 - b3*b5*g0*g2**2 + b1*b2**3*g3 - 2*b1**2*b2*b3*g3 + b1**3*b4*g3 + b2**2*b3*g0*g3 - 2*b1*b3**2*g0*g3 + b1**2*b5*g0*g3 - b3*b4*g0**2*g3 + b2*b5*g0**2*g3 - b2*b3**2*g1*g3 + b2**2*b4*g1*g3 + b1*b3*b4*g1*g3 - b1*b2*b5*g1*g3 + b4**2*g0*g1*g3 - b3*b5*g0*g1*g3 + b3**3*g2*g3 - 2*b2*b3*b4*g2*g3 + b1*b4**2*g2*g3 + b2**2*b5*g2*g3 - b1*b3*b5*g2*g3 - b2**4*g4 + 3*b1*b2**2*b3*g4 - b1**2*b3**2*g4 - 2*b1**2*b2*b4*g4 + b1**3*b5*g4 + 2*b2*b3**2*g0*g4 - 2*b2**2*b4*g0*g4 - 2*b1*b3*b4*g0*g4 + 2*b1*b2*b5*g0*g4 - b4**2*g0**2*g4 + b3*b5*g0**2*g4 - b3**3*g1*g4 + 2*b2*b3*b4*g1*g4 - b1*b4**2*g1*g4 - b2**2*b5*g1*g4 + b1*b3*b5*g1*g4)/den q5 <- (-b1**5 - 4*b1**3*b2*g0 - 3*b1*b2**2*g0**2 - 3*b1**2*b3*g0**2 - 2*b2*b3*g0**3 - 2*b1*b4*g0**3 - b5*g0**4 + 3*b1**2*b2**2*g1 - 3*b1**3*b3*g1 + 2*b2**3*g0*g1 + 2*b1*b2*b3*g0*g1 - 4*b1**2*b4*g0*g1 + b3**2*g0**2*g1 + 2*b2*b4*g0**2*g1 - 3*b1*b5*g0**2*g1 - b2**2*b3*g1**2 - 2*b1*b3**2*g1**2 + 4*b1*b2*b4*g1**2 - b1**2*b5*g1**2 - 2*b3*b4*g0*g1**2 + 2*b2*b5*g0*g1**2 + b4**2*g1**3 - b3*b5*g1**3 - 2*b1*b2**3*g2 + 4*b1**2*b2*b3*g2 - 2*b1**3*b4*g2 - 2*b2**2*b3*g0*g2 + 4*b1*b3**2*g0*g2 - 2*b1**2*b5*g0*g2 + 2*b3*b4*g0**2*g2 - 2*b2*b5*g0**2*g2 + 2*b2*b3**2*g1*g2 - 2*b2**2*b4*g1*g2 - 2*b1*b3*b4*g1*g2 + 2*b1*b2*b5*g1*g2 - 2*b4**2*g0*g1*g2 + 2*b3*b5*g0*g1*g2 - b3**3*g2**2 + 2*b2*b3*b4*g2**2 - b1*b4**2*g2**2 - b2**2*b5*g2**2 + b1*b3*b5*g2**2 + b2**4*g3 - 3*b1*b2**2*b3*g3 + b1**2*b3**2*g3 + 2*b1**2*b2*b4*g3 - b1**3*b5*g3 - 2*b2*b3**2*g0*g3 + 2*b2**2*b4*g0*g3 + 2*b1*b3*b4*g0*g3 - 2*b1*b2*b5*g0*g3 + b4**2*g0**2*g3 - b3*b5*g0**2*g3 + b3**3*g1*g3 - 2*b2*b3*b4*g1*g3 + b1*b4**2*g1*g3 + b2**2*b5*g1*g3 - b1*b3*b5*g1*g3)/den p1 <- b1 p2 <- b2 + b1*q1 p3 <- b3 + b1*q2 + b2*q1 p4 <- b4 + b3*q1 + b2*q2 + b1*q3 p5 <- g0*q5 y <- tau^al h.pade <- (p1 * y + p2 * y^2 + p3 * y^3 + p4 * y^4 + p5 * y^5)/(1 + q1 * y + q2 * y^2 + q3 * y^3 + q4 * y^4 + q5 * y^5) return(h.pade) } ######################################################################## h.Pade66 <- function (params) function(a, x) { H <- params$H rho <- params$rho nu <- params$nu al <- H + 1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp - (0 + (0+1i)) * rho * a aa <- sqrt(a * (a + (0 + (0+1i))) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa b1 <- -a * (a + (0+1i))/2 * 1/gamma(1 + al) b2 <- -b1 * lamTilde * nu * gamma(1 + al)/gamma(1 + 2 * al) b3 <- (-b2 * lamTilde * nu + nu^2 * b1^2/2) * gamma(1 + 2 * al)/gamma(1 + 3 * al) b4 <- (-b3 * lamTilde * nu + nu^2 * b1 * b2) * gamma(1 + 3 * al)/gamma(1 + 4 * al) b5 <- (-b4 * lamTilde * nu + nu^2 * (1/2 * b2 * b2 + b1 * b3)) * gamma(1 + 4 * al)/gamma(1 + 5 * al) b6 <- (-b5 * lamTilde * nu + nu^2 * (b2 * b3 + b1 * b4)) * gamma(1 + 5 * al)/gamma(1 + 6 * al) g0 <- rm/nu g1 <- ifelse(al == 1, 0, -1/(aa * nu) * 1/gamma(1 - al) * g0) g2 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - al)/gamma(1 - 2 * al) * g1 - 1/2 * nu^2 * g1 * g1)) g3 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - 2 * al)/gamma(1 - 3 * al) * g2 - nu^2 * g1 * g2)) g4 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - 3 * al)/gamma(1 - 4 * al) * g3 - nu^2 * (1/2 * g2 * g2 + g1 * g3))) g5 <- ifelse(al == 1, 0, -1/(aa * nu) * (gamma(1 - 4 * al)/gamma(1 - 5 * al) * g4 - nu^2 * (g2 * g3 + g1 * g4))) den <- ( -g0**6 - 5*b1*g0**4*g1 - 6*b1**2*g0**2*g1**2 + 4*b2*g0**3*g1**2 - b1**3*g1**3 + 6*b1*b2*g0*g1**3 - 3*b3*g0**2*g1**3 - b2**2*g1**4 - 2*b1*b3*g1**4 + 2*b4*g0*g1**4 - b5*g1**5 - 4*b1**2*g0**3*g2 - 4*b2*g0**4*g2 - 6*b1**3*g0*g1*g2 + 6*b3*g0**3*g1*g2 + 6*b1**2*b2*g1**2*g2 - 6*b4*g0**2*g1**2*g2 + 2*b2*b3*g1**3*g2 + 2*b1*b4*g1**3*g2 + 4*b5*g0*g1**3*g2 - b1**4*g2**2 - 4*b2**2*g0**2*g2**2 + 7*b1*b3*g0**2*g2**2 + 2*b4*g0**3*g2**2 - 4*b1*b2**2*g1*g2**2 + b1**2*b3*g1*g2**2 + 2*b2*b3*g0*g1*g2**2 - 8*b1*b4*g0*g1*g2**2 - 3*b5*g0**2*g1*g2**2 - b3**2*g1**2*g2**2 - 2*b2*b4*g1**2*g2**2 + 3*b1*b5*g1**2*g2**2 + 2*b1*b2*b3*g2**3 - 2*b1**2*b4*g2**3 - 2*b3**2*g0*g2**3 + 4*b2*b4*g0*g2**3 - 2*b1*b5*g0*g2**3 + 2*b3*b4*g1*g2**3 - 2*b2*b5*g1*g2**3 - b4**2*g2**4 + b3*b5*g2**4 - 3*b1**3*g0**2*g3 - 6*b1*b2*g0**3*g3 - 3*b3*g0**4*g3 - 2*b1**4*g1*g3 + 7*b2**2*g0**2*g1*g3 - b1*b3*g0**2*g1*g3 + 4*b4*g0**3*g1*g3 + b1*b2**2*g1**2*g3 - 4*b1**2*b3*g1**2*g3 - 8*b2*b3*g0*g1**2*g3 + 2*b1*b4*g0*g1**2*g3 - 3*b5*g0**2*g1**2*g3 + b3**2*g1**3*g3 + 2*b2*b4*g1**3*g3 - 3*b1*b5*g1**3*g3 + 2*b1**3*b2*g2*g3 - 2*b1*b2**2*g0*g2*g3 + 8*b1**2*b3*g0*g2*g3 - 2*b2*b3*g0**2*g2*g3 + 8*b1*b4*g0**2*g2*g3 + 2*b5*g0**3*g2*g3 + 2*b2**3*g1*g2*g3 - 2*b1*b2*b3*g1*g2*g3 + 6*b3**2*g0*g1*g2*g3 - 4*b2*b4*g0*g1*g2*g3 - 2*b1*b5*g0*g1*g2*g3 - 4*b3*b4*g1**2*g2*g3 + 4*b2*b5*g1**2*g2*g3 - b2**2*b3*g2**2*g3 - 2*b1*b3**2*g2**2*g3 + 4*b1*b2*b4*g2**2*g3 - b1**2*b5*g2**2*g3 - 2*b3*b4*g0*g2**2*g3 + 2*b2*b5*g0*g2**2*g3 + 3*b4**2*g1*g2**2*g3 - 3*b3*b5*g1*g2**2*g3 - b1**2*b2**2*g3**2 + b1**3*b3*g3**2 + 2*b2**3*g0*g3**2 - 6*b1*b2*b3*g0*g3**2 + 4*b1**2*b4*g0*g3**2 - 3*b3**2*g0**2*g3**2 + 2*b2*b4*g0**2*g3**2 + b1*b5*g0**2*g3**2 - 2*b2**2*b3*g1*g3**2 + 4*b1*b3**2*g1*g3**2 - 2*b1**2*b5*g1*g3**2 + 4*b3*b4*g0*g1*g3**2 - 4*b2*b5*g0*g1*g3**2 - b4**2*g1**2*g3**2 + b3*b5*g1**2*g3**2 + 2*b2*b3**2*g2*g3**2 - 2*b2**2*b4*g2*g3**2 - 2*b1*b3*b4*g2*g3**2 + 2*b1*b2*b5*g2*g3**2 - 2*b4**2*g0*g2*g3**2 + 2*b3*b5*g0*g2*g3**2 - b3**3*g3**3 + 2*b2*b3*b4*g3**3 - b1*b4**2*g3**3 - b2**2*b5*g3**3 + b1*b3*b5*g3**3 - 2*b1**4*g0*g4 - 6*b1**2*b2*g0**2*g4 - 2*b2**2*g0**3*g4 - 4*b1*b3*g0**3*g4 - 2*b4*g0**4*g4 + 2*b1**3*b2*g1*g4 + 8*b1*b2**2*g0*g1*g4 - 2*b1**2*b3*g0*g1*g4 + 8*b2*b3*g0**2*g1*g4 - 2*b1*b4*g0**2*g1*g4 + 2*b5*g0**3*g1*g4 - 2*b2**3*g1**2*g4 + 2*b1**2*b4*g1**2*g4 - 4*b3**2*g0*g1**2*g4 + 4*b1*b5*g0*g1**2*g4 + 2*b3*b4*g1**3*g4 - 2*b2*b5*g1**3*g4 - 2*b1**2*b2**2*g2*g4 + 2*b1**3*b3*g2*g4 - 4*b2**3*g0*g2*g4 + 4*b1*b2*b3*g0*g2*g4 + 2*b3**2*g0**2*g2*g4 - 4*b2*b4*g0**2*g2*g4 + 2*b1*b5*g0**2*g2*g4 + 4*b2**2*b3*g1*g2*g4 - 8*b1*b2*b4*g1*g2*g4 + 4*b1**2*b5*g1*g2*g4 - 2*b4**2*g1**2*g2*g4 + 2*b3*b5*g1**2*g2*g4 - 2*b2*b3**2*g2**2*g4 + 2*b2**2*b4*g2**2*g4 + 2*b1*b3*b4*g2**2*g4 - 2*b1*b2*b5*g2**2*g4 + 2*b4**2*g0*g2**2*g4 - 2*b3*b5*g0*g2**2*g4 + 2*b1*b2**3*g3*g4 - 4*b1**2*b2*b3*g3*g4 + 2*b1**3*b4*g3*g4 + 2*b2**2*b3*g0*g3*g4 - 4*b1*b3**2*g0*g3*g4 + 2*b1**2*b5*g0*g3*g4 - 2*b3*b4*g0**2*g3*g4 + 2*b2*b5*g0**2*g3*g4 - 2*b2*b3**2*g1*g3*g4 + 2*b2**2*b4*g1*g3*g4 + 2*b1*b3*b4*g1*g3*g4 - 2*b1*b2*b5*g1*g3*g4 + 2*b4**2*g0*g1*g3*g4 - 2*b3*b5*g0*g1*g3*g4 + 2*b3**3*g2*g3*g4 - 4*b2*b3*b4*g2*g3*g4 + 2*b1*b4**2*g2*g3*g4 + 2*b2**2*b5*g2*g3*g4 - 2*b1*b3*b5*g2*g3*g4 - b2**4*g4**2 + 3*b1*b2**2*b3*g4**2 - b1**2*b3**2*g4**2 - 2*b1**2*b2*b4*g4**2 + b1**3*b5*g4**2 + 2*b2*b3**2*g0*g4**2 - 2*b2**2*b4*g0*g4**2 - 2*b1*b3*b4*g0*g4**2 + 2*b1*b2*b5*g0*g4**2 - b4**2*g0**2*g4**2 + b3*b5*g0**2*g4**2 - b3**3*g1*g4**2 + 2*b2*b3*b4*g1*g4**2 - b1*b4**2*g1*g4**2 - b2**2*b5*g1*g4**2 + b1*b3*b5*g1*g4**2 - b1**5*g5 - 4*b1**3*b2*g0*g5 - 3*b1*b2**2*g0**2*g5 - 3*b1**2*b3*g0**2*g5 - 2*b2*b3*g0**3*g5 - 2*b1*b4*g0**3*g5 - b5*g0**4*g5 + 3*b1**2*b2**2*g1*g5 - 3*b1**3*b3*g1*g5 + 2*b2**3*g0*g1*g5 + 2*b1*b2*b3*g0*g1*g5 - 4*b1**2*b4*g0*g1*g5 + b3**2*g0**2*g1*g5 + 2*b2*b4*g0**2*g1*g5 - 3*b1*b5*g0**2*g1*g5 - b2**2*b3*g1**2*g5 - 2*b1*b3**2*g1**2*g5 + 4*b1*b2*b4*g1**2*g5 - b1**2*b5*g1**2*g5 - 2*b3*b4*g0*g1**2*g5 + 2*b2*b5*g0*g1**2*g5 + b4**2*g1**3*g5 - b3*b5*g1**3*g5 - 2*b1*b2**3*g2*g5 + 4*b1**2*b2*b3*g2*g5 - 2*b1**3*b4*g2*g5 - 2*b2**2*b3*g0*g2*g5 + 4*b1*b3**2*g0*g2*g5 - 2*b1**2*b5*g0*g2*g5 + 2*b3*b4*g0**2*g2*g5 - 2*b2*b5*g0**2*g2*g5 + 2*b2*b3**2*g1*g2*g5 - 2*b2**2*b4*g1*g2*g5 - 2*b1*b3*b4*g1*g2*g5 + 2*b1*b2*b5*g1*g2*g5 - 2*b4**2*g0*g1*g2*g5 + 2*b3*b5*g0*g1*g2*g5 - b3**3*g2**2*g5 + 2*b2*b3*b4*g2**2*g5 - b1*b4**2*g2**2*g5 - b2**2*b5*g2**2*g5 + b1*b3*b5*g2**2*g5 + b2**4*g3*g5 - 3*b1*b2**2*b3*g3*g5 + b1**2*b3**2*g3*g5 + 2*b1**2*b2*b4*g3*g5 - b1**3*b5*g3*g5 - 2*b2*b3**2*g0*g3*g5 + 2*b2**2*b4*g0*g3*g5 + 2*b1*b3*b4*g0*g3*g5 - 2*b1*b2*b5*g0*g3*g5 + b4**2*g0**2*g3*g5 - b3*b5*g0**2*g3*g5 + b3**3*g1*g3*g5 - 2*b2*b3*b4*g1*g3*g5 + b1*b4**2*g1*g3*g5 + b2**2*b5*g1*g3*g5 - b1*b3*b5*g1*g3*g5) q1 <- ( -(b1*g0**5) - 4*b1**2*g0**3*g1 + b2*g0**4*g1 - 3*b1**3*g0*g1**2 + 6*b1*b2*g0**2*g1**2 - b3*g0**3*g1**2 + 3*b1**2*b2*g1**3 - 2*b2**2*g0*g1**3 - 4*b1*b3*g0*g1**3 + b4*g0**2*g1**3 + 2*b2*b3*g1**4 + 2*b1*b4*g1**4 - b5*g0*g1**4 + b6*g1**5 - 3*b1**3*g0**2*g2 - 2*b1*b2*g0**3*g2 + b3*g0**4*g2 - 2*b1**4*g1*g2 + 6*b1**2*b2*g0*g1*g2 + b2**2*g0**2*g1*g2 + 5*b1*b3*g0**2*g1*g2 - 2*b4*g0**3*g1*g2 - 5*b1*b2**2*g1**2*g2 - 4*b1**2*b3*g1**2*g2 - 2*b2*b3*g0*g1**2*g2 - 4*b1*b4*g0*g1**2*g2 + 3*b5*g0**2*g1**2*g2 - b3**2*g1**3*g2 - 2*b2*b4*g1**3*g2 - b1*b5*g1**3*g2 - 4*b6*g0*g1**3*g2 + 2*b1**3*b2*g2**2 - 2*b1*b2**2*g0*g2**2 + 5*b1**2*b3*g0*g2**2 + b2*b3*g0**2*g2**2 - b1*b4*g0**2*g2**2 - b5*g0**3*g2**2 + 2*b2**3*g1*g2**2 + 4*b1*b2*b3*g1*g2**2 - 3*b1**2*b4*g1*g2**2 + 2*b2*b4*g0*g1*g2**2 + 4*b1*b5*g0*g1*g2**2 + 3*b6*g0**2*g1*g2**2 + 2*b3*b4*g1**2*g2**2 + b2*b5*g1**2*g2**2 - 3*b1*b6*g1**2*g2**2 - b2**2*b3*g2**3 - 2*b1*b3**2*g2**3 + 2*b1*b2*b4*g2**3 + b1**2*b5*g2**3 - 2*b2*b5*g0*g2**3 + 2*b1*b6*g0*g2**3 - b4**2*g1*g2**3 - b3*b5*g1*g2**3 + 2*b2*b6*g1*g2**3 + b4*b5*g2**4 - b3*b6*g2**4 - 2*b1**4*g0*g3 - 3*b1**2*b2*g0**2*g3 + b2**2*g0**3*g3 - b1*b3*g0**3*g3 + b4*g0**4*g3 + 4*b1**3*b2*g1*g3 + 6*b1*b2**2*g0*g1*g3 - 3*b2*b3*g0**2*g1*g3 + 3*b1*b4*g0**2*g1*g3 - 2*b5*g0**3*g1*g3 - b2**3*g1**2*g3 + 4*b1**2*b4*g1**2*g3 + 3*b3**2*g0*g1**2*g3 + 4*b2*b4*g0*g1**2*g3 - b1*b5*g0*g1**2*g3 + 3*b6*g0**2*g1**2*g3 - 2*b3*b4*g1**3*g3 - b2*b5*g1**3*g3 + 3*b1*b6*g1**3*g3 - 4*b1**2*b2**2*g2*g3 + 2*b1**3*b3*g2*g3 - 6*b1*b2*b3*g0*g2*g3 - 2*b2*b4*g0**2*g2*g3 - 4*b1*b5*g0**2*g2*g3 - 2*b6*g0**3*g2*g3 - 2*b2**2*b3*g1*g2*g3 + 4*b1*b3**2*g1*g2*g3 - 2*b1*b2*b4*g1*g2*g3 - 4*b3*b4*g0*g1*g2*g3 + 2*b2*b5*g0*g1*g2*g3 + 2*b1*b6*g0*g1*g2*g3 + 2*b4**2*g1**2*g2*g3 + 2*b3*b5*g1**2*g2*g3 - 4*b2*b6*g1**2*g2*g3 + 2*b2*b3**2*g2**2*g3 - b2**2*b4*g2**2*g3 - 2*b1*b2*b5*g2**2*g3 + b1**2*b6*g2**2*g3 + b4**2*g0*g2**2*g3 + b3*b5*g0*g2**2*g3 - 2*b2*b6*g0*g2**2*g3 - 3*b4*b5*g1*g2**2*g3 + 3*b3*b6*g1*g2**2*g3 + 2*b1*b2**3*g3**2 - 3*b1**2*b2*b3*g3**2 + b1**3*b4*g3**2 + 2*b1*b2*b4*g0*g3**2 - 2*b1**2*b5*g0*g3**2 + 2*b3*b4*g0**2*g3**2 - b2*b5*g0**2*g3**2 - b1*b6*g0**2*g3**2 + 2*b2**2*b4*g1*g3**2 - 4*b1*b3*b4*g1*g3**2 + 2*b1**2*b6*g1*g3**2 - 2*b4**2*g0*g1*g3**2 - 2*b3*b5*g0*g1*g3**2 + 4*b2*b6*g0*g1*g3**2 + b4*b5*g1**2*g3**2 - b3*b6*g1**2*g3**2 - b3**3*g2*g3**2 + b1*b4**2*g2*g3**2 + b2**2*b5*g2*g3**2 + b1*b3*b5*g2*g3**2 - 2*b1*b2*b6*g2*g3**2 + 2*b4*b5*g0*g2*g3**2 - 2*b3*b6*g0*g2*g3**2 + b3**2*b4*g3**3 - b2*b4**2*g3**3 - b2*b3*b5*g3**3 + b1*b4*b5*g3**3 + b2**2*b6*g3**3 - b1*b3*b6*g3**3 - b1**5*g4 - 2*b1**3*b2*g0*g4 + b1*b2**2*g0**2*g4 - b1**2*b3*g0**2*g4 + 2*b2*b3*g0**3*g4 + b5*g0**4*g4 + b1**2*b2**2*g1*g4 - 3*b1**3*b3*g1*g4 - 2*b2**3*g0*g1*g4 - 2*b1*b2*b3*g0*g1*g4 - 2*b1**2*b4*g0*g1*g4 - 3*b3**2*g0**2*g1*g4 - 4*b2*b4*g0**2*g1*g4 + b1*b5*g0**2*g1*g4 - 2*b6*g0**3*g1*g4 + 3*b2**2*b3*g1**2*g4 - 2*b1*b3**2*g1**2*g4 - b1**2*b5*g1**2*g4 + 4*b3*b4*g0*g1**2*g4 - 4*b1*b6*g0*g1**2*g4 - b4**2*g1**3*g4 - b3*b5*g1**3*g4 + 2*b2*b6*g1**3*g4 + 2*b1**2*b2*b3*g2*g4 - 2*b1**3*b4*g2*g4 + 4*b2**2*b3*g0*g2*g4 - 4*b1*b2*b4*g0*g2*g4 + 2*b2*b5*g0**2*g2*g4 - 2*b1*b6*g0**2*g2*g4 - 4*b2*b3**2*g1*g2*g4 + 4*b1*b3*b4*g1*g2*g4 + 4*b1*b2*b5*g1*g2*g4 - 4*b1**2*b6*g1*g2*g4 + 2*b4*b5*g1**2*g2*g4 - 2*b3*b6*g1**2*g2*g4 + b3**3*g2**2*g4 - b1*b4**2*g2**2*g4 - b2**2*b5*g2**2*g4 - b1*b3*b5*g2**2*g4 + 2*b1*b2*b6*g2**2*g4 - 2*b4*b5*g0*g2**2*g4 + 2*b3*b6*g0*g2**2*g4 - b2**4*g3*g4 + b1*b2**2*b3*g3*g4 + b1**2*b3**2*g3*g4 - b1**3*b5*g3*g4 - 2*b2**2*b4*g0*g3*g4 + 4*b1*b3*b4*g0*g3*g4 - 2*b1**2*b6*g0*g3*g4 + b4**2*g0**2*g3*g4 + b3*b5*g0**2*g3*g4 - 2*b2*b6*g0**2*g3*g4 + b3**3*g1*g3*g4 - b1*b4**2*g1*g3*g4 - b2**2*b5*g1*g3*g4 - b1*b3*b5*g1*g3*g4 + 2*b1*b2*b6*g1*g3*g4 - 2*b4*b5*g0*g1*g3*g4 + 2*b3*b6*g0*g1*g3*g4 - 2*b3**2*b4*g2*g3*g4 + 2*b2*b4**2*g2*g3*g4 + 2*b2*b3*b5*g2*g3*g4 - 2*b1*b4*b5*g2*g3*g4 - 2*b2**2*b6*g2*g3*g4 + 2*b1*b3*b6*g2*g3*g4 + b2**3*b3*g4**2 - 2*b1*b2*b3**2*g4**2 - b1*b2**2*b4*g4**2 + 2*b1**2*b3*b4*g4**2 + b1**2*b2*b5*g4**2 - b1**3*b6*g4**2 - b3**3*g0*g4**2 + b1*b4**2*g0*g4**2 + b2**2*b5*g0*g4**2 + b1*b3*b5*g0*g4**2 - 2*b1*b2*b6*g0*g4**2 + b4*b5*g0**2*g4**2 - b3*b6*g0**2*g4**2 + b3**2*b4*g1*g4**2 - b2*b4**2*g1*g4**2 - b2*b3*b5*g1*g4**2 + b1*b4*b5*g1*g4**2 + b2**2*b6*g1*g4**2 - b1*b3*b6*g1*g4**2 + b1**4*b2*g5 + 3*b1**2*b2**2*g0*g5 + b1**3*b3*g0*g5 + b2**3*g0**2*g5 + 4*b1*b2*b3*g0**2*g5 + b1**2*b4*g0**2*g5 + b3**2*g0**3*g5 + 2*b2*b4*g0**3*g5 + b1*b5*g0**3*g5 + b6*g0**4*g5 - 2*b1*b2**3*g1*g5 + b1**2*b2*b3*g1*g5 + b1**3*b4*g1*g5 - 4*b2**2*b3*g0*g1*g5 + 2*b1*b2*b4*g0*g1*g5 + 2*b1**2*b5*g0*g1*g5 - 2*b3*b4*g0**2*g1*g5 - b2*b5*g0**2*g1*g5 + 3*b1*b6*g0**2*g1*g5 + 2*b2*b3**2*g1**2*g5 - b2**2*b4*g1**2*g5 - 2*b1*b2*b5*g1**2*g5 + b1**2*b6*g1**2*g5 + b4**2*g0*g1**2*g5 + b3*b5*g0*g1**2*g5 - 2*b2*b6*g0*g1**2*g5 - b4*b5*g1**3*g5 + b3*b6*g1**3*g5 + b2**4*g2*g5 - b1*b2**2*b3*g2*g5 - b1**2*b3**2*g2*g5 + b1**3*b5*g2*g5 + 2*b2**2*b4*g0*g2*g5 - 4*b1*b3*b4*g0*g2*g5 + 2*b1**2*b6*g0*g2*g5 - b4**2*g0**2*g2*g5 - b3*b5*g0**2*g2*g5 + 2*b2*b6*g0**2*g2*g5 - b3**3*g1*g2*g5 + b1*b4**2*g1*g2*g5 + b2**2*b5*g1*g2*g5 + b1*b3*b5*g1*g2*g5 - 2*b1*b2*b6*g1*g2*g5 + 2*b4*b5*g0*g1*g2*g5 - 2*b3*b6*g0*g1*g2*g5 + b3**2*b4*g2**2*g5 - b2*b4**2*g2**2*g5 - b2*b3*b5*g2**2*g5 + b1*b4*b5*g2**2*g5 + b2**2*b6*g2**2*g5 - b1*b3*b6*g2**2*g5 - b2**3*b3*g3*g5 + 2*b1*b2*b3**2*g3*g5 + b1*b2**2*b4*g3*g5 - 2*b1**2*b3*b4*g3*g5 - b1**2*b2*b5*g3*g5 + b1**3*b6*g3*g5 + b3**3*g0*g3*g5 - b1*b4**2*g0*g3*g5 - b2**2*b5*g0*g3*g5 - b1*b3*b5*g0*g3*g5 + 2*b1*b2*b6*g0*g3*g5 - b4*b5*g0**2*g3*g5 + b3*b6*g0**2*g3*g5 - b3**2*b4*g1*g3*g5 + b2*b4**2*g1*g3*g5 + b2*b3*b5*g1*g3*g5 - b1*b4*b5*g1*g3*g5 - b2**2*b6*g1*g3*g5 + b1*b3*b6*g1*g3*g5)/den q2 <- (-(b1**2*g0**4) - b2*g0**5 - 3*b1**3*g0**2*g1 - 2*b1*b2*g0**3*g1 + b3*g0**4*g1 - b1**4*g1**2 + 3*b1**2*b2*g0*g1**2 + 2*b2**2*g0**2*g1**2 + b1*b3*g0**2*g1**2 - b4*g0**3*g1**2 - b1*b2**2*g1**3 - 2*b1**2*b3*g1**3 - 2*b2*b3*g0*g1**3 + b5*g0**2*g1**3 - b1*b5*g1**4 - b6*g0*g1**4 - 2*b1**4*g0*g2 - 3*b1**2*b2*g0**2*g2 - 3*b2**2*g0**3*g2 + 3*b1*b3*g0**3*g2 + b4*g0**4*g2 + 4*b1**3*b2*g1*g2 + 6*b1**2*b3*g0*g1*g2 + 3*b2*b3*g0**2*g1*g2 - 3*b1*b4*g0**2*g1*g2 - 2*b5*g0**3*g1*g2 + b2**3*g1**2*g2 + 2*b1*b2*b3*g1**2*g2 + b3**2*g0*g1**2*g2 + 5*b1*b5*g0*g1**2*g2 + 3*b6*g0**2*g1**2*g2 + b2*b5*g1**3*g2 - b1*b6*g1**3*g2 - 3*b1**2*b2**2*g2**2 + 2*b1**3*b3*g2**2 - 2*b2**3*g0*g2**2 + 2*b1*b2*b3*g0*g2**2 - 3*b1**2*b4*g0*g2**2 - 2*b3**2*g0**2*g2**2 + 3*b2*b4*g0**2*g2**2 - 4*b1*b5*g0**2*g2**2 - b6*g0**3*g2**2 - b2**2*b3*g1*g2**2 - 2*b1*b3**2*g1*g2**2 + 3*b1**2*b5*g1*g2**2 - 4*b2*b5*g0*g1*g2**2 + 4*b1*b6*g0*g1*g2**2 - b3*b5*g1**2*g2**2 + b2*b6*g1**2*g2**2 + b2**2*b4*g2**3 + 2*b1*b3*b4*g2**3 - 4*b1*b2*b5*g2**3 + b1**2*b6*g2**3 - b4**2*g0*g2**3 + 3*b3*b5*g0*g2**3 - 2*b2*b6*g0*g2**3 + b4*b5*g1*g2**3 - b3*b6*g1*g2**3 - b5**2*g2**4 + b4*b6*g2**4 - b1**5*g3 - 2*b1**3*b2*g0*g3 - 2*b1*b2**2*g0**2*g3 + 2*b1**2*b3*g0**2*g3 - b2*b3*g0**3*g3 + 3*b1*b4*g0**3*g3 + b5*g0**4*g3 - b1**2*b2**2*g1*g3 - b1**3*b3*g1*g3 + 2*b2**3*g0*g1*g3 - 8*b1*b2*b3*g0*g1*g3 - 3*b2*b4*g0**2*g1*g3 - 3*b1*b5*g0**2*g1*g3 - 2*b6*g0**3*g1*g3 - b2**2*b3*g1**2*g3 + 3*b1*b3**2*g1**2*g3 - 2*b1**2*b5*g1**2*g3 + b2*b5*g0*g1**2*g3 - b1*b6*g0*g1**2*g3 + b3*b5*g1**3*g3 - b2*b6*g1**3*g3 + 2*b1*b2**3*g2*g3 - 2*b1**3*b4*g2*g3 + 2*b2**2*b3*g0*g2*g3 - 2*b1*b3**2*g0*g2*g3 + 6*b1*b2*b4*g0*g2*g3 - 6*b1**2*b5*g0*g2*g3 + 4*b2*b5*g0**2*g2*g3 - 4*b1*b6*g0**2*g2*g3 + 2*b2*b3**2*g1*g2*g3 - 2*b2**2*b4*g1*g2*g3 - 4*b1*b3*b4*g1*g2*g3 + 4*b1*b2*b5*g1*g2*g3 + 2*b4**2*g0*g1*g2*g3 - 4*b3*b5*g0*g1*g2*g3 + 2*b2*b6*g0*g1*g2*g3 - 2*b4*b5*g1**2*g2*g3 + 2*b3*b6*g1**2*g2*g3 - 2*b2*b3*b4*g2**2*g3 - b1*b4**2*g2**2*g3 + 2*b2**2*b5*g2**2*g3 + 3*b1*b3*b5*g2**2*g3 - 2*b1*b2*b6*g2**2*g3 - b4*b5*g0*g2**2*g3 + b3*b6*g0*g2**2*g3 + 3*b5**2*g1*g2**2*g3 - 3*b4*b6*g1*g2**2*g3 - 2*b1*b2**2*b3*g3**2 + b1**2*b3**2*g3**2 + 3*b1**2*b2*b4*g3**2 - 2*b1**3*b5*g3**2 - 2*b2**2*b4*g0*g3**2 + 2*b1*b3*b4*g0*g3**2 + 2*b1*b2*b5*g0*g3**2 - 2*b1**2*b6*g0*g3**2 - b4**2*g0**2*g3**2 + 2*b3*b5*g0**2*g3**2 - b2*b6*g0**2*g3**2 - b3**3*g1*g3**2 + 2*b2*b3*b4*g1*g3**2 + b1*b4**2*g1*g3**2 - b2**2*b5*g1*g3**2 - b1*b3*b5*g1*g3**2 + 2*b4*b5*g0*g1*g3**2 - 2*b3*b6*g0*g1*g3**2 - b5**2*g1**2*g3**2 + b4*b6*g1**2*g3**2 + b3**2*b4*g2*g3**2 + b2*b4**2*g2*g3**2 - 3*b2*b3*b5*g2*g3**2 - b1*b4*b5*g2*g3**2 + b2**2*b6*g2*g3**2 + b1*b3*b6*g2*g3**2 - 2*b5**2*g0*g2*g3**2 + 2*b4*b6*g0*g2*g3**2 - b3*b4**2*g3**3 + b3**2*b5*g3**3 + b2*b4*b5*g3**3 - b1*b5**2*g3**3 - b2*b3*b6*g3**3 + b1*b4*b6*g3**3 + b1**4*b2*g4 + b1**2*b2**2*g0*g4 + 3*b1**3*b3*g0*g4 - b2**3*g0**2*g4 + 4*b1*b2*b3*g0**2*g4 + 3*b1**2*b4*g0**2*g4 + b3**2*g0**3*g4 + 3*b1*b5*g0**3*g4 + b6*g0**4*g4 - b1**2*b2*b3*g1*g4 + b1**3*b4*g1*g4 + 2*b2**2*b3*g0*g1*g4 - 6*b1*b2*b4*g0*g1*g4 + 4*b1**2*b5*g0*g1*g4 - b2*b5*g0**2*g1*g4 + b1*b6*g0**2*g1*g4 - 2*b2*b3**2*g1**2*g4 + b2**2*b4*g1**2*g4 + 2*b1*b3*b4*g1**2*g4 - b1**2*b6*g1**2*g4 - b4**2*g0*g1**2*g4 + b3*b5*g0*g1**2*g4 + b4*b5*g1**3*g4 - b3*b6*g1**3*g4 - b2**4*g2*g4 + 3*b1*b2**2*b3*g2*g4 - 3*b1**2*b3**2*g2*g4 + b1**3*b5*g2*g4 - 2*b2*b3**2*g0*g2*g4 + 2*b2**2*b4*g0*g2*g4 - 2*b1*b3*b4*g0*g2*g4 + 2*b1*b2*b5*g0*g2*g4 + b4**2*g0**2*g2*g4 - 3*b3*b5*g0**2*g2*g4 + 2*b2*b6*g0**2*g2*g4 + b3**3*g1*g2*g4 + 2*b2*b3*b4*g1*g2*g4 + b1*b4**2*g1*g2*g4 - 3*b2**2*b5*g1*g2*g4 - 5*b1*b3*b5*g1*g2*g4 + 4*b1*b2*b6*g1*g2*g4 - 2*b5**2*g1**2*g2*g4 + 2*b4*b6*g1**2*g2*g4 - b3**2*b4*g2**2*g4 - b2*b4**2*g2**2*g4 + 3*b2*b3*b5*g2**2*g4 + b1*b4*b5*g2**2*g4 - b2**2*b6*g2**2*g4 - b1*b3*b6*g2**2*g4 + 2*b5**2*g0*g2**2*g4 - 2*b4*b6*g0*g2**2*g4 + b2**3*b3*g3*g4 - 3*b1*b2**2*b4*g3*g4 + 3*b1**2*b2*b5*g3*g4 - b1**3*b6*g3*g4 + b3**3*g0*g3*g4 - 2*b2*b3*b4*g0*g3*g4 - b1*b4**2*g0*g3*g4 + b2**2*b5*g0*g3*g4 + b1*b3*b5*g0*g3*g4 - b4*b5*g0**2*g3*g4 + b3*b6*g0**2*g3*g4 - b3**2*b4*g1*g3*g4 - b2*b4**2*g1*g3*g4 + 3*b2*b3*b5*g1*g3*g4 + b1*b4*b5*g1*g3*g4 - b2**2*b6*g1*g3*g4 - b1*b3*b6*g1*g3*g4 + 2*b5**2*g0*g1*g3*g4 - 2*b4*b6*g0*g1*g3*g4 + 2*b3*b4**2*g2*g3*g4 - 2*b3**2*b5*g2*g3*g4 - 2*b2*b4*b5*g2*g3*g4 + 2*b1*b5**2*g2*g3*g4 + 2*b2*b3*b6*g2*g3*g4 - 2*b1*b4*b6*g2*g3*g4 - b2**2*b3**2*g4**2 + b1*b3**3*g4**2 + b2**3*b4*g4**2 - b1*b2**2*b5*g4**2 - b1**2*b3*b5*g4**2 + b1**2*b2*b6*g4**2 + b3**2*b4*g0*g4**2 + b2*b4**2*g0*g4**2 - 3*b2*b3*b5*g0*g4**2 - b1*b4*b5*g0*g4**2 + b2**2*b6*g0*g4**2 + b1*b3*b6*g0*g4**2 - b5**2*g0**2*g4**2 + b4*b6*g0**2*g4**2 - b3*b4**2*g1*g4**2 + b3**2*b5*g1*g4**2 + b2*b4*b5*g1*g4**2 - b1*b5**2*g1*g4**2 - b2*b3*b6*g1*g4**2 + b1*b4*b6*g1*g4**2 - b1**3*b2**2*g5 + b1**4*b3*g5 - 2*b1*b2**3*g0*g5 + b1**2*b2*b3*g0*g5 + b1**3*b4*g0*g5 - 2*b2**2*b3*g0**2*g5 + b1*b3**2*g0**2*g5 + b1**2*b5*g0**2*g5 - b2*b5*g0**3*g5 + b1*b6*g0**3*g5 + b2**4*g1*g5 - b1*b2**2*b3*g1*g5 + 2*b1**2*b3**2*g1*g5 - 3*b1**2*b2*b4*g1*g5 + b1**3*b5*g1*g5 + 2*b2*b3**2*g0*g1*g5 - 4*b1*b2*b5*g0*g1*g5 + 2*b1**2*b6*g0*g1*g5 + b3*b5*g0**2*g1*g5 - b2*b6*g0**2*g1*g5 - 2*b2*b3*b4*g1**2*g5 - b1*b4**2*g1**2*g5 + 2*b2**2*b5*g1**2*g5 + 3*b1*b3*b5*g1**2*g5 - 2*b1*b2*b6*g1**2*g5 - b4*b5*g0*g1**2*g5 + b3*b6*g0*g1**2*g5 + b5**2*g1**3*g5 - b4*b6*g1**3*g5 - b2**3*b3*g2*g5 + 3*b1*b2**2*b4*g2*g5 - 3*b1**2*b2*b5*g2*g5 + b1**3*b6*g2*g5 - b3**3*g0*g2*g5 + 2*b2*b3*b4*g0*g2*g5 + b1*b4**2*g0*g2*g5 - b2**2*b5*g0*g2*g5 - b1*b3*b5*g0*g2*g5 + b4*b5*g0**2*g2*g5 - b3*b6*g0**2*g2*g5 + b3**2*b4*g1*g2*g5 + b2*b4**2*g1*g2*g5 - 3*b2*b3*b5*g1*g2*g5 - b1*b4*b5*g1*g2*g5 + b2**2*b6*g1*g2*g5 + b1*b3*b6*g1*g2*g5 - 2*b5**2*g0*g1*g2*g5 + 2*b4*b6*g0*g1*g2*g5 - b3*b4**2*g2**2*g5 + b3**2*b5*g2**2*g5 + b2*b4*b5*g2**2*g5 - b1*b5**2*g2**2*g5 - b2*b3*b6*g2**2*g5 + b1*b4*b6*g2**2*g5 + b2**2*b3**2*g3*g5 - b1*b3**3*g3*g5 - b2**3*b4*g3*g5 + b1*b2**2*b5*g3*g5 + b1**2*b3*b5*g3*g5 - b1**2*b2*b6*g3*g5 - b3**2*b4*g0*g3*g5 - b2*b4**2*g0*g3*g5 + 3*b2*b3*b5*g0*g3*g5 + b1*b4*b5*g0*g3*g5 - b2**2*b6*g0*g3*g5 - b1*b3*b6*g0*g3*g5 + b5**2*g0**2*g3*g5 - b4*b6*g0**2*g3*g5 + b3*b4**2*g1*g3*g5 - b3**2*b5*g1*g3*g5 - b2*b4*b5*g1*g3*g5 + b1*b5**2*g1*g3*g5 + b2*b3*b6*g1*g3*g5 - b1*b4*b6*g1*g3*g5)/den q3 <- ( -(b1**3*g0**3) - 2*b1*b2*g0**4 - b3*g0**5 - 2*b1**4*g0*g1 - 3*b1**2*b2*g0**2*g1 + 2*b2**2*g0**3*g1 - 2*b1*b3*g0**3*g1 + b4*g0**4*g1 + 2*b1**3*b2*g1**2 + 5*b1*b2**2*g0*g1**2 - 2*b1**2*b3*g0*g1**2 - b2*b3*g0**2*g1**2 + b1*b4*g0**2*g1**2 - b5*g0**3*g1**2 - b2**3*g1**3 + 2*b1**2*b4*g1**3 + 2*b2*b4*g0*g1**3 + b6*g0**2*g1**3 - b2*b5*g1**4 + b1*b6*g1**4 - b1**5*g2 - 2*b1**3*b2*g0*g2 - 2*b1*b2**2*g0**2*g2 + 2*b1**2*b3*g0**2*g2 - b2*b3*g0**3*g2 + 3*b1*b4*g0**3*g2 + b5*g0**4*g2 - b1**2*b2**2*g1*g2 - b1**3*b3*g1*g2 - 4*b1*b2*b3*g0*g1*g2 - 2*b1**2*b4*g0*g1*g2 + 2*b3**2*g0**2*g1*g2 - 5*b2*b4*g0**2*g1*g2 - 3*b1*b5*g0**2*g1*g2 - 2*b6*g0**3*g1*g2 + b2**2*b3*g1**2*g2 + b1*b3**2*g1**2*g2 - 2*b1*b2*b4*g1**2*g2 - 2*b3*b4*g0*g1**2*g2 + 3*b2*b5*g0*g1**2*g2 - b1*b6*g0*g1**2*g2 + b3*b5*g1**3*g2 - b2*b6*g1**3*g2 + 2*b1**2*b2*b3*g2**2 - 2*b1**3*b4*g2**2 + b2**2*b3*g0*g2**2 + 2*b1*b2*b4*g0*g2**2 - 3*b1**2*b5*g0*g2**2 + 2*b2*b5*g0**2*g2**2 - 2*b1*b6*g0**2*g2**2 - b2**2*b4*g1*g2**2 + 2*b1*b2*b5*g1*g2**2 - b1**2*b6*g1*g2**2 + 2*b4**2*g0*g1*g2**2 - 2*b3*b5*g0*g1*g2**2 - b4*b5*g1**2*g2**2 + b3*b6*g1**2*g2**2 - b1*b4**2*g2**3 + b1*b3*b5*g2**3 - b4*b5*g0*g2**3 + b3*b6*g0*g2**3 + b5**2*g1*g2**3 - b4*b6*g1*g2**3 + b1**4*b2*g3 + b1**2*b2**2*g0*g3 + 3*b1**3*b3*g0*g3 + 2*b2**3*g0**2*g3 - 2*b1*b2*b3*g0**2*g3 + 6*b1**2*b4*g0**2*g3 - 2*b3**2*g0**3*g3 + 3*b2*b4*g0**3*g3 + 3*b1*b5*g0**3*g3 + b6*g0**4*g3 + 2*b1*b2**3*g1*g3 - 5*b1**2*b2*b3*g1*g3 + 3*b1**3*b4*g1*g3 - 2*b2**2*b3*g0*g1*g3 + 2*b1*b3**2*g0*g1*g3 + 4*b3*b4*g0**2*g1*g3 - 5*b2*b5*g0**2*g1*g3 + b1*b6*g0**2*g1*g3 - b2*b3**2*g1**2*g3 + 3*b2**2*b4*g1**2*g3 - 2*b1*b3*b4*g1**2*g3 - 2*b1*b2*b5*g1**2*g3 + 2*b1**2*b6*g1**2*g3 - 2*b4**2*g0*g1**2*g3 - b3*b5*g0*g1**2*g3 + 3*b2*b6*g0*g1**2*g3 + b4*b5*g1**3*g3 - b3*b6*g1**3*g3 - 2*b1*b2**2*b3*g2*g3 + 4*b1**2*b2*b4*g2*g3 - 2*b1**3*b5*g2*g3 - 2*b2**2*b4*g0*g2*g3 + 4*b1*b2*b5*g0*g2*g3 - 2*b1**2*b6*g0*g2*g3 - 2*b4**2*g0**2*g2*g3 + 2*b3*b5*g0**2*g2*g3 + 4*b1*b4**2*g1*g2*g3 - 4*b1*b3*b5*g1*g2*g3 + 4*b4*b5*g0*g1*g2*g3 - 4*b3*b6*g0*g1*g2*g3 - 2*b5**2*g1**2*g2*g3 + 2*b4*b6*g1**2*g2*g3 + b2*b4**2*g2**2*g3 - b2*b3*b5*g2**2*g3 - b1*b4*b5*g2**2*g3 + b1*b3*b6*g2**2*g3 - b5**2*g0*g2**2*g3 + b4*b6*g0*g2**2*g3 + 2*b1*b2*b3**2*g3**2 - 2*b1*b2**2*b4*g3**2 - 2*b1**2*b3*b4*g3**2 + 2*b1**2*b2*b5*g3**2 - b3**3*g0*g3**2 + 4*b2*b3*b4*g0*g3**2 - 5*b1*b4**2*g0*g3**2 - 3*b2**2*b5*g0*g3**2 + 3*b1*b3*b5*g0*g3**2 + 2*b1*b2*b6*g0*g3**2 - 2*b4*b5*g0**2*g3**2 + 2*b3*b6*g0**2*g3**2 + b3**2*b4*g1*g3**2 - 3*b2*b4**2*g1*g3**2 + b2*b3*b5*g1*g3**2 + 3*b1*b4*b5*g1*g3**2 + b2**2*b6*g1*g3**2 - 3*b1*b3*b6*g1*g3**2 + 2*b5**2*g0*g1*g3**2 - 2*b4*b6*g0*g1*g3**2 - b3*b4**2*g2*g3**2 + b3**2*b5*g2*g3**2 + b2*b4*b5*g2*g3**2 - b1*b5**2*g2*g3**2 - b2*b3*b6*g2*g3**2 + b1*b4*b6*g2*g3**2 + b4**3*g3**3 - 2*b3*b4*b5*g3**3 + b2*b5**2*g3**3 + b3**2*b6*g3**3 - b2*b4*b6*g3**3 - b1**3*b2**2*g4 + b1**4*b3*g4 - 3*b1**2*b2*b3*g0*g4 + 3*b1**3*b4*g0*g4 - 3*b1*b3**2*g0**2*g4 + 3*b1**2*b5*g0**2*g4 - 2*b3*b4*g0**3*g4 + b2*b5*g0**3*g4 + b1*b6*g0**3*g4 - b2**4*g1*g4 + 3*b1*b2**2*b3*g1*g4 + 2*b1**2*b3**2*g1*g4 - 5*b1**2*b2*b4*g1*g4 + b1**3*b5*g1*g4 + 2*b2*b3**2*g0*g1*g4 - 2*b2**2*b4*g0*g1*g4 + 4*b1*b3*b4*g0*g1*g4 - 4*b1*b2*b5*g0*g1*g4 + 2*b4**2*g0**2*g1*g4 + b3*b5*g0**2*g1*g4 - 3*b2*b6*g0**2*g1*g4 - 3*b1*b4**2*g1**2*g4 + 3*b1*b3*b5*g1**2*g4 - 3*b4*b5*g0*g1**2*g4 + 3*b3*b6*g0*g1**2*g4 + b5**2*g1**3*g4 - b4*b6*g1**3*g4 + b2**3*b3*g2*g4 - 2*b1*b2*b3**2*g2*g4 - b1*b2**2*b4*g2*g4 + 2*b1**2*b3*b4*g2*g4 + b1**2*b2*b5*g2*g4 - b1**3*b6*g2*g4 + b3**3*g0*g2*g4 - 4*b2*b3*b4*g0*g2*g4 + 3*b1*b4**2*g0*g2*g4 + 3*b2**2*b5*g0*g2*g4 - b1*b3*b5*g0*g2*g4 - 2*b1*b2*b6*g0*g2*g4 + b4*b5*g0**2*g2*g4 - b3*b6*g0**2*g2*g4 - b3**2*b4*g1*g2*g4 + b2*b4**2*g1*g2*g4 + b2*b3*b5*g1*g2*g4 - b1*b4*b5*g1*g2*g4 - b2**2*b6*g1*g2*g4 + b1*b3*b6*g1*g2*g4 + b3*b4**2*g2**2*g4 - b3**2*b5*g2**2*g4 - b2*b4*b5*g2**2*g4 + b1*b5**2*g2**2*g4 + b2*b3*b6*g2**2*g4 - b1*b4*b6*g2**2*g4 - b2**2*b3**2*g3*g4 - b1*b3**3*g3*g4 + b2**3*b4*g3*g4 + 4*b1*b2*b3*b4*g3*g4 - 2*b1**2*b4**2*g3*g4 - 3*b1*b2**2*b5*g3*g4 + b1**2*b3*b5*g3*g4 + b1**2*b2*b6*g3*g4 - b3**2*b4*g0*g3*g4 + 3*b2*b4**2*g0*g3*g4 - b2*b3*b5*g0*g3*g4 - 3*b1*b4*b5*g0*g3*g4 - b2**2*b6*g0*g3*g4 + 3*b1*b3*b6*g0*g3*g4 - b5**2*g0**2*g3*g4 + b4*b6*g0**2*g3*g4 + b3*b4**2*g1*g3*g4 - b3**2*b5*g1*g3*g4 - b2*b4*b5*g1*g3*g4 + b1*b5**2*g1*g3*g4 + b2*b3*b6*g1*g3*g4 - b1*b4*b6*g1*g3*g4 - 2*b4**3*g2*g3*g4 + 4*b3*b4*b5*g2*g3*g4 - 2*b2*b5**2*g2*g3*g4 - 2*b3**2*b6*g2*g3*g4 + 2*b2*b4*b6*g2*g3*g4 + b2*b3**3*g4**2 - 2*b2**2*b3*b4*g4**2 - b1*b3**2*b4*g4**2 + 2*b1*b2*b4**2*g4**2 + b2**3*b5*g4**2 - b1**2*b4*b5*g4**2 - b1*b2**2*b6*g4**2 + b1**2*b3*b6*g4**2 - b3*b4**2*g0*g4**2 + b3**2*b5*g0*g4**2 + b2*b4*b5*g0*g4**2 - b1*b5**2*g0*g4**2 - b2*b3*b6*g0*g4**2 + b1*b4*b6*g0*g4**2 + b4**3*g1*g4**2 - 2*b3*b4*b5*g1*g4**2 + b2*b5**2*g1*g4**2 + b3**2*b6*g1*g4**2 - b2*b4*b6*g1*g4**2 + b1**2*b2**3*g5 - 2*b1**3*b2*b3*g5 + b1**4*b4*g5 + b2**4*g0*g5 - b1*b2**2*b3*g0*g5 - 2*b1**2*b3**2*g0*g5 + b1**2*b2*b4*g0*g5 + b1**3*b5*g0*g5 - b2*b3**2*g0**2*g5 + 2*b2**2*b4*g0**2*g5 - 2*b1*b3*b4*g0**2*g5 + b1**2*b6*g0**2*g5 - b3*b5*g0**3*g5 + b2*b6*g0**3*g5 - b2**3*b3*g1*g5 + 3*b1*b2**2*b4*g1*g5 - 3*b1**2*b2*b5*g1*g5 + b1**3*b6*g1*g5 + 2*b1*b4**2*g0*g1*g5 - 2*b1*b3*b5*g0*g1*g5 + b4*b5*g0**2*g1*g5 - b3*b6*g0**2*g1*g5 + b2*b4**2*g1**2*g5 - b2*b3*b5*g1**2*g5 - b1*b4*b5*g1**2*g5 + b1*b3*b6*g1**2*g5 - b5**2*g0*g1**2*g5 + b4*b6*g0*g1**2*g5 + b2**2*b3**2*g2*g5 + b1*b3**3*g2*g5 - b2**3*b4*g2*g5 - 4*b1*b2*b3*b4*g2*g5 + 2*b1**2*b4**2*g2*g5 + 3*b1*b2**2*b5*g2*g5 - b1**2*b3*b5*g2*g5 - b1**2*b2*b6*g2*g5 + b3**2*b4*g0*g2*g5 - 3*b2*b4**2*g0*g2*g5 + b2*b3*b5*g0*g2*g5 + 3*b1*b4*b5*g0*g2*g5 + b2**2*b6*g0*g2*g5 - 3*b1*b3*b6*g0*g2*g5 + b5**2*g0**2*g2*g5 - b4*b6*g0**2*g2*g5 - b3*b4**2*g1*g2*g5 + b3**2*b5*g1*g2*g5 + b2*b4*b5*g1*g2*g5 - b1*b5**2*g1*g2*g5 - b2*b3*b6*g1*g2*g5 + b1*b4*b6*g1*g2*g5 + b4**3*g2**2*g5 - 2*b3*b4*b5*g2**2*g5 + b2*b5**2*g2**2*g5 + b3**2*b6*g2**2*g5 - b2*b4*b6*g2**2*g5 - b2*b3**3*g3*g5 + 2*b2**2*b3*b4*g3*g5 + b1*b3**2*b4*g3*g5 - 2*b1*b2*b4**2*g3*g5 - b2**3*b5*g3*g5 + b1**2*b4*b5*g3*g5 + b1*b2**2*b6*g3*g5 - b1**2*b3*b6*g3*g5 + b3*b4**2*g0*g3*g5 - b3**2*b5*g0*g3*g5 - b2*b4*b5*g0*g3*g5 + b1*b5**2*g0*g3*g5 + b2*b3*b6*g0*g3*g5 - b1*b4*b6*g0*g3*g5 - b4**3*g1*g3*g5 + 2*b3*b4*b5*g1*g3*g5 - b2*b5**2*g1*g3*g5 - b3**2*b6*g1*g3*g5 + b2*b4*b6*g1*g3*g5)/den q4 <- (-(b1**4*g0**2) - 3*b1**2*b2*g0**3 - b2**2*g0**4 - 2*b1*b3*g0**4 - b4*g0**5 - b1**5*g1 - 2*b1**3*b2*g0*g1 + 3*b1*b2**2*g0**2*g1 - 3*b1**2*b3*g0**2*g1 + 4*b2*b3*g0**3*g1 - 2*b1*b4*g0**3*g1 + b5*g0**4*g1 + 2*b1**2*b2**2*g1**2 - 3*b1**3*b3*g1**2 - b2**3*g0*g1**2 - 2*b1**2*b4*g0*g1**2 - 3*b3**2*g0**2*g1**2 - b2*b4*g0**2*g1**2 + b1*b5*g0**2*g1**2 - b6*g0**3*g1**2 + b2**2*b3*g1**3 - 2*b1*b3**2*g1**3 + 2*b1*b2*b4*g1**3 - b1**2*b5*g1**3 + 2*b3*b4*g0*g1**3 - 2*b1*b6*g0*g1**3 - b3*b5*g1**4 + b2*b6*g1**4 + b1**4*b2*g2 + b1**2*b2**2*g0*g2 + 3*b1**3*b3*g0*g2 - 2*b2**3*g0**2*g2 + 6*b1*b2*b3*g0**2*g2 + 2*b1**2*b4*g0**2*g2 + 2*b3**2*g0**3*g2 - b2*b4*g0**3*g2 + 3*b1*b5*g0**3*g2 + b6*g0**4*g2 - 2*b1*b2**3*g1*g2 + 3*b1**2*b2*b3*g1*g2 - b1**3*b4*g1*g2 + 2*b2**2*b3*g0*g1*g2 + 2*b1*b3**2*g0*g1*g2 - 8*b1*b2*b4*g0*g1*g2 + 4*b1**2*b5*g0*g1*g2 - b2*b5*g0**2*g1*g2 + b1*b6*g0**2*g1*g2 - b2*b3**2*g1**2*g2 - b2**2*b4*g1**2*g2 + 2*b1*b3*b4*g1**2*g2 + 2*b1*b2*b5*g1**2*g2 - 2*b1**2*b6*g1**2*g2 - 2*b4**2*g0*g1**2*g2 + 3*b3*b5*g0*g1**2*g2 - b2*b6*g0*g1**2*g2 + b4*b5*g1**3*g2 - b3*b6*g1**3*g2 + b1*b2**2*b3*g2**2 - 2*b1**2*b3**2*g2**2 + b1**3*b5*g2**2 - 2*b2*b3**2*g0*g2**2 + 3*b2**2*b4*g0*g2**2 - 2*b1*b3*b4*g0*g2**2 + b1**2*b6*g0*g2**2 + b4**2*g0**2*g2**2 - 3*b3*b5*g0**2*g2**2 + 2*b2*b6*g0**2*g2**2 + 2*b2*b3*b4*g1*g2**2 - 2*b2**2*b5*g1*g2**2 - 2*b1*b3*b5*g1*g2**2 + 2*b1*b2*b6*g1*g2**2 - b5**2*g1**2*g2**2 + b4*b6*g1**2*g2**2 - b2*b4**2*g2**3 + b2*b3*b5*g2**3 + b1*b4*b5*g2**3 - b1*b3*b6*g2**3 + b5**2*g0*g2**3 - b4*b6*g0*g2**3 - b1**3*b2**2*g3 + b1**4*b3*g3 - 3*b1**2*b2*b3*g0*g3 + 3*b1**3*b4*g0*g3 - 3*b1*b3**2*g0**2*g3 + 3*b1**2*b5*g0**2*g3 - 2*b3*b4*g0**3*g3 + b2*b5*g0**3*g3 + b1*b6*g0**3*g3 + 3*b1**2*b3**2*g1*g3 - 3*b1**2*b2*b4*g1*g3 + 6*b1*b3*b4*g0*g1*g3 - 6*b1*b2*b5*g0*g1*g3 + 3*b4**2*g0**2*g1*g3 - 3*b2*b6*g0**2*g1*g3 + b3**3*g1**2*g3 - 2*b2*b3*b4*g1**2*g3 - 2*b1*b4**2*g1**2*g3 + b2**2*b5*g1**2*g3 + 2*b1*b3*b5*g1**2*g3 - 3*b4*b5*g0*g1**2*g3 + 3*b3*b6*g0*g1**2*g3 + b5**2*g1**3*g3 - b4*b6*g1**3*g3 + 2*b3**3*g0*g2*g3 - 4*b2*b3*b4*g0*g2*g3 + 2*b1*b4**2*g0*g2*g3 + 2*b2**2*b5*g0*g2*g3 - 2*b1*b3*b5*g0*g2*g3 - 2*b3**2*b4*g1*g2*g3 + 2*b2*b4**2*g1*g2*g3 + 2*b2*b3*b5*g1*g2*g3 - 2*b1*b4*b5*g1*g2*g3 - 2*b2**2*b6*g1*g2*g3 + 2*b1*b3*b6*g1*g2*g3 + b3*b4**2*g2**2*g3 - b3**2*b5*g2**2*g3 - b2*b4*b5*g2**2*g3 + b1*b5**2*g2**2*g3 + b2*b3*b6*g2**2*g3 - b1*b4*b6*g2**2*g3 - b1*b3**3*g3**2 + 2*b1*b2*b3*b4*g3**2 - b1**2*b4**2*g3**2 - b1*b2**2*b5*g3**2 + b1**2*b3*b5*g3**2 - b3**2*b4*g0*g3**2 + b2*b4**2*g0*g3**2 + b2*b3*b5*g0*g3**2 - b1*b4*b5*g0*g3**2 - b2**2*b6*g0*g3**2 + b1*b3*b6*g0*g3**2 + b3*b4**2*g1*g3**2 - b3**2*b5*g1*g3**2 - b2*b4*b5*g1*g3**2 + b1*b5**2*g1*g3**2 + b2*b3*b6*g1*g3**2 - b1*b4*b6*g1*g3**2 - b4**3*g2*g3**2 + 2*b3*b4*b5*g2*g3**2 - b2*b5**2*g2*g3**2 - b3**2*b6*g2*g3**2 + b2*b4*b6*g2*g3**2 + b1**2*b2**3*g4 - 2*b1**3*b2*b3*g4 + b1**4*b4*g4 - b2**4*g0*g4 + 5*b1*b2**2*b3*g0*g4 - 4*b1**2*b3**2*g0*g4 - 3*b1**2*b2*b4*g0*g4 + 3*b1**3*b5*g0*g4 + 3*b2*b3**2*g0**2*g4 - 2*b2**2*b4*g0**2*g4 - 6*b1*b3*b4*g0**2*g4 + 4*b1*b2*b5*g0**2*g4 + b1**2*b6*g0**2*g4 - 2*b4**2*g0**3*g4 + b3*b5*g0**3*g4 + b2*b6*g0**3*g4 + b2**3*b3*g1*g4 - 4*b1*b2*b3**2*g1*g4 + b1*b2**2*b4*g1*g4 + 4*b1**2*b3*b4*g1*g4 - b1**2*b2*b5*g1*g4 - b1**3*b6*g1*g4 - 4*b3**3*g0*g1*g4 + 4*b2*b3*b4*g0*g1*g4 + 2*b1*b4**2*g0*g1*g4 + 2*b1*b3*b5*g0*g1*g4 - 4*b1*b2*b6*g0*g1*g4 + 3*b4*b5*g0**2*g1*g4 - 3*b3*b6*g0**2*g1*g4 + 2*b3**2*b4*g1**2*g4 - b2*b4**2*g1**2*g4 - 3*b2*b3*b5*g1**2*g4 + b1*b4*b5*g1**2*g4 + 2*b2**2*b6*g1**2*g4 - b1*b3*b6*g1**2*g4 - b5**2*g0*g1**2*g4 + b4*b6*g0*g1**2*g4 - b2**2*b3**2*g2*g4 + 3*b1*b3**3*g2*g4 + b2**3*b4*g2*g4 - 4*b1*b2*b3*b4*g2*g4 + 2*b1**2*b4**2*g2*g4 + b1*b2**2*b5*g2*g4 - 3*b1**2*b3*b5*g2*g4 + b1**2*b2*b6*g2*g4 + 3*b3**2*b4*g0*g2*g4 - b2*b4**2*g0*g2*g4 - 5*b2*b3*b5*g0*g2*g4 + b1*b4*b5*g0*g2*g4 + 3*b2**2*b6*g0*g2*g4 - b1*b3*b6*g0*g2*g4 - b5**2*g0**2*g2*g4 + b4*b6*g0**2*g2*g4 - 3*b3*b4**2*g1*g2*g4 + 3*b3**2*b5*g1*g2*g4 + 3*b2*b4*b5*g1*g2*g4 - 3*b1*b5**2*g1*g2*g4 - 3*b2*b3*b6*g1*g2*g4 + 3*b1*b4*b6*g1*g2*g4 + b4**3*g2**2*g4 - 2*b3*b4*b5*g2**2*g4 + b2*b5**2*g2**2*g4 + b3**2*b6*g2**2*g4 - b2*b4*b6*g2**2*g4 + b2*b3**3*g3*g4 - 2*b2**2*b3*b4*g3*g4 - b1*b3**2*b4*g3*g4 + 2*b1*b2*b4**2*g3*g4 + b2**3*b5*g3*g4 - b1**2*b4*b5*g3*g4 - b1*b2**2*b6*g3*g4 + b1**2*b3*b6*g3*g4 - b3*b4**2*g0*g3*g4 + b3**2*b5*g0*g3*g4 + b2*b4*b5*g0*g3*g4 - b1*b5**2*g0*g3*g4 - b2*b3*b6*g0*g3*g4 + b1*b4*b6*g0*g3*g4 + b4**3*g1*g3*g4 - 2*b3*b4*b5*g1*g3*g4 + b2*b5**2*g1*g3*g4 + b3**2*b6*g1*g3*g4 - b2*b4*b6*g1*g3*g4 - b3**4*g4**2 + 3*b2*b3**2*b4*g4**2 - b2**2*b4**2*g4**2 - 2*b1*b3*b4**2*g4**2 - 2*b2**2*b3*b5*g4**2 + 2*b1*b3**2*b5*g4**2 + 2*b1*b2*b4*b5*g4**2 - b1**2*b5**2*g4**2 + b2**3*b6*g4**2 - 2*b1*b2*b3*b6*g4**2 + b1**2*b4*b6*g4**2 - b4**3*g0*g4**2 + 2*b3*b4*b5*g0*g4**2 - b2*b5**2*g0*g4**2 - b3**2*b6*g0*g4**2 + b2*b4*b6*g0*g4**2 - b1*b2**4*g5 + 3*b1**2*b2**2*b3*g5 - b1**3*b3**2*g5 - 2*b1**3*b2*b4*g5 + b1**4*b5*g5 - b2**3*b3*g0*g5 + 4*b1*b2*b3**2*g0*g5 - b1*b2**2*b4*g0*g5 - 4*b1**2*b3*b4*g0*g5 + b1**2*b2*b5*g0*g5 + b1**3*b6*g0*g5 + b3**3*g0**2*g5 - 2*b1*b4**2*g0**2*g5 - b2**2*b5*g0**2*g5 + 2*b1*b2*b6*g0**2*g5 - b4*b5*g0**3*g5 + b3*b6*g0**3*g5 + b2**2*b3**2*g1*g5 - 2*b1*b3**3*g1*g5 - b2**3*b4*g1*g5 + 2*b1*b2*b3*b4*g1*g5 - b1**2*b4**2*g1*g5 + 2*b1**2*b3*b5*g1*g5 - b1**2*b2*b6*g1*g5 - 2*b3**2*b4*g0*g1*g5 + 4*b2*b3*b5*g0*g1*g5 - 2*b2**2*b6*g0*g1*g5 + b5**2*g0**2*g1*g5 - b4*b6*g0**2*g1*g5 + b3*b4**2*g1**2*g5 - b3**2*b5*g1**2*g5 - b2*b4*b5*g1**2*g5 + b1*b5**2*g1**2*g5 + b2*b3*b6*g1**2*g5 - b1*b4*b6*g1**2*g5 - b2*b3**3*g2*g5 + 2*b2**2*b3*b4*g2*g5 + b1*b3**2*b4*g2*g5 - 2*b1*b2*b4**2*g2*g5 - b2**3*b5*g2*g5 + b1**2*b4*b5*g2*g5 + b1*b2**2*b6*g2*g5 - b1**2*b3*b6*g2*g5 + b3*b4**2*g0*g2*g5 - b3**2*b5*g0*g2*g5 - b2*b4*b5*g0*g2*g5 + b1*b5**2*g0*g2*g5 + b2*b3*b6*g0*g2*g5 - b1*b4*b6*g0*g2*g5 - b4**3*g1*g2*g5 + 2*b3*b4*b5*g1*g2*g5 - b2*b5**2*g1*g2*g5 - b3**2*b6*g1*g2*g5 + b2*b4*b6*g1*g2*g5 + b3**4*g3*g5 - 3*b2*b3**2*b4*g3*g5 + b2**2*b4**2*g3*g5 + 2*b1*b3*b4**2*g3*g5 + 2*b2**2*b3*b5*g3*g5 - 2*b1*b3**2*b5*g3*g5 - 2*b1*b2*b4*b5*g3*g5 + b1**2*b5**2*g3*g5 - b2**3*b6*g3*g5 + 2*b1*b2*b3*b6*g3*g5 - b1**2*b4*b6*g3*g5 + b4**3*g0*g3*g5 - 2*b3*b4*b5*g0*g3*g5 + b2*b5**2*g0*g3*g5 + b3**2*b6*g0*g3*g5 - b2*b4*b6*g0*g3*g5)/den q5 <- (-(b1**5*g0) - 4*b1**3*b2*g0**2 - 3*b1*b2**2*g0**3 - 3*b1**2*b3*g0**3 - 2*b2*b3*g0**4 - 2*b1*b4*g0**4 - b5*g0**5 + b1**4*b2*g1 + 6*b1**2*b2**2*g0*g1 - 2*b1**3*b3*g0*g1 + 3*b2**3*g0**2*g1 + 6*b1*b2*b3*g0**2*g1 - 3*b1**2*b4*g0**2*g1 + 2*b3**2*g0**3*g1 + 4*b2*b4*g0**3*g1 - 2*b1*b5*g0**3*g1 + b6*g0**4*g1 - 2*b1*b2**3*g1**2 + b1**2*b2*b3*g1**2 + b1**3*b4*g1**2 - 5*b2**2*b3*g0*g1**2 - 2*b1*b3**2*g0*g1**2 + 6*b1*b2*b4*g0*g1**2 + b1**2*b5*g0*g1**2 - 4*b3*b4*g0**2*g1**2 + b2*b5*g0**2*g1**2 + 3*b1*b6*g0**2*g1**2 + 2*b2*b3**2*g1**3 - b2**2*b4*g1**3 - 2*b1*b2*b5*g1**3 + b1**2*b6*g1**3 + 2*b4**2*g0*g1**3 - 2*b2*b6*g0*g1**3 - b4*b5*g1**4 + b3*b6*g1**4 - b1**3*b2**2*g2 + b1**4*b3*g2 - 4*b1*b2**3*g0*g2 + 5*b1**2*b2*b3*g0*g2 - b1**3*b4*g0*g2 - 4*b2**2*b3*g0**2*g2 + 5*b1*b3**2*g0**2*g2 - b1**2*b5*g0**2*g2 + 2*b3*b4*g0**3*g2 - 3*b2*b5*g0**3*g2 + b1*b6*g0**3*g2 + 2*b2**4*g1*g2 - 2*b1*b2**2*b3*g1*g2 + b1**2*b3**2*g1*g2 - 3*b1**2*b2*b4*g1*g2 + 2*b1**3*b5*g1*g2 + 4*b2*b3**2*g0*g1*g2 - 6*b1*b3*b4*g0*g1*g2 - 2*b1*b2*b5*g0*g1*g2 + 4*b1**2*b6*g0*g1*g2 - 3*b4**2*g0**2*g1*g2 + 2*b3*b5*g0**2*g1*g2 + b2*b6*g0**2*g1*g2 - b3**3*g1**2*g2 - 2*b2*b3*b4*g1**2*g2 + 3*b2**2*b5*g1**2*g2 + 4*b1*b3*b5*g1**2*g2 - 4*b1*b2*b6*g1**2*g2 + b4*b5*g0*g1**2*g2 - b3*b6*g0*g1**2*g2 + b5**2*g1**3*g2 - b4*b6*g1**3*g2 - b2**3*b3*g2**2 + 3*b1*b2**2*b4*g2**2 - 3*b1**2*b2*b5*g2**2 + b1**3*b6*g2**2 - 2*b3**3*g0*g2**2 + 4*b2*b3*b4*g0*g2**2 - 2*b2**2*b5*g0*g2**2 + b4*b5*g0**2*g2**2 - b3*b6*g0**2*g2**2 + 2*b3**2*b4*g1*g2**2 - 4*b2*b3*b5*g1*g2**2 + 2*b2**2*b6*g1*g2**2 - 2*b5**2*g0*g1*g2**2 + 2*b4*b6*g0*g1*g2**2 - b3*b4**2*g2**3 + b3**2*b5*g2**3 + b2*b4*b5*g2**3 - b1*b5**2*g2**3 - b2*b3*b6*g2**3 + b1*b4*b6*g2**3 + b1**2*b2**3*g3 - 2*b1**3*b2*b3*g3 + b1**4*b4*g3 + 2*b2**4*g0*g3 - 4*b1*b2**2*b3*g0*g3 - b1**2*b3**2*g0*g3 + 3*b1**2*b2*b4*g0*g3 - 3*b2*b3**2*g0**2*g3 + 4*b2**2*b4*g0**2*g3 - 2*b1*b2*b5*g0**2*g3 + b1**2*b6*g0**2*g3 + b4**2*g0**3*g3 - 2*b3*b5*g0**3*g3 + b2*b6*g0**3*g3 - 2*b2**3*b3*g1*g3 + 2*b1*b2*b3**2*g1*g3 + 4*b1*b2**2*b4*g1*g3 - 2*b1**2*b3*b4*g1*g3 - 4*b1**2*b2*b5*g1*g3 + 2*b1**3*b6*g1*g3 + 2*b3**3*g0*g1*g3 - 2*b2*b3*b4*g0*g1*g3 + 2*b1*b4**2*g0*g1*g3 - 4*b1*b3*b5*g0*g1*g3 + 2*b1*b2*b6*g0*g1*g3 - b3**2*b4*g1**2*g3 + 2*b2*b4**2*g1**2*g3 - 2*b1*b4*b5*g1**2*g3 - b2**2*b6*g1**2*g3 + 2*b1*b3*b6*g1**2*g3 - b5**2*g0*g1**2*g3 + b4*b6*g0*g1**2*g3 + 2*b2**2*b3**2*g2*g3 - 2*b2**3*b4*g2*g3 - 4*b1*b2*b3*b4*g2*g3 + 2*b1**2*b4**2*g2*g3 + 4*b1*b2**2*b5*g2*g3 - 2*b1**2*b2*b6*g2*g3 - 4*b2*b4**2*g0*g2*g3 + 4*b2*b3*b5*g0*g2*g3 + 4*b1*b4*b5*g0*g2*g3 - 4*b1*b3*b6*g0*g2*g3 + 2*b5**2*g0**2*g2*g3 - 2*b4*b6*g0**2*g2*g3 + b4**3*g2**2*g3 - 2*b3*b4*b5*g2**2*g3 + b2*b5**2*g2**2*g3 + b3**2*b6*g2**2*g3 - b2*b4*b6*g2**2*g3 - b2*b3**3*g3**2 + 2*b2**2*b3*b4*g3**2 + b1*b3**2*b4*g3**2 - 2*b1*b2*b4**2*g3**2 - b2**3*b5*g3**2 + b1**2*b4*b5*g3**2 + b1*b2**2*b6*g3**2 - b1**2*b3*b6*g3**2 + b3*b4**2*g0*g3**2 - b3**2*b5*g0*g3**2 - b2*b4*b5*g0*g3**2 + b1*b5**2*g0*g3**2 + b2*b3*b6*g0*g3**2 - b1*b4*b6*g0*g3**2 - b4**3*g1*g3**2 + 2*b3*b4*b5*g1*g3**2 - b2*b5**2*g1*g3**2 - b3**2*b6*g1*g3**2 + b2*b4*b6*g1*g3**2 - b1*b2**4*g4 + 3*b1**2*b2**2*b3*g4 - b1**3*b3**2*g4 - 2*b1**3*b2*b4*g4 + b1**4*b5*g4 - b2**3*b3*g0*g4 + 4*b1*b2*b3**2*g0*g4 - b1*b2**2*b4*g0*g4 - 4*b1**2*b3*b4*g0*g4 + b1**2*b2*b5*g0*g4 + b1**3*b6*g0*g4 + b3**3*g0**2*g4 - 2*b1*b4**2*g0**2*g4 - b2**2*b5*g0**2*g4 + 2*b1*b2*b6*g0**2*g4 - b4*b5*g0**3*g4 + b3*b6*g0**3*g4 + b2**2*b3**2*g1*g4 - 2*b1*b3**3*g1*g4 - b2**3*b4*g1*g4 + 2*b1*b2*b3*b4*g1*g4 - b1**2*b4**2*g1*g4 + 2*b1**2*b3*b5*g1*g4 - b1**2*b2*b6*g1*g4 - 2*b3**2*b4*g0*g1*g4 + 4*b2*b3*b5*g0*g1*g4 - 2*b2**2*b6*g0*g1*g4 + b5**2*g0**2*g1*g4 - b4*b6*g0**2*g1*g4 + b3*b4**2*g1**2*g4 - b3**2*b5*g1**2*g4 - b2*b4*b5*g1**2*g4 + b1*b5**2*g1**2*g4 + b2*b3*b6*g1**2*g4 - b1*b4*b6*g1**2*g4 - b2*b3**3*g2*g4 + 2*b2**2*b3*b4*g2*g4 + b1*b3**2*b4*g2*g4 - 2*b1*b2*b4**2*g2*g4 - b2**3*b5*g2*g4 + b1**2*b4*b5*g2*g4 + b1*b2**2*b6*g2*g4 - b1**2*b3*b6*g2*g4 + b3*b4**2*g0*g2*g4 - b3**2*b5*g0*g2*g4 - b2*b4*b5*g0*g2*g4 + b1*b5**2*g0*g2*g4 + b2*b3*b6*g0*g2*g4 - b1*b4*b6*g0*g2*g4 - b4**3*g1*g2*g4 + 2*b3*b4*b5*g1*g2*g4 - b2*b5**2*g1*g2*g4 - b3**2*b6*g1*g2*g4 + b2*b4*b6*g1*g2*g4 + b3**4*g3*g4 - 3*b2*b3**2*b4*g3*g4 + b2**2*b4**2*g3*g4 + 2*b1*b3*b4**2*g3*g4 + 2*b2**2*b3*b5*g3*g4 - 2*b1*b3**2*b5*g3*g4 - 2*b1*b2*b4*b5*g3*g4 + b1**2*b5**2*g3*g4 - b2**3*b6*g3*g4 + 2*b1*b2*b3*b6*g3*g4 - b1**2*b4*b6*g3*g4 + b4**3*g0*g3*g4 - 2*b3*b4*b5*g0*g3*g4 + b2*b5**2*g0*g3*g4 + b3**2*b6*g0*g3*g4 - b2*b4*b6*g0*g3*g4 + b2**5*g5 - 4*b1*b2**3*b3*g5 + 3*b1**2*b2*b3**2*g5 + 3*b1**2*b2**2*b4*g5 - 2*b1**3*b3*b4*g5 - 2*b1**3*b2*b5*g5 + b1**4*b6*g5 - 3*b2**2*b3**2*g0*g5 + 2*b1*b3**3*g0*g5 + 3*b2**3*b4*g0*g5 + 2*b1*b2*b3*b4*g0*g5 - b1**2*b4**2*g0*g5 - 4*b1*b2**2*b5*g0*g5 - 2*b1**2*b3*b5*g0*g5 + 3*b1**2*b2*b6*g0*g5 + b3**2*b4*g0**2*g5 + 2*b2*b4**2*g0**2*g5 - 4*b2*b3*b5*g0**2*g5 - 2*b1*b4*b5*g0**2*g5 + b2**2*b6*g0**2*g5 + 2*b1*b3*b6*g0**2*g5 - b5**2*g0**3*g5 + b4*b6*g0**3*g5 + 2*b2*b3**3*g1*g5 - 4*b2**2*b3*b4*g1*g5 - 2*b1*b3**2*b4*g1*g5 + 4*b1*b2*b4**2*g1*g5 + 2*b2**3*b5*g1*g5 - 2*b1**2*b4*b5*g1*g5 - 2*b1*b2**2*b6*g1*g5 + 2*b1**2*b3*b6*g1*g5 - 2*b3*b4**2*g0*g1*g5 + 2*b3**2*b5*g0*g1*g5 + 2*b2*b4*b5*g0*g1*g5 - 2*b1*b5**2*g0*g1*g5 - 2*b2*b3*b6*g0*g1*g5 + 2*b1*b4*b6*g0*g1*g5 + b4**3*g1**2*g5 - 2*b3*b4*b5*g1**2*g5 + b2*b5**2*g1**2*g5 + b3**2*b6*g1**2*g5 - b2*b4*b6*g1**2*g5 - b3**4*g2*g5 + 3*b2*b3**2*b4*g2*g5 - b2**2*b4**2*g2*g5 - 2*b1*b3*b4**2*g2*g5 - 2*b2**2*b3*b5*g2*g5 + 2*b1*b3**2*b5*g2*g5 + 2*b1*b2*b4*b5*g2*g5 - b1**2*b5**2*g2*g5 + b2**3*b6*g2*g5 - 2*b1*b2*b3*b6*g2*g5 + b1**2*b4*b6*g2*g5 - b4**3*g0*g2*g5 + 2*b3*b4*b5*g0*g2*g5 - b2*b5**2*g0*g2*g5 - b3**2*b6*g0*g2*g5 + b2*b4*b6*g0*g2*g5)/den q6 <- (-b1**6 - 5*b1**4*b2*g0 - 6*b1**2*b2**2*g0**2 - 4*b1**3*b3*g0**2 - b2**3*g0**3 - 6*b1*b2*b3*g0**3 - 3*b1**2*b4*g0**3 - b3**2*g0**4 - 2*b2*b4*g0**4 - 2*b1*b5*g0**4 - b6*g0**5 + 4*b1**3*b2**2*g1 - 4*b1**4*b3*g1 + 6*b1*b2**3*g0*g1 - 6*b1**3*b4*g0*g1 + 6*b2**2*b3*g0**2*g1 - 6*b1**2*b5*g0**2*g1 + 2*b3*b4*g0**3*g1 + 2*b2*b5*g0**3*g1 - 4*b1*b6*g0**3*g1 - b2**4*g1**2 - 4*b1**2*b3**2*g1**2 + 7*b1**2*b2*b4*g1**2 - 2*b1**3*b5*g1**2 - 4*b2*b3**2*g0*g1**2 + b2**2*b4*g0*g1**2 - 2*b1*b3*b4*g0*g1**2 + 8*b1*b2*b5*g0*g1**2 - 3*b1**2*b6*g0*g1**2 - b4**2*g0**2*g1**2 - 2*b3*b5*g0**2*g1**2 + 3*b2*b6*g0**2*g1**2 + 2*b2*b3*b4*g1**3 + 2*b1*b4**2*g1**3 - 2*b2**2*b5*g1**3 - 4*b1*b3*b5*g1**3 + 2*b1*b2*b6*g1**3 + 2*b4*b5*g0*g1**3 - 2*b3*b6*g0*g1**3 - b5**2*g1**4 + b4*b6*g1**4 - 3*b1**2*b2**3*g2 + 6*b1**3*b2*b3*g2 - 3*b1**4*b4*g2 - 2*b2**4*g0*g2 + 7*b1**2*b3**2*g0*g2 - b1**2*b2*b4*g0*g2 - 4*b1**3*b5*g0*g2 + b2*b3**2*g0**2*g2 - 4*b2**2*b4*g0**2*g2 + 8*b1*b3*b4*g0**2*g2 - 2*b1*b2*b5*g0**2*g2 - 3*b1**2*b6*g0**2*g2 + b4**2*g0**3*g2 + 2*b3*b5*g0**3*g2 - 3*b2*b6*g0**3*g2 + 2*b2**3*b3*g1*g2 + 2*b1*b2*b3**2*g1*g2 - 8*b1*b2**2*b4*g1*g2 - 2*b1**2*b3*b4*g1*g2 + 8*b1**2*b2*b5*g1*g2 - 2*b1**3*b6*g1*g2 + 2*b3**3*g0*g1*g2 - 2*b2*b3*b4*g0*g1*g2 - 6*b1*b4**2*g0*g1*g2 + 4*b1*b3*b5*g0*g1*g2 + 2*b1*b2*b6*g0*g1*g2 - 4*b4*b5*g0**2*g1*g2 + 4*b3*b6*g0**2*g1*g2 - b3**2*b4*g1**2*g2 - 2*b2*b4**2*g1**2*g2 + 4*b2*b3*b5*g1**2*g2 + 2*b1*b4*b5*g1**2*g2 - b2**2*b6*g1**2*g2 - 2*b1*b3*b6*g1**2*g2 + 3*b5**2*g0*g1**2*g2 - 3*b4*b6*g0*g1**2*g2 - b2**2*b3**2*g2**2 - 2*b1*b3**3*g2**2 + b2**3*b4*g2**2 + 6*b1*b2*b3*b4*g2**2 - 3*b1**2*b4**2*g2**2 - 4*b1*b2**2*b5*g2**2 + 2*b1**2*b3*b5*g2**2 + b1**2*b2*b6*g2**2 - 2*b3**2*b4*g0*g2**2 + 4*b2*b4**2*g0*g2**2 - 4*b1*b4*b5*g0*g2**2 - 2*b2**2*b6*g0*g2**2 + 4*b1*b3*b6*g0*g2**2 - b5**2*g0**2*g2**2 + b4*b6*g0**2*g2**2 + 2*b3*b4**2*g1*g2**2 - 2*b3**2*b5*g1*g2**2 - 2*b2*b4*b5*g1*g2**2 + 2*b1*b5**2*g1*g2**2 + 2*b2*b3*b6*g1*g2**2 - 2*b1*b4*b6*g1*g2**2 - b4**3*g2**3 + 2*b3*b4*b5*g2**3 - b2*b5**2*g2**3 - b3**2*b6*g2**3 + b2*b4*b6*g2**3 + 2*b1*b2**4*g3 - 6*b1**2*b2**2*b3*g3 + 2*b1**3*b3**2*g3 + 4*b1**3*b2*b4*g3 - 2*b1**4*b5*g3 + 2*b2**3*b3*g0*g3 - 8*b1*b2*b3**2*g0*g3 + 2*b1*b2**2*b4*g0*g3 + 8*b1**2*b3*b4*g0*g3 - 2*b1**2*b2*b5*g0*g3 - 2*b1**3*b6*g0*g3 - 2*b3**3*g0**2*g3 + 4*b1*b4**2*g0**2*g3 + 2*b2**2*b5*g0**2*g3 - 4*b1*b2*b6*g0**2*g3 + 2*b4*b5*g0**3*g3 - 2*b3*b6*g0**3*g3 - 2*b2**2*b3**2*g1*g3 + 4*b1*b3**3*g1*g3 + 2*b2**3*b4*g1*g3 - 4*b1*b2*b3*b4*g1*g3 + 2*b1**2*b4**2*g1*g3 - 4*b1**2*b3*b5*g1*g3 + 2*b1**2*b2*b6*g1*g3 + 4*b3**2*b4*g0*g1*g3 - 8*b2*b3*b5*g0*g1*g3 + 4*b2**2*b6*g0*g1*g3 - 2*b5**2*g0**2*g1*g3 + 2*b4*b6*g0**2*g1*g3 - 2*b3*b4**2*g1**2*g3 + 2*b3**2*b5*g1**2*g3 + 2*b2*b4*b5*g1**2*g3 - 2*b1*b5**2*g1**2*g3 - 2*b2*b3*b6*g1**2*g3 + 2*b1*b4*b6*g1**2*g3 + 2*b2*b3**3*g2*g3 - 4*b2**2*b3*b4*g2*g3 - 2*b1*b3**2*b4*g2*g3 + 4*b1*b2*b4**2*g2*g3 + 2*b2**3*b5*g2*g3 - 2*b1**2*b4*b5*g2*g3 - 2*b1*b2**2*b6*g2*g3 + 2*b1**2*b3*b6*g2*g3 - 2*b3*b4**2*g0*g2*g3 + 2*b3**2*b5*g0*g2*g3 + 2*b2*b4*b5*g0*g2*g3 - 2*b1*b5**2*g0*g2*g3 - 2*b2*b3*b6*g0*g2*g3 + 2*b1*b4*b6*g0*g2*g3 + 2*b4**3*g1*g2*g3 - 4*b3*b4*b5*g1*g2*g3 + 2*b2*b5**2*g1*g2*g3 + 2*b3**2*b6*g1*g2*g3 - 2*b2*b4*b6*g1*g2*g3 - b3**4*g3**2 + 3*b2*b3**2*b4*g3**2 - b2**2*b4**2*g3**2 - 2*b1*b3*b4**2*g3**2 - 2*b2**2*b3*b5*g3**2 + 2*b1*b3**2*b5*g3**2 + 2*b1*b2*b4*b5*g3**2 - b1**2*b5**2*g3**2 + b2**3*b6*g3**2 - 2*b1*b2*b3*b6*g3**2 + b1**2*b4*b6*g3**2 - b4**3*g0*g3**2 + 2*b3*b4*b5*g0*g3**2 - b2*b5**2*g0*g3**2 - b3**2*b6*g0*g3**2 + b2*b4*b6*g0*g3**2 - b2**5*g4 + 4*b1*b2**3*b3*g4 - 3*b1**2*b2*b3**2*g4 - 3*b1**2*b2**2*b4*g4 + 2*b1**3*b3*b4*g4 + 2*b1**3*b2*b5*g4 - b1**4*b6*g4 + 3*b2**2*b3**2*g0*g4 - 2*b1*b3**3*g0*g4 - 3*b2**3*b4*g0*g4 - 2*b1*b2*b3*b4*g0*g4 + b1**2*b4**2*g0*g4 + 4*b1*b2**2*b5*g0*g4 + 2*b1**2*b3*b5*g0*g4 - 3*b1**2*b2*b6*g0*g4 - b3**2*b4*g0**2*g4 - 2*b2*b4**2*g0**2*g4 + 4*b2*b3*b5*g0**2*g4 + 2*b1*b4*b5*g0**2*g4 - b2**2*b6*g0**2*g4 - 2*b1*b3*b6*g0**2*g4 + b5**2*g0**3*g4 - b4*b6*g0**3*g4 - 2*b2*b3**3*g1*g4 + 4*b2**2*b3*b4*g1*g4 + 2*b1*b3**2*b4*g1*g4 - 4*b1*b2*b4**2*g1*g4 - 2*b2**3*b5*g1*g4 + 2*b1**2*b4*b5*g1*g4 + 2*b1*b2**2*b6*g1*g4 - 2*b1**2*b3*b6*g1*g4 + 2*b3*b4**2*g0*g1*g4 - 2*b3**2*b5*g0*g1*g4 - 2*b2*b4*b5*g0*g1*g4 + 2*b1*b5**2*g0*g1*g4 + 2*b2*b3*b6*g0*g1*g4 - 2*b1*b4*b6*g0*g1*g4 - b4**3*g1**2*g4 + 2*b3*b4*b5*g1**2*g4 - b2*b5**2*g1**2*g4 - b3**2*b6*g1**2*g4 + b2*b4*b6*g1**2*g4 + b3**4*g2*g4 - 3*b2*b3**2*b4*g2*g4 + b2**2*b4**2*g2*g4 + 2*b1*b3*b4**2*g2*g4 + 2*b2**2*b3*b5*g2*g4 - 2*b1*b3**2*b5*g2*g4 - 2*b1*b2*b4*b5*g2*g4 + b1**2*b5**2*g2*g4 - b2**3*b6*g2*g4 + 2*b1*b2*b3*b6*g2*g4 - b1**2*b4*b6*g2*g4 + b4**3*g0*g2*g4 - 2*b3*b4*b5*g0*g2*g4 + b2*b5**2*g0*g2*g4 + b3**2*b6*g0*g2*g4 - b2*b4*b6*g0*g2*g4)/den p1 <- b1 p2 <- b2 + b1*q1 p3 <- b3 + b1*q2 + b2*q1 p4 <- b4 + b3*q1 + b2*q2 + b1*q3 p5 <- b5 + b4*q1 + b3*q2 + b2*q3 + b1*q4 p6 <- g0*q6 y <- x^al h.pade <- (p1 * y + p2 * y^2 + p3 * y^3 + p4 * y^4 + p5 * y^5 + p6 * y^6)/ (1 + q1 * y + q2 * y^2 + q3 * y^3 + q4 * y^4 + q5 * y^5 + q6 * y^6) return(h.pade) } ######################################################################## # Padé approximations to D^\alpha h(a,t) ######################################################################## d.h.Pade22 <- function(params)function(a,tau){ H <- params$H rho <- params$rho nu <- params$nu al <- H+1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp -(0+1i) * rho * a aa <- sqrt(a * (a + (0+1i)) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa h.pade <- h.Pade22(params)(a,tau) res <- 1/2*(nu*h.pade-rm)*(nu*h.pade-rp) return(res) } ######################################################################## d.h.Pade33 <- function(params)function(a,tau){ rho <- params$rho nu <- params$nu lam <- params$lam lamp <- lam/nu lamTilde <- lamp -(0+1i) * rho * a aa <- sqrt(a * (a + (0+1i)) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa h.pade <- h.Pade33(params)(a,tau) res <- 1/2*(nu*h.pade-rm)*(nu*h.pade-rp) return(res) } ######################################################################## d.h.Pade44 <- function(params)function(a,tau){ H <- params$H rho <- params$rho nu <- params$nu al <- H+1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp -(0+1i) * rho * a aa <- sqrt(a * (a + (0+1i)) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa h.pade <- h.Pade44(params)(a,tau) res <- 1/2*(nu*h.pade-rm)*(nu*h.pade-rp) return(res) } ######################################################################## d.h.Pade55 <- function(params)function(a,tau){ H <- params$H rho <- params$rho nu <- params$nu al <- H+1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp -(0+1i) * rho * a aa <- sqrt(a * (a + (0+1i)) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa h.pade <- h.Pade55(params)(a,tau) res <- 1/2*(nu*h.pade-rm)*(nu*h.pade-rp) return(res) } ######################################################################## d.h.Pade66 <- function (params) function(a, tau) { H <- params$H rho <- params$rho nu <- params$nu al <- H + 1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp - (0 + (0+1i)) * rho * a aa <- sqrt(a * (a + (0 + (0+1i))) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa h.pade <- h.Pade66(params)(a, tau) res <- 1/2 * (nu * h.pade - rm) * (nu * h.pade - rp) return(res) } ######################################################################## # # Characteristic function using Padé approximation # ######################################################################## phiRoughHestonRational.raw <- function(params, xiCurve, h.approx, n=100) function(a, tau) { H <- params$H rho <- params$rho nu <- params$nu al <- H+1/2 lam <- params$lam lamp <- lam/nu lamTilde <- lamp -(0+1i) * rho * a aa <- sqrt(a * (a + (0+1i)) + lamTilde^2) rm <- lamTilde - aa rp <- lamTilde + aa ti <- (0:n)/n * tau h <- h.approx(params)(a,ti) dah <- 1/2*(nu*h-rm)*(nu*h-rp) g <- dah + lam * h xi <- xiCurve(ti) conv <- t(g) %*% rev(xi) * tau cor <- (g[1]*xi[n+1]+g[n+1]*xi[1])/2*tau psi <- (conv-cor)/n return(exp(psi)) } phiRoughHestonRational <- function(params, xiCurve, h.approx, n=100) function(a, tau){ phi1 <- function(u){ifelse(u==0,1,phiRoughHestonRational.raw(params, xiCurve, h.approx, n)(u,tau))} return(sapply(a,phi1)) }