Re: How to compute functions defined by recurrences?

And indeed

-5/36 * prodnumrat(1-5/24/k/(3*(k-1)/2),2)

is about -0.1205, so my empirical k was (approximately) correct.

I recognize that this doesn't address the part of your question about computing the whole polynomial for large n. My code generates D_250(x) in about 8 seconds, which isn't bad considering the output is around 1 MB (a bit more printed, a bit less internal). I'm sure there are clever ways to do it better but if you need a bit output it's never going to be super fast.

On Mon, Sep 23, 2024 at 12:36 PM Charles Greathouse <crgreathouse@gmail.com> wrote:

I'm not afraid of inflicting bad code (or embarrassing myself) on the list, so I'll give

D(n)=if(n==0, return(1)); my(d=D(n-1)); x^2*(1-x)^2*(d')/2 + intformal((1-5*x^2)*d)/8

which I think is the naive transcription of the definition. Let's use f(n) as the leading coefficient of D_n(x). Then, speaking purely numerically, it looks like
f(n)/f(n-1) = 1.5n + C
where the constant looks suspiciously like -1.5. If that's right, you should have a factorial-like function as the leading coefficient, something like
f(n) = -k * (n-1)! * (3/2)^n
Numerically it seems to work pretty well with a constant k near 0.120.

OK, so how did we do? From the defining formula, it looks like

f(n) = prod(k=1,n, (3*k-3)/2-5/24/k)

which indeed has the (n-1)! * (3/2)^n but I missed the correction term.

On Mon, Sep 23, 2024 at 11:06 AM kevin lucas <lucaskevin296@gmail.com> wrote:
In some recent work I found myself having to investigate the asymptotics of the Debye polynomials, which I will here take as *defined* by the following integrodifferential recurrence:

$D_{n+1}(x)= \frac{1}{2}x^2(1-x)^2D'_n(x) + \frac{1}{8}\int_0^x(1-5t^2)D_n(t)dt$
with $D_0(x)=1$

These are polynomials of degree 3n. I mostly care about the leading term as n grows, but I'd also like to know how to compute these polynomials to high degree assuming only this recurrence. I have a horribly slow and ugly GP script which I am too embarrassed by to inflict on this list, and it doesn't go high enough to get the digits to guess at asymptotics.

More generally, I can't find from the usual sources how to accelerate computing families of *functions* given by recurrences. If for example all one had was this recurrence for Bernoulli polynomials, could one easily compute B_1000(x) (or at least an arbitrary coefficient) in PARI/GP?

As usual any help, comments or references are appreciated.