(50g) Bernoulli polynomials

[peacecalc]

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

[DavidM]

(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.

[John Keith]

(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.

[DavidM]

(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!

[ttw]

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.

[peacecalc]

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

[ijabbott]

(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)!).

[peacecalc]

Hello ijabbott,

thank you for answer!

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?

[John Keith]

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.

[peacecalc]

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

[Albert Chan]

(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

[DavidM]

(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.

[John Keith]

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)'.

[peacecalc]

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

[ttw]

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 >> >>

[peacecalc]

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.

Thanx in advance
peacecalc