(HP71B) Bernoulli Polynomials

[Namir]

Bernoulli Polynomials for HP-17B
================================

Memory Map
==========

N: Order (Integer)
N0:Last N (Integer)
X: X argument (Floating point)
B(): Bernoulli number array (Floating point)
B0: Bernoulli polynomial Bern(N,X) (Floating point)
I, J, K, M: Loop control variables (Integer)


Usage
=====

1) Press RUN key.
2) The program displays the prompt "N? ".
3) Enter N and press the END LINE key.
4) The program displays the prompt "X? ".
5) Enter X and press the END LINE key.
6) The program displays N, X, and the value of the Bernoulli polynomial.
7) Press the [f][+] keys and go to step 2.

Example
=======

Examples of B(N,X) are:

B(8,2) = 7.96666666751
B(6,3) = 198.023809524
B(4,4) = 143.966666667

Listing
=======

Code:


Line     Step
10    REM BERNOULLI POLYNOMIAL
20    DESTROY ALL
30    N0=0
40    INPUT "N? ";N
50    INPUT X? ";X
60    IF N=N0 THEN GOTO 100
70    DIM B(N)
80    CALL BERN(B(),N)
90    N0=N
100    B0=0
110    FOR K=0 TO N
120    B0=B0+FACT(N)/FACT(K)/FACT(N-K)*B(K)*X^(N-K)
130    NEXT K
140    DISP "B(";N;",";X;"=";B0
150    PAUSE
160    GOTO 40
170    END
180    SUB BERN(B(),N)
190    B(0)=1 @ B(1)=-.5
200    FOR M=2 TO N
210    S=0.5-1/(M+1)
220    FOR K=2 TO M-1
230    R=1
240    FOR J=2 TO K
250    R=R*(J+M-K)/J
260    NEXT J
270    S=S-R*B(K)
280    NEXT K
290    B(M) = S
300    NEXT M
310    FOR M=3 TO N-1 STEP 2
320    B(M) = 0
330    NEXT M
340    END SUB

The above code checks if the value of N you enter is the same one you entered the last time. If so, the code skips calling the subroutine BERN since it already has the values of the Bernoulli numbers from the last round of calculations.

[C.Ret]

Hello,

Introduction:

Calculating Bernoulli polynomials is an historic operation where Bernoulli numbers are calculated mechanically as Ada Lovalace did on Babbage's Analytical Machine a long long time ago.
This mechanical machine only processed numbers and was very capable in computations. In comparison the HP-71B is a monster of power and ease of use.

