(HP15C)(HP67)(HP41C) Bernoulli Polynomials
08-28-2023, 01:10 PM (This post was last modified: 08-28-2023 10:28 PM by Namir.)
Post: #1
 Namir Senior Member Posts: 1,083 Joined: Dec 2013
(HP15C)(HP67)(HP41C) Bernoulli Polynomials
Program calculates Bernoulli polynomials using nested summations (see Wikipedia).

The memory map is:
Code:
 R0 = x R1 = m R2 = n R3 = k R4 = Bm(x) R5 = inner sum

The listing for HP15C is:
Code:
 01        LBL A         02        STO 0        # store x 03        X<>Y          04        STO 1        # store m 05        0         06        STO 2        # set n = 0 07        STO 4        # set Bm(x) = 0 08        LBL 1        # start of outer summation 09        0         10        STO 3        # set k=0 11        STO 5        # set inner sum = 0 12        LBL 2        # start of inner summation 13        RCL 2         14        RCL 3         15        Cy,x         16        RCL 0         17        RCL 3         18        +         19        RCL 1         20        Y^X         21        *         22        1         23        CHS         24        RCL 3         25        Y^X         26        *         27        STO+ 5        # update inner sum 28        1         29        STO+ 3         30        RCL 2         31        RCL 3         32        X<=Y?        # test end of inner loop 33        GTO 2         34        RCL 5         35        RCL 2         36        1         37        +         38        /         39        STO+ 4        # Update Bm(x) 40        1         41        STO+ 2         42        RCL 1         43        RCL 2         44        X<=Y?        # test end of outer loop 45        GTO 1         46        RCL 4         47        RTN

The listing for HP67 is:
Code:
 01        LBL A 02        STO 0 03        X<>Y  04        STO 1 05        0 06        STO 2 07        STO 4 08        LBL 1 09        0 10        STO 3 11        STO 5 12        LBL 2 13        RCL 2 14        N! 15        RCL 3 16        N! 17        / 18        RCL 2 19        RCL 3 20        - 21        N! 22        / 23        RCL 0 24        RCL 3 25        + 26        RCL 1 27        Y^X 28        * 29        1 30        CHS 31        RCL 3 32        Y^X 33        * 34        STO+ 5 35        1 36        STO+ 3 37        RCL 2 38        RCL 3 39        X<=Y? 40        GTO 2 41        RCL 5 42        RCL 2 43        1 44        + 45        / 46        STO+ 4 47        1 48        STO+ 2 49        RCL 1 50        RCL 2 51        X<=Y? 52        GTO 1 53        RCL 4 54        RTN

Listing for HP-41C is:

Code:
 01        LBL "BERNL" 02        LBL A 03        STO 00 04        X<>Y  05        STO 01 06        0 07        STO 02 08        STO 04 09        LBL 01 10        0 11        STO 03 12        STO 05 13        LBL 02 14        RCL 02 15        FACT 16        RCL 03 17        FACT 18        / 19        RCL 02 20        RCL 03 21        - 22        FACT 23        / 24        RCL 00 25        RCL 03 26        + 27        RCL 01 28        Y^X 29        * 30        1 31        CHS 32        RCL 03 33        Y^X 34        * 35        STO+ 05 36        1 37        STO+ 03 38        RCL 02 39        RCL 03 40        X<=Y? 41        GTO 02 42        RCL 05 43        RCL 02 44        1 45        + 46        / 47        STO+ 04 48        1 49        STO+ 02 50        RCL 01 51        RCL 02 52        X<=Y? 53        GTO 01 54        RCL 04 55        RTN

Example 1

5
ENTER
3.5
f A

output is 220.9375

Example 2

3
ENTER
5.5
f A

output is 123.7500

Note: I was thinking about writing a second version of the listing for the HP-14C that uses the ISG command for registers 2 and 3. The problem is whenever the listing uses RCL 2 and RCL 3 it must be followed by the INT command. This adds more steps IMHO. So I abandoned the ISG-using version.
08-29-2023, 01:57 AM
Post: #2
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
Code:
B(m,x) := sum(1/(n+1) * sum((-1)^k * comb(n,k) * (x+k)^m, k = 0 .. n), n = 0 .. m)

Another way to get Bm(x) is by synthetic division, for falling factorial form.

We can think of B as "derivative" of sum of power function, symbol Bm(x) == Bm

