+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (49g 50g) Bernoulli, Euler and Bessel Polynomials (/thread-17824.html)
(49g 50g) Bernoulli, Euler and Bessel Polynomials - John Keith - 12-17-2021 09:35 PM
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
\>>
\>>
ENDRE: (49g 50g) Bernoulli, Euler and Bessel Polynomials - John Keith - 12-26-2021 09:01 PM
Updated previous post and thread title to add a program for Bessel polynomials, and to clean up the code for BERPOLY and EUPOLY.