(HP71B) Bernoulli Polynomials

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