Σ(k^m, k=0 .. x-1) = ∫(Bm) = Bm+1 / (m+1)       + C

(08-05-2019 12:21 AM)Albert Chan Wrote:  From problem 12, Σ(xm, x=0 to n-1) = nm+1 / (m+1)

Except for symbol names, the 2 sums look very similar!
We can easily get x^m to its falling factorial form, via synthetic division.

The only issue is to get the integration constant, Bm(0) = B(m)
B(1) = -1/2, B(2k+1) = 0, we only need to worry about B(2k)

(07-28-2019 12:02 AM)Albert Chan Wrote:  Sum of k^6 formula = $$1\binom{n}{1}+63\binom{n}{2}+602\binom{n}{3}+2100\binom{n}{4}+3360\binom{n}{5}​+2520\binom{n}{6}+720\binom{n}{7}$$

B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42

Example, lets do formula for B6(x)

Synthetic Division, x^m to falling factorial form, for Bm formula
Code:
     1   0   0   0   0   0   0  // x^6 1>  1   1   1   1   1   1   1 2>  1   3   7  15  31  63 3>  1   6  25  90 301 4>  1  10  65 350 5>  1  15 140 6>  1  21

Constant of integration, adjusted for using falling factorial numbers.

B(6) = 1*0!/1 - 63*1!/2 + 301*2!/3 - 350*3!/4 + 140*4!/5 - 21*5!/6 + 1*6!/7 = 1/42

x^5 = x^4 * x, we pick numbers 3rd entries from the right.

x^5 = [1, 15, 25, 10, 1] • [x^1, x^2, x^3, x^4, x^5]

(B^6-B(6))/6 = [1/2, 15/3, 25/4, 10/5, 1/6] • [x^2, x^3, x^4, x^5, x^6]      // "integrate"

B^6 = 1/42 + [3, 30, 75/2, 12, 1] • [x^2, x^3, x^4, x^5, x^6]

Or, in efficient "horner" form:

B6(x) = 1/42 + x*(x-1)*(3 + (x-2)*(30 + (x-3)*(75/2 + (x-4)*(12 + (x-5)))))
08-29-2023, 11:47 AM
Post: #3
 Namir Senior Member Posts: 1,083 Joined: Dec 2013
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
A few weeks ago I published here a version of the HP-17B that evaluates Bernoulli polynomials. The BASIC code had a subroutine that calculates Bernoulli numbers and stores then in an array. The program allowed the user to repeatedly enter values for the polynomial order and for x. It remembers that last order used. If that is case, then it simply uses the current values in the Bernoulii numbers array.

I started thinking about using the same approach but realized the storage scheme is a bit quirky (and limited the order of the Bernoulli polynomials) since odd orders (except for 1) of the Bernoulli numbers are zero. When I saw the nested summation formula in Wikipedia, I decided that it would be an easier approach to code (albeit at the expense of some extra CPU time).

Namir
08-29-2023, 03:49 PM
Post: #4
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
Hi, Namir

We could trade speed for space getting Stirling number of the 2nd kind.
Code:
SS(n,k) := sum((-1)^(k-j) * comb(k,j)*j^n, j = 0 .. k); S(n, k) := SS(n,k) / k!;  /* Stiring number 2nd kind */

XCas> [S(5,k) for k in range(1,6)] /* x^5 falling factorial coefficients */

[1, 15, 25, 10, 1]

Bernoulli numbers from Stirling numbers 2nd kind
Code:
B0(m) := {   local k, t:=1;   if (m==1) return -1/2;     if (odd(m)) return 0;     m += 1;   for(k:=2; k<=m; k+=1) t = SS(m,k)/k^2 - t;   return t; }

XCas> map(B0, range(10))

[1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0]

Code:
sum_xm(x,m) := {local k, t := 1/(m+1);   for(k:=m; k>1; k-=1) { t := (x-k)*t + SS(m,k-1)/k!; }   return x*(x-1)*t; }; B2(m,x) := sum_xm(x,m-1) * m + B0(m);

XCas> sum_xm(101, 1);                /* = sum(k, k=0 .. 100) */

5050

XCas> sqrt(sum_xm(101, 3));       /* sum of cubes = square of sum */

5050

XCaS> B2(5, 3.5), B2(3, 5.5);       /* OP example */

