Variable order in large-order expansions (at least in `besselk`)

I finally got around to cleaning up the uniform asymptotic expansion code in `BesselK.jl` and I ended up with a solution that is within a few `ns` of `Bessels.besselk_large_orders` while still allowing a user to pick the order of the expansion. I'm wondering if any part of this new work would be interesting as a PR here. There are a few parts to it:

1. Auto-generating the U_k polynomials, which I do with [this file](https://github.com/cgeoga/BesselK.jl/blob/main/src/genUk/generate_uk_polys_hardcoded.jl). The result is that you automatically generate a file that looks like [this](https://github.com/cgeoga/BesselK.jl/blob/main/src/uk_polys.jl), and you can see that there's a struct for each polynomial that trims out zeros and has a slightly specialized eval method. The generation uses `Polynomials.jl`, but once the file is generated you don't need it, and so it would not be a new dependency.

2. The actual equivalent to `Bessels.besselk_large_orders` now looks like this ([file ref](https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_asv.jl)):
```julia
@generated function _besselk_vz_asv(v, z, ::Val{N}) where{N}
quote
ez = sqrt(1+z*z)+log(z/(1+sqrt(1+z*z)))
pz = inv(sqrt(1+z*z))
out = sqrt(pi/(2*v))/sqrt(sqrt(1+z*z))
(ser, sgn, _v) = (zero(z), one(z), one(v))
evaled_polys = tuple($([:($(Symbol(:uk_, j, :_poly))(pz)) for j in 0:(N-1)]...))
Base.Cartesian.@nexprs $N j -> begin
ser += sgn*evaled_polys[j]/_v
sgn *= -one(z)
_v *= v
end
exp(-v*ez)*out*ser
end
end

_besselk_asv(v, z, maxorder) = _besselk_vz_asv(v, z/v, maxorder)
```
Which maybe is objectionable because it is generated function that automatically generates and uses function names. But it does work, and the `@nexprs` is a way to be completely sure the loop is unrolling and stuff. I have a feeling that with some smarter inner-loop tweaks this could at least match the current performance.

There are a few potential advantages to this:
1. Maybe there is some performance on the table where one could more gradually lower the order of the approximation or something.
2. The function that generates the `U_k` file uses `Polynomials.jl`, and so one could actually generate exact rational representations for all of the `U_k` polynomial coefficients and potentially use something like that to get more accurate coefficients for big floats or whtatever.

I haven't really explored either of those options. But if they're potentially interesting, I'm happy to prepare a PR that puts in the generated U_k polynomials and stuff. I'm sure other people here will have thoughts about the custom struct and methods and stuff, so I'm happy to discuss all aspects of it if there is any interest in the first place.