(50g) Bernoulli polynomials
12-17-2018, 05:31 PM (This post was last modified: 12-17-2018 10:10 PM by peacecalc.)
Post: #1
 peacecalc Member Posts: 187 Joined: Dec 2013
(50g) Bernoulli polynomials
Hello HP50g-fans,

I made a program for Bernoulli-polynomials ( = BP), it works with three equations:

$B_o(x) = 1 ~~~~~~~ \text{start}$

$B_n(x) = n \cdot \int B_{n-1}(t)dt ~~~~~~~ \text{recursive definition}$

$\int_0^1 B_n(t)dt = 0 ~~~~~~~~ \text{ calculate the constant from second equation}$

Therefore we have polynomials it is easy to program and make use of CAS of the HP50g.
The program "BCALC" has to possibilities to work:

a) with one parameter "N", then it will calculate all BP from 0 to N, but it only stores the last N'th BP.

b) with two parameters "K" to "N". Let us say you calculated BP for N = 5, then you decide to calculate further on to N = 7. If you don't delete the global variables "BN1" and "BN2", you can input 6 and 7 and you get BP for N=7 from N = 5 as first equation.
That is a advantage in time, because the CAS operations are very slow.

To let the program know how many parameters I input, I use a list with one or two elements.

The program uses only one integration because the second one (stored in "BN2") prepares the BP for "N+1" (again time saving). I only make comments for the first part, I hope the rest of the program is understandable for you.

Code:
 %%HP: T(2)A(R)F(.); « OBJ->     IF 2. ==   THEN \-> K N     « K N FOR J              BN2 J * 'C' + EVAL DUP 'BN1' STO   @@RCL the prepared BP                                                              @@ and multiplied with J                                                              @@ add the constant "C" and stored                                                             @@in BN1 temporally             INTVX EVAL DUP 'BN2' STO            @@ integrated BP and stored                                                             @@in BN2 temporally             'X=1' SUBST EVAL 0 = 'C' ISOL DUP @@ set the upper limit 1 and calculate                                                             @@ the constant "C"             BN1 SWAP SUBST EVAL FDISTRIB   @@ RCL the BN1 und substitute the "C"             'BN1' STO                                   @@ with its number and store it in BN1             BN2 SWAP SUBST EVAL 'BN2' STO   @@ RCL the BN2 und substitute the "C"                                                             @@ with its number and store it in BN2                                                             @@ as the prepared BP of J + 1             NEXT     »   ELSE \-> N     « B0 INTVX 'C' + DUP 'BN1' STO        INTVX DUP 'BN2' STO        'X=1' SUBST EVAL 0 = 'C' ISOL DUP       BN1 SWAP SUBST EVAL 'BN1' STO        BN2 SWAP SUBST EVAL 'BN2' STO       IF N 1 >        THEN            2 N FOR J                 BN2 J * 'C' + EVAL DUP                 'BN1' STO INTVX EVAL DUP                 'BN2' STO 'X=1' SUBST EVAL 0 = 'C' ISOL DUP                 BN1 SWAP SUBST EVAL FDISTRIB                 'BN1' STO                 BN2 SWAP SUBST EVAL                 'BN2' STO         NEXT       END     »   END »

Enjoy it and tell me your improvements and critics.

SORRY SORRY SORRY, dear Mr. Admin move this thread to "General Software Libary".

greetings
peacecalc
12-17-2018, 06:07 PM (This post was last modified: 12-17-2018 06:29 PM by DavidM.)
Post: #2
 DavidM Senior Member Posts: 848 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
(12-17-2018 05:31 PM)peacecalc Wrote:  Enjoy it and tell me your improvements and critics.

Hi peacecalc!

I won't pretend to know anything about the math side of what this does (or is appropriate for), but reading through your description sparked a thought about a possible alternative means to deal with a varying number of stack arguments.

In my RPL programs, I've always tried to avoid situations where the stack can't have results from previous operations on it, since I feel like anything on the stack above the needed arguments for a given program is purely up to the user -- why limit the user's ability to have the results from other calculations present? One of the best things about RPN operations (IMHO) is the ability to retain the results of intermediate calculations while you are experimenting. Forcing the user to have to clear the stack before executing a program can be very disruptive.

One very easy way to accommodate the possibility of a varying number of arguments is to place them in a list { } before invoking the program. Placing a list of items on the stack is actually very easy (only one shifted key more than individually loading them). Your program then simply needs to execute LIST→ instead of DEPTH to determine how many arguments were given.