220.9375, 123.75
08-29-2023, 05:47 PM (This post was last modified: 08-29-2023 08:24 PM by Albert Chan.)
Post: #5
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-29-2023 11:47 AM)Namir Wrote:  When I saw the nested summation formula in Wikipedia, I decided that it would be an easier approach to code
(albeit at the expense of some extra CPU time).

That's quite true.
But noted that formula had huge cancellation issue with factor (x+k)^m, when x and/or m is big.

For example, if m=6, n=m (for last sum), it may have huge cancellation errors.

$$\displaystyle \frac{\Delta^6(x^6)}{7} = \frac{x^6-6\,(x\!+\!1)^6+15\,(x\!+\!2)^6-20\,(x\!+\!3)^6+15\,(x\!+\!4)^6-6\,(x\!+\!5)^6+(x\!+\!6)^6}{7} = \frac{720}{7}$$

Example, tested on XCas 1.9.0, with 48 bits float:

XCas> B(m,x) := sum(1/(n+1) * sum((-1)^k * comb(n,k) * (x+k)^m, k = 0 .. n), n = 0 .. m)

XCas> float(B(16, pi))      → 1462871.932         /* reference */
XCas> B(16, float(pi))      → -423965.265625    /* massive cancellations */
XCas> B2(16, float(pi))    → 1462871.93177      /* falling factorial form more stable */
08-29-2023, 07:23 PM (This post was last modified: 08-29-2023 07:24 PM by Namir.)
Post: #6
 Namir Senior Member Posts: 1,083 Joined: Dec 2013
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
Albert, you are right. Using the nested summation I wrote in Matlab, B16(pi) gave me 1.543898027587891e+06, a result that was 5.5% higher than the correct answer. Matlab's bernoulli function gave the correct answer. But running the program on two different HP15C emulators on my iPhone gave me -4.5426E11! I am disapointed at the credibility of the nested summation equation.
08-29-2023, 08:21 PM (This post was last modified: 08-29-2023 08:38 PM by Albert Chan.)
Post: #7
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
I just noticed B2(16, float(pi)) cheated with good B0(16) fraction.

Again XCas 1.9.0 default float used 48 bits, biased truncate rounding!
I wished XCas switch back to IEEE double default ...

XCas> float(B0(16))      /* good, exact B0(16) = -3617/510 */

-7.09215686275

XCas> B0(float(16))      /* bad, suffered massive cancellations */

-733244.429688

Perhaps B0(m) need asymptotic series formula ?
Ideas on getting accurate Bernoulli's number (with float, not exact) welcome.

Here's a test that does not use flawed B0(m).

XCas> float(B(15, pi))      → 640433.197324      /* reference */
XCas> B(15, float(pi))      → 594743.757812      /* massive cancellations */
XCas> B2(15, float(pi))    → 640433.197329      /* falling factorial form */

Update:

Last B2() still cheated from perfect Stirling's numbers. Redo with float m:

XCas> B2(15.0, float(pi))  → 640433.053861      /* falling factorial form */
08-29-2023, 09:16 PM
Post: #8
 John Keith Senior Member Posts: 1,033 Joined: Dec 2013
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
The problem here is that all useful methods of calculating Bernoulli numbers involve fairly large numbers, and the classic hp's can only represent numbers < 10^10 exactly. Albert's method using Stirling numbers may be better than nested summation but will still see cancellation errors for n > 16.

One of my favorite methods uses Euler zigzag numbers. The numbers involved are smaller than factorials, and computing the zigzag numbers requires only addition. Also, this method requires only one division at the end, which prevents rounding errors from accumulating. However, i have not compared the two methods directly and there may be little or no improvement in practice.
08-30-2023, 09:58 AM
Post: #9
 Namir Senior Member Posts: 1,083 Joined: Dec 2013
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-29-2023 09:16 PM)John Keith Wrote:  The problem here is that all useful methods of calculating Bernoulli numbers involve fairly large numbers, and the classic hp's can only represent numbers < 10^10 exactly. Albert's method using Stirling numbers may be better than nested summation but will still see cancellation errors for n > 16.

One of my favorite methods uses Euler zigzag numbers. The numbers involved are smaller than factorials, and computing the zigzag numbers requires only addition. Also, this method requires only one division at the end, which prevents rounding errors from accumulating. However, i have not compared the two methods directly and there may be little or no improvement in practice.

