(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[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