The code of Namir inspire me a new version where display capacity of this powerful (infinitely powered compared to the historic Babbage's Analytical Machinery) Pocket Computer is used to express Bernoulli Polynomial as symbolic expression.

So I come here to share with you my adaptation of NAMIR's code which consists of three parts:

Initialization:
From the binary code stored in the line of DATA, lines 30 to 50 create a CHARSET so that the power exponents are easier to read. This set of character also save precious space in the reduced display of our HP-71B.

Symbolic determination and display of the Bernoulli polynomial
Lines 100 to 190 symbolically display the polynomial as soon as its coefficients are determined. For this the Bernoulli Numbers are compute one after the other using the technique derive from the Pascal's triangle.

As Bernoulli Numbers are rational numbers, so are polynomial coefficients. This is why the program uses 4 arrays; A() and B() are respectively numerators and denominators of the rational Bernoulli Numbers. In the same way P() and Q() are respectively numerators and denominators of the polynomial's coefficients.

A sub-routine SUB RED(n,d) makes it possible to reduce the fractions n/d and also to adapt the numerical representation in case of decimal values in order to lighten the computations and to avoid miscalculations by overflows. This subroutine is called in several places as the calculations proceed in order to keep as many exact rational numbers as possible.


Numerical evaluation as rational and approximation value:
At lines 200 to 240, the user can enter a value for x. The numerical evaluation is performed and the result given, as far as possible, as an irreducible fraction. The variable E and F are used to memorized the rational sub-results of the computation, using again the same subroutine RED(E,F)
From this fraction E/F, the approximate value of the result Bn(x) is given.

Usage:
Launch the program using the RUN key or by any other appropriate means.
The title "Bernoulli Polynomials" remains displayed for a few seconds while generating the characters of the CHARSET allowing the display of the exponents.

At the B_ prompt, the user must enter the degree of the polynomial of interest.
The HP-71B symbolically determines this polynomial and displays it as the monomials are determined.

Then the x=_ prompt appears to the right of the polynomial. Scroll the screen if necessary to display the entire expression of the polynomial. For very high powers, remember to adjust the DELAY of the display in order to have time to read or copy it as it is determined.

At the x=_ prompt the user can enter a value and get the value of the polynomial at that point. When possible, the exact rational value is displayed as well as the approximate numerical value.

The user can also enter four commands at the x=_ prompt:
Q or q to quit and exit the program
N or n to change the degree n of the polynomial
X or x to re-display the symbolic expression of the current polynomial
B or b to obtain the expression of Bernoulli Numbers up to a(n)/b(n).


Program listing :

   

Exemples of outputs:

   



P.S.: As I copy-paste the program in this message, I realize that I am using the COMB(n,k) instruction in several places in the program.
This instruction which calculates the number of combinations of k elements among a set of n elements comes from the JPC ROM.
If someone wants to use this program on an HP-71B without this ROM, they will have to adapt by defining a quick and dirty user function:
5 DEF FNC0(N,K)=FACT(N)/FACT(K)/FACT(N-K)
and replacing every use of COMB(n,k) with FNC0(n,k) in the code.
Of course a better way to compute the COMB(n,k) may be used here, see following posts.

For example, a more efficient and stable may be:
500 DEF FNC0(N0,K0)
510 C0=1 @ K0=MIN(K0,N0-K0) @ IF K0=0 THEN 530
520 FOR I0=1 TO K0 @ C0=C0*N0/I0 @ N0=N0-1 @ NEXT I0
530 FNC0=C0 @ END DEF

A few words about the cryptic DATA on line 1:
A more readable version of this line would be the following. But unfortunately this line is too long to fit into only one entry.

1 DATA CHR$(31)&CHR$(17)&CHR$(31), CHR$(0)&CHR$(31), CHR$(29)&CHR$(21)&CHR$(23), CHR$(21)&CHR$(21)&CHR$(31), CHR$(4)&CHR$(7)&CHR$(4), CHR$(23)&CHR$(21)&CHR$(29), CHR$(31)&CHR$(20)&CHR$(28), CHR$(1)&CHR$(29)&CHR$(3), CHR$(31)&CHR$(21)&CHR$(31), CHR$(23)&CHR$(21)&CHR$(31)

The data are coding the pixels grid to draw each user defined character. Each exponents only use a fraction of the grid in order to have up to two exponent figures per character cell. That the aim of the second loop in lines 30-50 which combine pairs of figures to count from 0 up to 29. Each user defined character have to use exactly six bytes, padding CHR$(0) are insert to fit size.

Here an example of a character coding for ^7:
■ ■ ■ ▫ ▫ ▫ 1
□ □ ■ ▫ ▫ ▫ 2
□ ■ □ ▫ ▫ ▫ 4
□ ■ □ ▫ ▫ ▫ 8
□ ■ □ ▫ ▫ ▫ 16
□ □ □ ▫ ▫ ▫ 32
□ □ □ ▫ ▫ ▫ 64
□ □ □ ▫ ▫ ▫ 128
1 2 3 . . .
. 9 . . . .
CHR$(1) g CTRL A "°"
& CHR$(29) g CTRL g ] "≠"
& CHR$(3) g CTRL C "←"

