Help with implementing asymptotic expansion for fractional-power polynomials in R

I've been developing an R script to compute asymptotic expansions of a
polynomial raised to a fractional power in order to approximate that
integral.

he implementation is my own, but I'm encountering issues with higher-order
terms and fractional exponents. I'd like guidance on making the code:

Robust for any polynomial degree
Accurate with fractional powers
Efficient for high-degree polynomials

library(partitions)

asymptotic_expansion <- function(coeffs, alpha, max_order = 4) {
  degree <- length(coeffs) - 1
  r_lead <- coeffs[1]
  lower_coeffs <- coeffs[-1] / r_lead
  n <- length(lower_coeffs)
  
  result <- numeric(max_order + 1)
  result[1] <- 1
  names(result) <- paste0("x^", degree*alpha - 0:max_order)
  
  for (k in 1:max_order) {
    combos <- compositions(k, n)
    for (col in 1:ncol(combos)) {
      j <- combos[, col]
      power_y <- sum(j * (1:n))
      if (power_y <= max_order) {
        multinom_coef <- factorial(k) / prod(factorial(j))
        term <- multinom_coef * prod(lower_coeffs^j)
        result[power_y + 1] <- result[power_y + 1] + choose(alpha, k) * term
      }
    }
  }
  
  result <- result * r_lead^alpha
  return(result)
}

# Example polynomial and parameters
coeffs <- c(1, -17, 80, -100)
alpha <- 0.01
max_order <- 4

# Compute expansion
coefs <- asymptotic_expansion(coeffs, alpha, max_order)
print(coefs)

x^0.03   x^-0.97   x^-1.97   x^-2.97   x^-3.97 
1.0         -0.17         0.798    -3.667    ...

The coefficient for the fourth term (-3.667) seems inconsistent with
expected results from higher-precision calculations or reference tools
(should be roughly -1.656). I suspect this may be related to fractional
exponents or truncation of higher-order terms.

I'd appreciate guidance on:

Correctly handling fractional exponents in the expansion
Improving accuracy for higher-order terms
Making the code general for any polynomial degree

Thank you for your help!

Thomas Subia