Post Reply 
(HP15C)(HP67)(HP41C) Bernoulli Polynomials
09-01-2023, 06:58 PM (This post was last modified: 09-01-2023 09:43 PM by Albert Chan.)
Post: #18
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
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-01-2023 06:58 PM



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