Fortunately, main of the CHR$(nnn) character are directly available through normal keystroke on the HP-71B mainly using g CTRL sequences.
For example CHR$(29) is display on HP-71B screen as ≠and can be type with the sequence g CTRL g ].
The CHR$(3) display an ← and is available through g CTRL C sequence.

There is a few exceptions, such as CHR$(31) (string "▒") which need to affect a dedicated user-key since no sequence exist. I used DEF KEY "V",CHR$(31); such that I may enter the ▒ character by typing g 1USER V.

Here is a table summarised the codes used for each single exponent with string représentations.

Code:

^0: CHR$(31)&CHR$(17)&CHR$(31) "▒Ω▒"     ^1: CHR$( 0)&CHR$(31)        CHR$(0)&"▒"        
^2: CHR$(29)&CHR$(21)&CHR$(23) "≠ÄÖ"     ^3: CHR$(21)&CHR$(21)&CHR$(31) "ÄÄ▒"
^4: CHR$( 4)&CHR$( 7)&CHR$(30) "♠α£"     ^5: CHR$(23)&CHR$(21)&CHR$(29) "ÖÄ≠"
^6: CHR$(31)&CHR$(20)&CHR$(28) "▒πΣ"     ^7: CHR$( 1)&CHR$(29)&CHR$( 3) "°≠←"
^8: CHR$(31)&CHR$(21)&CHR$(31) "▒Ä▒"     ^9: CHR$(23)&CHR$(21)&CHR$(31) "ÖÄ▒"

To type all these data into the HP-71C, one may prefer enter several DATA lines:

1 DATA CHR$(31)&CHR$(17)&CHR$(31), CHR$(0)&CHR$(31), CHR$(29)&CHR$(21)&CHR$(23), CHR$(21)&CHR$(21)&CHR$(31)
2 DATA CHR$(4)&CHR$(7)&CHR$(30), CHR$(23)&CHR$(21)&CHR$(29), CHR$(31)&CHR$(20)&CHR$(28)
3 DATA CHR$(1)&CHR$(29)&CHR$(3), CHR$(31)&CHR$(21)&CHR$(31), CHR$(23)&CHR$(21)&CHR$(31)



EDITED:
Correct Ada Lovegrace forename.
Add a link to the useful comment by John Keith from concerning the way of computing FN C0(n,k)
Remove any criticism about Namir's code which I find great and that have me inspire with great excitement this new version. The aim of this post was not to offense any one. I was rushed and very awkward.

[rprosperi]

(07-17-2023 07:16 PM)C.Ret Wrote:  [snip]

P.S.: As I copy-paste the program in this message, I realize that I am using the COMB(n,k) instruction in several places in the program.
This instruction which calculates the number of combinations of k elements among a set of n elements comes from the PPC ROM.
If someone wants to use this program on an HP-71B without this ROM, they will have to adapt by defining a user function:
5 DEF FNC0(N,K)=FACT(N)/FACT(K)/FACT(N-K)
and replacing every use of COMB(n,k) with FNC0(n,k) in the code.

I'm pretty sure you mean the JPC ROM, which is available here:

http://www.jeffcalc.hp41.eu/emu71/jpcrom.html

[John Keith]

Quite an impressive program. Congratulations for representing Bernoulli numbers as rational numbers on a system that does not support rationals. However, ...

(07-17-2023 07:16 PM)C.Ret Wrote:  If someone wants to use this program on an HP-71B without this ROM, they will have to adapt by defining a user function:
5 DEF FNC0(N,K)=FACT(N)/FACT(K)/FACT(N-K)
and replacing every use of COMB(n,k) with FNC0(n,k) in the code.

Computing binomial coefficients this way is numerically unstable and will likely undo the advantages of maintaining Bernoulli numbers as integer ratios. This section of the Wikipedia page on binomial coefficients gives several useful methods.

Finally, not to nit-pick, but it's Ada, not Eva.