I believe this is a small adjustment to make in comparison to the benefits to the user of not having to clear the stack first.

Just a thought!

Edit:

One other possibility... since you are looking for either 1 or 2 arguments, simply mandate that the user place both of them on the stack, using ENTER to create a second copy if only one is actually needed. The program could then check to see if the two stack arguments are equal before deciding which course of action to take. That's just one more keystroke (ENTER) in the case of only one argument.
12-17-2018, 07:03 PM
Post: #3 John Keith Senior Member Posts: 751 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
(12-17-2018 06:07 PM)DavidM Wrote:  One very easy way to accommodate the possibility of a varying number of arguments is to place them in a list { } before invoking the program. Placing a list of items on the stack is actually very easy (only one shifted key more than individually loading them). Your program then simply needs to execute LIST→ instead of DEPTH to determine how many arguments were given.

I believe this is a small adjustment to make in comparison to the benefits to the user of not having to clear the stack first.

I second this method. One can also begin the program with

Code:
 IF DUP TYPE 5. SAME THEN OBJ\->... (do whatever with multiple args) ELSE... (use program-supplied defaults)

to allow either a single argument with the program supplying defaults, or some or all arguments supplied by the user.

For example compare UTPN on the HP 50 with NORMALD_CDF, its rough equivalent on the Prime.
12-17-2018, 08:37 PM
Post: #4
 DavidM Senior Member Posts: 848 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
(12-17-2018 07:03 PM)John Keith Wrote:  One can also begin the program with

Code:
 IF DUP TYPE 5. SAME THEN OBJ\->... (do whatever with multiple args) ELSE... (use program-supplied defaults)

to allow either a single argument with the program supplying defaults, or some or all arguments supplied by the user.

That's even better, especially if the single argument use case is more common than the two argument one. Nice job!
12-17-2018, 09:17 PM
Post: #5
 ttw Member Posts: 248 Joined: Jun 2014
