(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[Namir]

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.

[Albert Chan]

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 B_m(x) is by synthetic division, for falling factorial form.

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

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

(08-05-2019 12:21 AM)Albert Chan Wrote:  From problem 12, Σ(x^m, x=0 to n-1) = n^m+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, B_m(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 B₆(x)

Synthetic Division, x^m to falling factorial form, for B_m 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:

B₆(x) = 1/42 + x*(x-1)*(3 + (x-2)*(30 + (x-3)*(75/2 + (x-4)*(12 + (x-5)))))

[Namir]

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

[Albert Chan]

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

[Albert Chan]

(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 */

[Namir]

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.

[Albert Chan]

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

[John Keith]

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.

[Namir]

(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

[Albert Chan]

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

[Albert Chan]

(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

[Albert Chan]

(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

[John Keith]

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

[Albert Chan]

(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

[Albert Chan]

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

[John Keith]

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.

[Albert Chan]

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.

[Albert Chan]

(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

[Albert Chan]

(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

[Albert Chan]

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