Thanks John for the hint and the link. I am off today for a few weeks vacation in Europe (which kicks off with attending the wedding of a nephew in Greece!), so I may be slow in digesting the Euler Sizgzag numbers. They do look interesting!

Namir
08-30-2023, 12:10 PM
Post: #10
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
Falling factorial form Bernoulli functions, version 2
• Stirling number (2nd kind) rounded to integer, because it always is.
• if m is float, it pollute XCas calculations to approximate ... no cheating!
• B0(even m) sum a pair (of alternate signs) at a time, improve accuracy

Code:
S(m,k) := round(sum((-1)^(k-j) * comb(k,j)*j^m, j = 0 .. k) / k!); B0(m) := {   local k, t:=1;   if (m==1) return 1/2 - m;     if (remain(m,2)) return 0;     m += 1;   for(k:=m-m+2; k<=m; k+=2) t += (k*k*S(m,k+1) - (k+1)*S(m,k)) * (k-1)! / (k*k+k);    return t; }; sum_xm(x,m) := {local k, t := 1/(m+1);   for(k:=m-m+m; k>1; k--) t := (x-k)*t + S(m,k-1)/k;   return x*(x-1)*t; }; B2(m,x) := sum_xm(x,m-1) * m + B0(m);

Note: XCas 1.9.0 default was 48 bits float, truncated rounding

XCas> float(B0(16)), B0(float(16))

-7.09215686275, -5.25

XCas> float(B2(16, pi)), B2(float(16), pi)

1462871.93156, 1462873.77393
08-30-2023, 04:30 PM
Post: #11
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-29-2023 01:57 AM)Albert Chan Wrote:  Constant of integration, adjusted for using falling factorial numbers.

B(6) = 1*0!/1 - 63*1!/2 + 301*2!/3 - 350*3!/4 + 140*4!/5 - 21*5!/6 + 1*6!/7 = 1/42

It is not wrong, but we overkill a bit with this formula.
Coefficients are lower factorial form of x^7, not x^6

XCas> f(m) := factor(simplify(x!/(x-m)!));      /* = x^m */
XCas> simplify([1, 63, 301, 350, 140, 21, 1] * map(f, range(1,8)))           → x^7

All we need is lower factorial form of x^6, integrate it (just pretend)

