Post Reply 
(HP15C)(HP67)(HP41C) Bernoulli Polynomials
09-01-2023, 06:18 PM (This post was last modified: 09-01-2023 07:11 PM by Albert Chan.)
Post: #17
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
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.
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:18 PM



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