(50g) Bernoulli polynomials

12172018, 05:31 PM
(This post was last modified: 12172018 10:10 PM by peacecalc.)
Post: #1




(50g) Bernoulli polynomials
Hello HP50gfans,
I made a program for Bernoullipolynomials ( = BP), it works with three equations: \[ B_o(x) = 1 ~~~~~~~ \text{start}\] \[ B_n(x) = n \cdot \int B_{n1}(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:
Enjoy it and tell me your improvements and critics. SORRY SORRY SORRY, dear Mr. Admin move this thread to "General Software Libary". greetings peacecalc 

12172018, 06:07 PM
(This post was last modified: 12172018 06:29 PM by DavidM.)
Post: #2




RE: (50g) Bernoulli polynomials
(12172018 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. 

12172018, 07:03 PM
Post: #3




RE: (50g) Bernoulli polynomials
(12172018 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 second this method. One can also begin the program with Code:
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. 

12172018, 08:37 PM
Post: #4




RE: (50g) Bernoulli polynomials
(12172018 07:03 PM)John Keith Wrote: One can also begin the program with That's even better, especially if the single argument use case is more common than the two argument one. Nice job! 

12172018, 09:17 PM
Post: #5




RE: (50g) Bernoulli polynomials
There are some direct formulas without integrations (although I always use integration because it's nonnumerical 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. 

12172018, 10:03 PM
Post: #6




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 

12172018, 11:58 PM
(This post was last modified: 12172018 11:58 PM by ijabbott.)
Post: #7




RE: (50g) Bernoulli polynomials
(12172018 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!(nk)!). — Ian Abbott 

12222018, 11:06 PM
Post: #8




RE: (50g) Bernoulli polynomials
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 psifunction appears? I'm not shure, but is there a bug known for the summation command? 

12232018, 07:58 PM
Post: #9




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. 

12292018, 02:18 PM
Post: #10




RE: (50g) Bernoulli polynomials
Hello Keith,
of course here is the code: Code:
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:
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:
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 

12292018, 03:52 PM
Post: #11




RE: (50g) Bernoulli polynomials
(12292018 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 ... It is the summation, more specifically, the term 1/12*(55*9^11) 55*9^11 = 1 725958 278495 > roundingup > 1 725958 278500 Introduced constant term error = 5/12 ~ 0.41667 

12292018, 04:40 PM
Post: #12




RE: (50g) Bernoulli polynomials
(12292018 03:52 PM)Albert Chan Wrote: It is the summation, more specifically, the term 1/12*(55*9^11) 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. 

12292018, 05:54 PM
Post: #13




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


12302018, 01:44 PM
Post: #14




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:
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 CASfunctions "INTVX", "ISOL" and "SUBTL")... The summations goes quatric with n... greetings peacecalc 

01022019, 11:13 PM
Post: #15




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 numeratordenominator 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, nonumeric, 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 >> >> 

01032019, 10:20 AM
Post: #16




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. Thanx in advance peacecalc 

« Next Oldest  Next Newest »

User(s) browsing this thread: