(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[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