(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[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