RE: (50g) Bernoulli polynomials
There are some direct formulas without integrations (although I always use integration because it's non-numerical and easier to remember.)

B_m(X) = Sum(n:0:m)[1/(n+1)Sum(k:0:n)[(-1)^k*C(n,k)*(x+k)^m]

It's easier to start with 1 and integrate then scale the leading coefficient to 1 and add the constant to get the average to be zero. That's probably equivalent to what you are doing. It can all be done symbolically.

It's useful for periodizing functions before doing some types of integrations. The results are much "smoother" (fewer high frequency left over terms) than most periodization methods.
12-17-2018, 10:03 PM
Post: #6
 peacecalc Member Posts: 187 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Hello all,

wow, what a resonance! Thank you all!

@DavidM,@JohnKeith: Yes, the input is (better now: was) quit a problem. I like to change the code for inputting a list.

@ttw, I've never seen such a double sum for the BP's. What does C(n,k) mean? It is the Stirling numbers or a kind of permutations?

greetings
peaceglue
12-17-2018, 11:58 PM (This post was last modified: 12-17-2018 11:58 PM by ijabbott.)
Post: #7 ijabbott Senior Member Posts: 1,090 Joined: Jul 2015
RE: (50g) Bernoulli polynomials
(12-17-2018 10:03 PM)peacecalc Wrote:  What does C(n,k) mean? It is the Stirling numbers or a kind of permutations?

Combinations, binomial coefficient, nCk, n!/(k!(n-k)!).

— Ian Abbott
12-22-2018, 11:06 PM
Post: #8
 peacecalc Member Posts: 187 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Hello ijabbott,

In general I tried to write a program for the hp50g, that uses the formular from ttw and another which operates with three sums, but the results are strange...

The calc needs the complex mode and and for higher B_n(x) (n >= 11) I get wrong results (for n is odd B_n(0) = Bernoulli Number has to be zero, but it isn't), and more worse: the resulting coefficients are outputed as floating point numbers, even when the input are integers (and same in program)? And sometime the psi-function appears? I'm not shure, but is there a bug known for the summation command?
12-23-2018, 07:58 PM
Post: #9 John Keith Senior Member Posts: 751 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Show what you are inputting that is giving you the strange results. Also calculator modes (exact vs approx etc.).

If you are getting Psi as part of the result then it is returning symbolic results.
12-29-2018, 02:18 PM
Post: #10
 peacecalc Member Posts: 187 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Hello Keith,

of course here is the code:
Code:
 « \-> M X '\GS(N=0,M,1/(N+1)*\GS(K=0,N,(-1)^K*COMB(N,K)*(X+K)^M))' »

That is ttw's expression: when I input 11 and 'X', at the machine forced me to switch to approximative mode (If I refuse to switch in approximative mode, the program interrupts and ends...) I get after some minutes:

Code:
 'X^11.+-5.5*X^10.+9.1666666666*X^9.+0.*X^8.+-11.*X^7.+-.00000002*X^6.+10.99999999*X^5.+-.000003*X^4.+-5.500011*X^3.+-.416696676667*X^2.+.416646633333*X-.41678'

Most of the constants in front of the powers of X are approx. correct, but not the last one: .41678, that one have to be 0 or near by, because every Bernoulli number with odd n greater or equal 3 is zero. And the approximative numbers of X^4 or X^6 are drifting away from zero (it should be like the zero by X^8)

My program working with the integral returns:

Code:
 'X^11-11*X^10/2+55*X^9/6-11*X^7+11*X^5-11*X^3/2+5*X/6'

The coeffcients are exact, because the machine works in exact mode.
So what is the problem: the summation function \GS or the combination function COMB?

After some work on the programs workig with summations the PSI function dont't appear any longer, but don't ask me why.

greetings
peacecalc
12-29-2018, 03:52 PM
Post: #11
 Albert Chan Senior Member Posts: 1,791 Joined: Jul 2018
RE: (50g) Bernoulli polynomials
(12-29-2018 02:18 PM)peacecalc Wrote:  Most of the constants in front of the powers of X are approx. correct, but not the last one: .41678 ...
So what is the problem: the summation function \GS or the combination function COMB?

It is the summation, more specifically, the term 1/12*(-55*9^11)

55*9^11 = 1 725958 278495 --> rounding-up --> 1 725958 278500

Introduced constant term error = -5/12 ~ -0.41667
12-29-2018, 04:40 PM
Post: #12
 DavidM Senior Member Posts: 848 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
(12-29-2018 03:52 PM)Albert Chan Wrote:  It is the summation, more specifically, the term 1/12*(-55*9^11)

55*9^11 = 1 725958 278495 --> rounding-up --> 1 725958 278500

Introduced constant term error = -5/12 ~ -0.41667

Not long ago I discovered that Σ will always attempt to apply a special form of CAS simplification to its argument; perhaps the method it uses is causing this scenario to unfold.
12-29-2018, 05:54 PM
Post: #13 John Keith Senior Member Posts: 751 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Indeed, I discovered this while comparing the speed of various programs to calculate exact (rational) values of harmonic numbers. Using the \GS (Sigma) function almost immediately returns a symbolic result like 'Psi(n) - Psi(1)'.
12-30-2018, 01:44 PM
Post: #14
 peacecalc Member Posts: 187 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Hello folks,

now I have the summation function \GS replaced by FOR NEXT loops:

Input is n and the variable for the polynomials 'X'

Code:
 « 0 DUP \-> M X RSI RSO   « 0 M     FOR N 1 N + INV 0 N 0 'RSI' STO       FOR K -1 K ^ N K COMB * X K + M ^ * 'RSI' STO+       NEXT RSI * 'RSO' STO+     NEXT RSO EVAL FDISTRIB    » »

For getting a symbolic calculation it is necessary to use integer values for the loops (although floatings has a faster performance).

Nice comparison, the "integral" code is slower until n = 18 or 19 and the "sum" code is faster. For n greater or equal 20 the tide turned, the "integral" code is faster (although it uses the CAS-functions "INTVX", "ISOL" and "SUBTL")...
The summations goes quatric with n...

greetings
peacecalc
01-02-2019, 11:13 PM
Post: #15
 ttw Member Posts: 248 Joined: Jun 2014
RE: (50g) Bernoulli polynomials
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 >> >>
01-03-2019, 10:20 AM
Post: #16
 peacecalc Member Posts: 187 Joined: Dec 2013
RE: (50g) Bernoulli polynomials
Hello ttw,

I keyed in your proggies, but it doesn't work. Can you help me? I tried to take as a input
the B_1(x) = X - 1/2, but I don't get the expression for B_2(x). I got a error message:
ƩLIST Error: Invalid dimension.