(48G 49 50g) Stirling and Stirling-Bernoulli Transforms

[John Keith]

First, the Stirling transform, also known as the Stirling S2 transform. More information here and here.

Code:


\<< DUP SIZE \-> a n
  \<< a HEAD { 1 } 2. n
    FOR k 0 SWAP + 2.
      \<< NSUB R\->I * +
      \>> DOSUBS 1 + a 1. k SUB OVER * \GSLIST SWAP
    NEXT DROP n \->LIST
  \>>
\>>

Next the inverse Stirling transform, also known as the Stirling S1 transform:

Code:


\<< DUP SIZE \-> b n
  \<< b HEAD { 1 } 1 n 1. - R\->I
    FOR k 0 SWAP + 2.
      \<< k * -
      \>> DOSUBS 1 + b 1 k 1 + SUB OVER * \GSLIST SWAP
    NEXT DROP n \->LIST
  \>>
\>>

Next, the Stirling-Bernoulli transform. While this transform is defined in terms of products of Stirling numbers and factorials, the program uses an algorithm described in this paper which is about three times as fast.

Code:


\<< DUP SIZE \-> n
  \<< DUP HEAD SWAP 2. n
    START 2.
      \<< - NSUB R\->I * EVAL
      \>> DOSUBS DUP HEAD NEWOB SWAP
    NEXT DROP n \->LIST
  \>>
\>>

Finally, the inverse Stirling-Bernoulli transform, similar to the program above.

Code:


\<< DUP SIZE 1. - R\->I \-> n
  \<< DUP HEAD SWAP 1 n
    FOR k 2.
      \<< k / - EVAL
      \>> DOSUBS DUP HEAD NEWOB SWAP
    NEXT DROP n 1 + \->LIST
  \>>
\>>

The EVAL commands inside the program arguments for DOSUBS are required because the terms of sequences returned by these programs may not be integers.

All of these programs will work on the 48G as long as the I->R commands are removed. The EVAL commands and decimal points following integers can also be removed. However, the results may not be accurate for some sequences due to the 12-digit limit of the 48G.

Edited 9/20/2019 to add NEWOB to both Stirling-Bernoulli transform programs. This speeds execution of the programs and allows them to process longer lists without running out of memory.