XCas> simplify([1, 31, 90, 65, 15, 1] * map(f, range(1,7)))                       → x^6
XCas> expand([1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * map(f, range(2,8)))

x^7/7 - x^6/2 + x^5/2 - x^3/6 + x/42

XCas> expand(bernoulli(7,x)/7)            /* B^7/7 matched, since B(7) = 0 */

x^7/7 - x^6/2 + x^5/2 - x^3/6 + x/42

Linear term coefficient is B(6) = 1/42, what we wanted.
We differentiate, evaluate it at x = 0, to recover it. (again, pretend)

(x^m)' = (x * (x-1)^(m-1))' = (x-1)^(m-1) + x * ((x-1)^(m-1))'

if x = 0, last term goes away, RHS = (-1)^(m-1) = (-1)^m * m!

This is why we have factorial with alterntating sign factor.

B(6) = [1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * [-1!, 2!, -3!, 4!, -5!, 6!] = 1/42

We have slightly smaller numbers, and 1 less term.
Code:
B0(m) := {   local k, t:=0;   if (m<=1) return 1 - 3/2*m;   if (remain(m,2)) return 0;     for(k:=m-m+2; k<=m; k+=2) t += (k*k*S(m,k) - (k+1)*S(m,k-1)) * (k-1)! / (k*k+k);    return t; };

XCas> map(B0, range(10))

[1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0]

XCas> float(B0(16)), B0(float(16))

-7.09215686275 , -6.5625
08-30-2023, 11:59 PM
Post: #12
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-30-2023 04:30 PM)Albert Chan Wrote:  B(6) = [1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * [-1!, 2!, -3!, 4!, -5!, 6!] = 1/42

We could use horner's rule, and remove factorial function.

B(6) = -(1/2 - 2*(31/3 - 3*(90/4 - 4*(65/5 - 5*(15/6 - 6*(1/7))))))

Code is simpler, calculations faster, and slightly more accurate.

Code:
B0(m) := {   local k, t:=1/(m+1);   if (m<=1) return 1 - 3/2*m;   if (remain(m,2)) return 0;     for(k:=m; k>2; k-=1) t := S(m,k-1)/k - k*t;   return 2*t - 1/2; };

XCas> map(B0, range(10))

[1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0]

XCas> float(B0(16)), B0(float(16))

-7.09215686275, -7.43446095788
08-31-2023, 01:08 PM (This post was last modified: 08-31-2023 01:40 PM by John Keith.)
Post: #13
 John Keith Senior Member Posts: 1,033 Joined: Dec 2013
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-30-2023 11:59 PM)Albert Chan Wrote:
(08-30-2023 04:30 PM)Albert Chan Wrote:  B(6) = [1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * [-1!, 2!, -3!, 4!, -5!, 6!] = 1/42

We could use horner's rule, and remove factorial function.

B(6) = -(1/2 - 2*(31/3 - 3*(90/4 - 4*(65/5 - 5*(15/6 - 6*(1/7))))))

Nice use of Horner's method.

The expression that you are computing, (-1)^k*k!*Stirling2(n+1, k+1), is A163626 which is actually simpler to compute than the Stirling numbers of the second kind. From the linked page, T(n, k) = (k+1)*T(n-1,k) - k*T(n-1,k-1). Thus each term in your expression can be computed with just one subtraction and two multiplications, after which each term must be divided by (k+1) as above.

Here is my RPL implementation, which computes the whole triangle but can be easily modified to do one row at a time. \GDLIST is deltaList and LSEQ returns a list of integers 1..n, equivalent to range(1, n+1).

Code:
 \<< I\->R \-> n   \<< { 1 } n R\->I LSEQ 1. n     FOR k DUP2 1. k SUB * 0 + \GDLIST 1 SWAP + SWAP     NEXT DROP n 1. + \->LIST   \>> \>>

Edit: modified version, returns row n.

Code:
 \<< I\->R \-> n   \<< { 1 } n R\->I LSEQ 1. n               @ Row 0, list of 1..n     FOR k SWAP OVER 1. k SUB * 0 + \GDLIST  @ (k+1)*T(n-1,k) - k*T(n-1,k-1)       1 SWAP + SWAP                         @ Append last term = 1     NEXT DROP                               @ Drop list of 1..n   \>> \>>
08-31-2023, 06:43 PM (This post was last modified: 08-31-2023 06:56 PM by Albert Chan.)
Post: #14
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-31-2023 01:08 PM)John Keith Wrote:
(08-30-2023 11:59 PM)Albert Chan Wrote:  We could use horner's rule, and remove factorial function.

B(6) = -(1/2 - 2*(31/3 - 3*(90/4 - 4*(65/5 - 5*(15/6 - 6*(1/7))))))

The expression that you are computing, (-1)^k*k!*Stirling2(n+1, k+1), is A163626 which is actually simpler to compute than the Stirling numbers of the second kind. From the linked page, T(n, k) = (k+1)*T(n-1,k) - k*T(n-1,k-1).

Let's try B(6) with T(5th row), instead of S(6th row)

XCas> T5 := [1, -31, 180, -390, 360, -120]
XCas> T5 / range(2, len(T5)+2) - T5         → [-1/2, 62/3, -135, 312, -300, 720/7]
XCas> sum(ans())                                    → 1/42

Let's try S (Stirling 2nd kind) 5th row too!

x^5 = [1, 15, 25, 10, 1] • [x^1, x^2, x^3, x^4, x^5]

We have this recurrence relation: S(m,k) = S(m-1,k) * k + S(m-1,k-1)

Code:
B(6) = -(1/2 - 2*(       31/3 - 3*(        90/4 - 4*(        65/5 - 5*(      15/6 - 6*(1/7      = -(1/2 - 2*( (15*2+1)/3 - 3*( (25*3+15)/4 - 4*( (10*4+25)/5 - 5*((1*5+10)/6 - 6*(1/7      = -(1/2 - 2*(15+(1-15)/3 - 3*(25+(15-25)/4 - 4*(10+(25-10)/5 - 5*(1+(10-1)/6 - 6*(1/7      = -(1/2 - 2*(15+(1-15)/3 - 3*(25+(15-25)/4 - 4*(10+(25-10)/5 - 5*(1 + 270/420      = -(1/2 - 2*(15+(1-15)/3 - 3*(25+(15-25)/4 - 4*(5 - 90/420      = -(1/2 - 2*(15+(1-15)/3 - 3*(5 - 690/420)      = -(1/2 - 2*(0 + 110/420      = +1/42

Fractional part scaled by denominator LCM, to make calculations exact, within limits of machine precision.

Code:
B0(m) := {   local k, s, t, h := round(m);   if (m<=1) return 1 - 3/2*m;   if (odd(h)) return 0;     h := lcm(range(2,h+2));   s := [1, 1];          t := [0, h/(m+1)];     /* = [IP, h*FP] */   for(k:=m; k>2; k-=1) {     s[1] := s[0];        /* = S(m-1,k-1) */     s[0] := S(m-1,k-2);       t := [s[1], h/k*(s[0]-s[1])] - k*t;   }   t := 2*t - [0, h/2];   return t[0] + t[1]/h; };

For XCas, both 48 bits and IEEE double version, B0(float(16)) give correctly rounded result.

XCas> float(B0(16)), B0(float(16))

-7.09215686275, -7.09215686275

Just before B(m) sum collapsed to a number, integer part always goes to zero.

x^5 falling factorial form, evaluated at x=-1

[1, 15, 25, 10, 1] * [-1!, 2!, -3!, 4!, -5!] = (-1)^5 = -1

--> [15, 25, 10, 1] * [2!, -3!, 4!, -5!] = 0

With horner's rule, intermediate IP calculations will reach a peak, then fall back down.
Example, for B0(16), there is no overflow. Here's intermediate t values.

t[0] = IP      t[1] = FP * h
[        1,          68108040]
[       90,        2609126520]
[     3290,       52668535920]
[    63700,      609341612880]
[   715078,     4106142920880]
[  4796792,    15498659095920]
[ 19160570,    28183600645200]
[ 44182710,     7883219023680]
[ 55279653,   -44760997521240]
[ 33735702,   -54111311534280]
[  8352708,   -19234067973120]
[   592410,    -1800651604752]
[     5461,      -22288338072]
[        0,         -40384344]

h = lcm(2 .. 17) = 2^4*3^2*5*7*11*13*17 = 12252240
B(16) = 2 * -40384344 / 12252240 - 1/2 = -3617/510
08-31-2023, 07:30 PM
Post: #15
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(08-31-2023 06:43 PM)Albert Chan Wrote:  For XCas, both 48 bits and IEEE double version, B0(float(16)) give correctly rounded result.

Not in XCas with higher precision, even though MPFR guaranteed correctly rounded result.

XCas> Digits := 16
XCas> 21 / float(7) - 3      → ≈ 2.220e-16      ???

XCas> Digits := 100
XCas> 21 / float(7) - 3      → ≈ 2.286e-100     ???

When I run B0(16) with Digits := 16, which has 54-48 = 6 bits more room, I get:

XCas> Digits := 16;
XCas> float(B0(16)), B0(float(16))

-7.092156862745098, -7.101414149575914

Same result for XCas 1.5.0 and 1.9.0. I think this is a serious bug.
09-01-2023, 11:53 AM
Post: #16
 John Keith Senior Member Posts: 1,033 Joined: Dec 2013
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
For B(16) using my program listed above in approximate mode (all floating-point numbers) I get exactly 7.1. The largest term in the 15th row of A163626 is ~10^14 so some rounding is inevitable, and would certainly be worse for 10-digit calculators.

The largest row with all numbers < 10^10 is row 11 which should allow accurate computation of B(12). The value of B(12) to 12 digits is -0.253113553114, while the program returns -0.25312, which has only four correct digits. The exact value of B(12) is -691/2730 so we should expect more accurate results.

Your idea of scaling the terms by their LCM is interesting but I can't see how one could compute the LCM on calculators limited to 10-digit floats.
09-01-2023, 06:18 PM (This post was last modified: 09-01-2023 07:11 PM by Albert Chan.)
Post: #17
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
Hi, John Keith

You are right fixed precision calculator is the not right machine for Bernoulli numbers.
Practically, we should just stored them in an array, not waste time calculating them.

I just considered it a challenge to push it as far as I can.

We started from lower factorial form of x^7, then to x^6:

B(6) = [1,31,90,65,15,1] * [-1!/2, 2!/3, -3!/4, 4!/5, -5!/6, 6!/7]
= -1/2 + 62/3 - 135 + 312 - 300 + 720/7 = 1/42

Here is a trick to reduce size of intermediate numbers, with factorial form of x^5:

B(6) = [1,15,25,10,1] * [1!/(2*3), -2!/(3*4), 3!/(4*5), -4!/(5*6), 5!/(6*7)]
= 1/6 - 5/2 + 15/2 - 8 + 20/7 = 1/42

Again, to push calculations to reach higher B(m), we split the numbers.

Code:
B(6) = 1/6 * (1 - (15 - 3*3/5*(25 - 4*4/6*(10 - 5*5/7*(1      = 1/6 * (1 - (15 - 3*3/5*(25 - 4*4/6*(5 + 600/420      = 1/6 * (1 - (15 - 3*3/5*(5 + 1200/420      = 1/6 * (1 - ( 0 + 360/420      = 1/42

IP(t) = S(m-1,k-1) - k * IP(previous t).
IP from the left,             10 - 5*1 = 5,     25-4*5 = 5,     15-3*5 = 0

15 - 3*(25 - 4*(10 - 5*(1))) = [15,25,10,1] * [2!,-3!,4!,-5!] / 2 = 0

IP will reach a peak, then fall back down (to 0).
We had proven this previously, and will use it into code.
(08-31-2023 06:43 PM)Albert Chan Wrote:  x^5 falling factorial form, evaluated at x=-1

[1, 15, 25, 10, 1] * [-1!, 2!, -3!, 4!, -5!] = (-1)^5 = -1

--> [15, 25, 10, 1] * [2!, -3!, 4!, -5!] = 0

Code:
B0(m) := {   local k, t, h := round(m);   if (m<=1) return 1 - 3/2*m;   if (odd(h)) return 0;     h := lcm(range(4, h+2));   t := [1, 0];                  /* = [IP, h*FP] */   for(k:=m-1; k>2; k-=1)      t := [S(m-1,k-1),-(t[0]*h+t[1])/(k+2)*(k*k)] - t[0]*k*[1,-h];       return (h - t[1]) / (6*h);    /* t[0] = 0, guaranteed */ };

XCas> float(B0(16)), B0(float(16))

-7.09215686275, -7.09215686275

Normally we pre-compute LCM (to machine limits), and simply hard coded into program.
I use LCM functions here only because I wanted to try XCas MPFR.
09-01-2023, 06:58 PM (This post was last modified: 09-01-2023 09:43 PM by Albert Chan.)
Post: #18
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(09-01-2023 06:18 PM)Albert Chan Wrote:  Here is a trick to reduce size of intermediate numbers, with factorial form of x^5:

B(6) = [1,15,25,10,1] * [1!/(2*3), -2!/(3*4), 3!/(4*5), -4!/(5*6), 5!/(6*7)]
= 1/6 - 5/2 + 15/2 - 8 + 20/7 = 1/42

We can keep going! From falling factorial of x^4:

B(6) = [1,7,6,1] * [0*1!/(2*3*4), -1*2!/(3*4*5), 2*3!/(4*5*6), -3*4!/(5*6*7)]
= 0 - 7/30 + 3/5 - 12/35 = 1/42

Or, with T numbers, 3rd row, remove RHS sign and factorial

B(6) = [1,-7,12,-6] * [0*1/(2*3*4), 1*2/(3*4*5), 2*3/(4*5*6), 3*4/(5*6*7)]
= 0 - 7/30 + 3/5 - 12/35 = 1/42

Simply lower Striling 2nd kind by 2 rows work well.

Code:
B0(m) := {   local k, t := m*0+1;   if (m<=1) return 1 - 3/2*m;   if (odd(round(m))) return 0;     for(k:=m-2; k>2; k-=1) t := S(m-2,k-1) - (k-1)*k*k/((k-2)*(k+3))*t;   return t / when(m<3, 6, -30); }

XCas 1.9.0 48 bits float, round by truncation

XCas> for (m=2; m<=20; m+=2) print(m, float(B0(m)), B0(float(m)));

2,     0.166666666667,    0.166666666667
4,    -0.0333333333333,  -0.0333333333333
6,     0.0238095238095,   0.0238095238095
8,    -0.0333333333333,  -0.0333333333333
10,    0.0757575757576,   0.0757575757513
12,   -0.253113553114,   -0.253113554362
14,    1.16666666667,     1.1666665062
16,   -7.09215686275,    -7.09217076724
18,   54.9711779449,     54.971246923
20, -529.124242424,    -523.164965208
09-02-2023, 05:44 PM (This post was last modified: 09-03-2023 01:51 PM by Albert Chan.)
Post: #19
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(09-01-2023 11:53 AM)John Keith Wrote:  The largest row with all numbers < 10^10 is row 11 which should allow accurate computation of B(12). The value of B(12) to 12 digits is -0.253113553114, while the program returns -0.25312, which has only four correct digits. The exact value of B(12) is -691/2730 so we should expect more accurate results.

We may be able to correct for tiny errors. B(m) largest denominator is 2*(2^m-1).

d = 2*(2^12-1) = 8190 (largest d)
d * -0.25312 = -2073.0528            → B(12) = -2073/8190 = -691/2730

Denominator likely smaller, but required work. d = product(p, for (p-1) | even m)
This is why d is divisible by 6. 2-1=1, 3-1=2, both divides even m.

d = 2*3*5*7*13 = 2730
d * -0.25312 = -691.0176             → B(12) = -691/2730

Amazingly, we get correct B(m = 12), without correction (algorithm based on above link)
If m divisible by 4, B(m) is negative. (positive otherwise)

lua> m = 12
lua> d = 2*(2^m-1)
lua> b = 4*d/(d+2) * fac(m)/pi^m
lua> 1+b, d
2073.4899880820685      8190
lua> _ + b*2^-m
2073.995967083065

--> B(12) = -2073/8190 = -691/2730

lua> m = 16
lua> d = 2*(2^m-1)
lua> b = 4*d/(d+2) * fac(m)/pi^m
lua> 1+b, d
929555.7943024995      131070
lua> _ + b*2^-m
929569.9781830278
lua> _ + b*3^-m
929569.9997771184

--> B(16) = -929569/131070 = -3617/510

Rounding errors might cause off-by-1 error with B(m) numerator.
It may be safer to drop numbers by 1/2, equivalent to rounding without 1+

Code:
gcd = require'factor'.gcd fac = fn'x: tgamma(x+1)' roughB = fn'm: I.real(-2*fac(m)/(2*pi*I)^m) * (1+2^-m+3^-m)' ratioB = fn'm,n,d: d=2*(2^m-1); n=round(roughB(m)*d); m=gcd(n,d); n/m, d/m'

lua> for m=2, 26, 2 do print(m, ratioB(m)) end

2           1     6
4          -1    30
6           1    42
8          -1    30
10          5    66
12       -691  2730
14          7     6
16      -3617   510
18      43867   798
20    -174611   330
22     854513   138
24 -236364091  2730
26    8553103     6

IEEE double (53 bits mantissa) is able to get upto B(26) exact fractions.
For B(28), n reached 54 bits, always even, but numerator should be odd.

roughB(big m) still good, just not enough precision to get back fractions.

lua> roughB(32)      → -15116315767.092175
lua> roughB(36)      → -13711655205088.352
09-03-2023, 01:20 AM (This post was last modified: 09-11-2023 09:05 PM by Albert Chan.)
Post: #20
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
The reason we sidetrack for Bernoulli numbers is we wanted to use this:

With previous post for bernoulli numbers (O(1), speed and storage), we have:
Code:
function b0(m)     if m%2==1 then return m==1 and -1/2 or 0 end     if m == 0 then return 1 end     local b = (-1)^(m/2+1) * 2*tgamma(m+1)/pi^m     local d = 2^m-1     local n = floor(b + b*3^-m)     return (n+0.5)/d, b/(d+1) -- ratio = zeta(m) end
Code:
function b(m, x)     local c, t = 1, 1     for k = 1, m do         c = c*(m-k+1)/k         t = t*x + c*b0(k)     end     return t end

lua> b(3, 5.5)
123.75
lua> b(5, 3.5)
220.9375

lua> b(15, pi)
640433.1973240786
lua> b(16, pi)
1462871.9320025162

Update 1: halved d, and use floor instead of round

b0(m) = (n+0.5)/d = (2n+1)/(2d) = odd/even, as expected

Update 2: zeta = (1+3^-m)/(1-2^-m) estimate is better than (1+2^-m+3^-m)
We could implement this without cost (actually, cheaper!)
Also, I remove complex number math for sign of b0(m)

Update 3: return bad b0(m) too, before zeta correction. Ratio is zeta(m)

lua> good, bad = b0(2)
lua> good, bad, good/bad -- = zeta(2)
0.16666666666666666      0.10132118364233778      1.6449340668482262
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)