(49g 50g) Bernoulli, Euler and Bessel Polynomials

[John Keith]

This group of programs computes the values of the Bernoulli polynomials, Euler polynomials and reverse Bessel polynomials. Given an integer n on level 2 and a number x on level 1, the program BERPOLY returns the value of the nth Bernoulli polynomial at x. The program EUPOLY similarly computes the value of the nth Euler polynomial at x.

The program BESSEL similarly returns the value of the nth reverse Bessel polynomial at x. If regular Bessel polynomials are desired, insert REVLIST after the / on line 5 of the program.

The programs require ListExt 1.3

For both programs, x can be an approximate (floating-point) number, an exact integer, a rational number or a symbolic expression (e. g. 'X' ). If the level 1 argument is 'X', the programs will return the polynomial as an expression, similar to the built-in polynomial commands. The programs should be used in exact mode if exact integer or rational results are desired.

The two remaining programs are subprograms called by the main polynomial programs. PASCROW is the first program in this post. It returns the nth row of Pascal's triangle given n on the stack. EUNUM returns terms 0 through n of A155585, the signed Euler tangent numbers.

BERPOLY and EUPOLY have been checked against the Mathematica functions BernoulliB and EulerE respectively. Both return the same results for all values that I have checked, although the form of symbolic expressions may be different.

Code:


DIR
  BERPOLY
  \<< SWAP DUP I\->R 1.
  IF >                        @ Main body of program if n > 1
  THEN DUP PASCROW            @ Row n of Pascal's triangle
    0 ROT LSEQR IBERNOULLI *  @ Multiply by Bernoulli numbers 0..n
    { EVAL EVAL } LMAP        @ Fully simplify rational numbers
    AXL SWAP PEVAL            @ Evaluate polynomial at x
  ELSE I\->R 1. SAME
    { 2 INV - }               @ If n = 1 then return x-1/2
    { DROP 1 } IFTE           @ If n = 0 then return 1
  END EVAL                    @ Simplify result
\>>
  EUPOLY
  \<< SWAP DUP I\->R 1.
  IF >                         @ Main body of program if n > 1
  THEN DUP EUNUM               @ Signed Euler tangent numbers
    OVER PASCROW *             @ Multiply by row n of Pascal's triangle
    SWAP 1 -2 ROT 1 + LMSEQ /  @ Divide by powers of -2
    AXL SWAP PEVAL             @ Evaluate polynomial at x
  ELSE I\->R 1. SAME
    { 2 INV - }                @ If n = 1 then return x-1/2
    { DROP 1 } IFTE            @ If n = 0 then return 1
  END EVAL                     @ Simplify result
\>>
  BESSEL
  \<< SWAP DUP I\->R 1.
  IF >
  THEN DUP PASCROW                              @ Row n of Pascal's triangle
    OVER 1 + DUPDUP 2 - + LSEQR 1 :: * LSCAN *  @ List of rising factorials
    SWAP 1 2 ROT 1 + LMSEQ /                    @ Divide by powers of 2
    AXL SWAP PEVAL                              @ Evaluate polynomial at x
  ELSE I\->R 1. SAME
    { 1 + }                                     @ If n = 1 then return x+1
    { DROP 1 } IFTE                             @ If n = 0 then return 1
  END EVAL                                      @ Simplify result
\>>
  PASCROW
  \<< DUP I\->R
    IF DUP 1. >
    THEN
      IF DUP 3. >
      THEN DUP 2. / IP SWAP 1. + 2. / IP  @ Indices for left half of row
      ELSE 0.
      END \-> m t
      \<< 1 SWAP LSEQ                     @ List of 1..n
        DUP 1. m SUB SWAP REVLIST 1. m SUB 2.
        \<< PICK3 * SWAP /
        \>> DOLIST +                      @ ConvOffs transform
        DUP 1. t SUB REVLIST +            @ Reverse and append right half
      \>>
    ELSE 1. SAME { DUP 2. \->LIST } { DROP { 1 } } IFTE
    END
\>>
  EUNUM
  \<< I\->R \-> n
  \<< 1 { 1 } 1. n                                   @ Save a(0), set up temp. list
    FOR k REVLIST 0 :: + LSCAN                       @ Accumulate list with 0
      DUP k 2. MOD { LPOPR NIP } :: HEAD IFTE NEWOB  @ Pick either first or last term
      k 4. MOD 3. SAME :: NEG IFT SWAP               @ Negate every 4th term and save
    NEXT DROP n 1. + \->LIST                         @ Drop temp. list, create output list
  \>>
\>>
END

[John Keith]

Updated previous post and thread title to add a program for Bessel polynomials, and to clean up the code for BERPOLY and EUPOLY.