(50g) Bernoulli polynomials

[ttw]

I thought to try a slightly different approach to the Bernoulli Polynomials. These have several nice properties that make the easy to computer. I just keep the coefficients (as fractions in the code below, but these could be reals or numerator-denominator pairs) rather than the whole polynomial.

The procedure is based on the fact that each Bernoulli Polynomial has a simple integral (or differential from the other direction) relation with the next one in the series. The program operates in three steps. First the coefficients of the integrated polynomial are computed. Second these are normalized to have a leading coefficient of one. Third; the constant of integration is added such that the polynomial's integral is zero.

I use an auxiliary program that yields a list {1,2,3...N} when given N as an input. This is actually a program that I have written for every programmable calculator I have; it seems to occur rather commonly. I could just put this inline.

Program KSEQ: input N, output a list from 1 to N. (No check on illegal input.)

<< 1 SWAP FOR I I NEXT DUP ->LIST >>

The next program takes a Bernoulli Polynomial as input and updates it to the next higher polynomial. The flags are set for exact, non-umeric, and polynomials being 1+x... rather then ...x+1. Thereare three local variables (though these could be eliminated by clever stack usage; this didn't make things run faster though and was a bit harder to read.) Variable X is the list of coefficients; N is the length of the current polynomial (1 more than the index) and E is the index of powers of the polynomial.

<< DUP SIZE XQ 0 -> X N E <<
X N KSEQ DUP 'E' STO / N * EXPAND DUP E 1 ADD / {0} ƩLIST NEG SWAP + EXPAND >> >>