(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[John Keith]

(08-30-2023 11:59 PM)Albert Chan Wrote:  

(08-30-2023 04:30 PM)Albert Chan Wrote:  B(6) = [1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * [-1!, 2!, -3!, 4!, -5!, 6!] = 1/42

We could use horner's rule, and remove factorial function.

B(6) = -(1/2 - 2*(31/3 - 3*(90/4 - 4*(65/5 - 5*(15/6 - 6*(1/7))))))

Nice use of Horner's method.

The expression that you are computing, (-1)^k*k!*Stirling2(n+1, k+1), is A163626 which is actually simpler to compute than the Stirling numbers of the second kind. From the linked page, T(n, k) = (k+1)*T(n-1,k) - k*T(n-1,k-1). Thus each term in your expression can be computed with just one subtraction and two multiplications, after which each term must be divided by (k+1) as above.

Here is my RPL implementation, which computes the whole triangle but can be easily modified to do one row at a time. \GDLIST is deltaList and LSEQ returns a list of integers 1..n, equivalent to range(1, n+1).

Code:


\<< I\->R \-> n
  \<< { 1 } n R\->I LSEQ 1. n
    FOR k DUP2 1. k SUB * 0 + \GDLIST 1 SWAP + SWAP
    NEXT DROP n 1. + \->LIST
  \>>
\>>

Edit: modified version, returns row n.

Code:


\<< I\->R \-> n
  \<< { 1 } n R\->I LSEQ 1. n               @ Row 0, list of 1..n
    FOR k SWAP OVER 1. k SUB * 0 + \GDLIST  @ (k+1)*T(n-1,k) - k*T(n-1,k-1)
      1 SWAP + SWAP                         @ Append last term = 1
    NEXT DROP                               @ Drop list of 1..n
  \>>
\>>