(HP15C)(HP67)(HP41C) Bernoulli Polynomials - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (HP15C)(HP67)(HP41C) Bernoulli Polynomials (/thread-20416.html) Pages: 1 2 |
(HP15C)(HP67)(HP41C) Bernoulli Polynomials - Namir - 08-28-2023 01:10 PM Program calculates Bernoulli polynomials using nested summations (see Wikipedia). The memory map is: Code:
The listing for HP15C is: Code:
The listing for HP67 is: Code:
Listing for HP-41C is: Code:
Example 1 5 ENTER 3.5 f A output is 220.9375 Example 2 3 ENTER 5.5 f A output is 123.7500 Note: I was thinking about writing a second version of the listing for the HP-14C that uses the ISG command for registers 2 and 3. The problem is whenever the listing uses RCL 2 and RCL 3 it must be followed by the INT command. This adds more steps IMHO. So I abandoned the ISG-using version. RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-29-2023 01:57 AM Code: B(m,x) := sum(1/(n+1) * sum((-1)^k * comb(n,k) * (x+k)^m, k = 0 .. n), n = 0 .. m) Another way to get Bm(x) is by synthetic division, for falling factorial form. We can think of B as "derivative" of sum of power function, symbol Bm(x) == Bm Σ(k^m, k=0 .. x-1) = ∫(Bm) = Bm+1 / (m+1) + C (08-05-2019 12:21 AM)Albert Chan Wrote: From problem 12, Σ(xm, x=0 to n-1) = nm+1 / (m+1) Except for symbol names, the 2 sums look very similar! We can easily get x^m to its falling factorial form, via synthetic division. The only issue is to get the integration constant, Bm(0) = B(m) B(1) = -1/2, B(2k+1) = 0, we only need to worry about B(2k) (07-28-2019 12:02 AM)Albert Chan Wrote: Sum of k^6 formula = \(1\binom{n}{1}+63\binom{n}{2}+602\binom{n}{3}+2100\binom{n}{4}+3360\binom{n}{5}+2520\binom{n}{6}+720\binom{n}{7}\) Example, lets do formula for B6(x) Synthetic Division, x^m to falling factorial form, for Bm formula Code:
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 x^5 = x^4 * x, we pick numbers 3rd entries from the right. x^5 = [1, 15, 25, 10, 1] • [x^1, x^2, x^3, x^4, x^5] (B^6-B(6))/6 = [1/2, 15/3, 25/4, 10/5, 1/6] • [x^2, x^3, x^4, x^5, x^6] // "integrate" B^6 = 1/42 + [3, 30, 75/2, 12, 1] • [x^2, x^3, x^4, x^5, x^6] Or, in efficient "horner" form: B6(x) = 1/42 + x*(x-1)*(3 + (x-2)*(30 + (x-3)*(75/2 + (x-4)*(12 + (x-5))))) RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Namir - 08-29-2023 11:47 AM A few weeks ago I published here a version of the HP-17B that evaluates Bernoulli polynomials. The BASIC code had a subroutine that calculates Bernoulli numbers and stores then in an array. The program allowed the user to repeatedly enter values for the polynomial order and for x. It remembers that last order used. If that is case, then it simply uses the current values in the Bernoulii numbers array. I started thinking about using the same approach but realized the storage scheme is a bit quirky (and limited the order of the Bernoulli polynomials) since odd orders (except for 1) of the Bernoulli numbers are zero. When I saw the nested summation formula in Wikipedia, I decided that it would be an easier approach to code (albeit at the expense of some extra CPU time). Namir RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-29-2023 03:49 PM Hi, Namir We could trade speed for space getting Stirling number of the 2nd kind. Code: SS(n,k) := sum((-1)^(k-j) * comb(k,j)*j^n, j = 0 .. k); XCas> [S(5,k) for k in range(1,6)] /* x^5 falling factorial coefficients */ [1, 15, 25, 10, 1] Bernoulli numbers from Stirling numbers 2nd kind Code: B0(m) := { XCas> map(B0, range(10)) [1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0] Code: sum_xm(x,m) := {local k, t := 1/(m+1); XCas> sum_xm(101, 1); /* = sum(k, k=0 .. 100) */ 5050 XCas> sqrt(sum_xm(101, 3)); /* sum of cubes = square of sum */ 5050 XCaS> B2(5, 3.5), B2(3, 5.5); /* OP example */ 220.9375, 123.75 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-29-2023 05:47 PM (08-29-2023 11:47 AM)Namir Wrote: When I saw the nested summation formula in Wikipedia, I decided that it would be an easier approach to code That's quite true. But noted that formula had huge cancellation issue with factor (x+k)^m, when x and/or m is big. For example, if m=6, n=m (for last sum), it may have huge cancellation errors. \(\displaystyle \frac{\Delta^6(x^6)}{7} = \frac{x^6-6\,(x\!+\!1)^6+15\,(x\!+\!2)^6-20\,(x\!+\!3)^6+15\,(x\!+\!4)^6-6\,(x\!+\!5)^6+(x\!+\!6)^6}{7} = \frac{720}{7} \) Example, tested on XCas 1.9.0, with 48 bits float: XCas> B(m,x) := sum(1/(n+1) * sum((-1)^k * comb(n,k) * (x+k)^m, k = 0 .. n), n = 0 .. m) XCas> float(B(16, pi)) → 1462871.932 /* reference */ XCas> B(16, float(pi)) → -423965.265625 /* massive cancellations */ XCas> B2(16, float(pi)) → 1462871.93177 /* falling factorial form more stable */ RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Namir - 08-29-2023 07:23 PM Albert, you are right. Using the nested summation I wrote in Matlab, B16(pi) gave me 1.543898027587891e+06, a result that was 5.5% higher than the correct answer. Matlab's bernoulli function gave the correct answer. But running the program on two different HP15C emulators on my iPhone gave me -4.5426E11! I am disapointed at the credibility of the nested summation equation. RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-29-2023 08:21 PM I just noticed B2(16, float(pi)) cheated with good B0(16) fraction. Again XCas 1.9.0 default float used 48 bits, biased truncate rounding! I wished XCas switch back to IEEE double default ... XCas> float(B0(16)) /* good, exact B0(16) = -3617/510 */ -7.09215686275 XCas> B0(float(16)) /* bad, suffered massive cancellations */ -733244.429688 Perhaps B0(m) need asymptotic series formula ? Ideas on getting accurate Bernoulli's number (with float, not exact) welcome. Here's a test that does not use flawed B0(m). XCas> float(B(15, pi)) → 640433.197324 /* reference */ XCas> B(15, float(pi)) → 594743.757812 /* massive cancellations */ XCas> B2(15, float(pi)) → 640433.197329 /* falling factorial form */ Update: Last B2() still cheated from perfect Stirling's numbers. Redo with float m: XCas> B2(15.0, float(pi)) → 640433.053861 /* falling factorial form */ RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 08-29-2023 09:16 PM The problem here is that all useful methods of calculating Bernoulli numbers involve fairly large numbers, and the classic hp's can only represent numbers < 10^10 exactly. Albert's method using Stirling numbers may be better than nested summation but will still see cancellation errors for n > 16. One of my favorite methods uses Euler zigzag numbers. The numbers involved are smaller than factorials, and computing the zigzag numbers requires only addition. Also, this method requires only one division at the end, which prevents rounding errors from accumulating. However, i have not compared the two methods directly and there may be little or no improvement in practice. RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Namir - 08-30-2023 09:58 AM (08-29-2023 09:16 PM)John Keith Wrote: The problem here is that all useful methods of calculating Bernoulli numbers involve fairly large numbers, and the classic hp's can only represent numbers < 10^10 exactly. Albert's method using Stirling numbers may be better than nested summation but will still see cancellation errors for n > 16. Thanks John for the hint and the link. I am off today for a few weeks vacation in Europe (which kicks off with attending the wedding of a nephew in Greece!), so I may be slow in digesting the Euler Sizgzag numbers. They do look interesting! Namir RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-30-2023 12:10 PM Falling factorial form Bernoulli functions, version 2
Code: S(m,k) := round(sum((-1)^(k-j) * comb(k,j)*j^m, j = 0 .. k) / k!); 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 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-30-2023 04:30 PM (08-29-2023 01:57 AM)Albert Chan Wrote: Constant of integration, adjusted for using falling factorial numbers. 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) := { 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 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-30-2023 11:59 PM (08-30-2023 04:30 PM)Albert Chan Wrote: B(6) = [1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * [-1!, 2!, -3!, 4!, -5!, 6!] = 1/42 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)))))) Code is simpler, calculations faster, and slightly more accurate. Code: B0(m) := { 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, -7.43446095788 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 08-31-2023 01:08 PM (08-30-2023 11:59 PM)Albert Chan Wrote:(08-30-2023 04:30 PM)Albert Chan Wrote: B(6) = [1/2, 31/3, 90/4, 65/5, 15/6, 1/7] * [-1!, 2!, -3!, 4!, -5!, 6!] = 1/42 Nice use of Horner's method. 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). Thus each term in your expression can be computed with just one subtraction and two multiplications, after which each term must be divided by (k+1) as above. Here is my RPL implementation, which computes the whole triangle but can be easily modified to do one row at a time. \GDLIST is deltaList and LSEQ returns a list of integers 1..n, equivalent to range(1, n+1). Code:
Edit: modified version, returns row n. Code:
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-31-2023 06:43 PM (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. 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 Fractional part scaled by denominator LCM, to make calculations exact, within limits of machine precision. Code: B0(m) := { 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 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 08-31-2023 07:30 PM (08-31-2023 06:43 PM)Albert Chan Wrote: For XCas, both 48 bits and IEEE double version, B0(float(16)) give correctly rounded result. Not in XCas with higher precision, even though MPFR guaranteed correctly rounded result. XCas> Digits := 16 XCas> 21 / float(7) - 3 → ≈ 2.220e-16 ??? XCas> Digits := 100 XCas> 21 / float(7) - 3 → ≈ 2.286e-100 ??? When I run B0(16) with Digits := 16, which has 54-48 = 6 bits more room, I get: XCas> Digits := 16; XCas> float(B0(16)), B0(float(16)) -7.092156862745098, -7.101414149575914 Same result for XCas 1.5.0 and 1.9.0. I think this is a serious bug. RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 09-01-2023 11:53 AM For B(16) using my program listed above in approximate mode (all floating-point numbers) I get exactly 7.1. The largest term in the 15th row of A163626 is ~10^14 so some rounding is inevitable, and would certainly be worse for 10-digit calculators. The largest row with all numbers < 10^10 is row 11 which should allow accurate computation of B(12). The value of B(12) to 12 digits is -0.253113553114, while the program returns -0.25312, which has only four correct digits. The exact value of B(12) is -691/2730 so we should expect more accurate results. Your idea of scaling the terms by their LCM is interesting but I can't see how one could compute the LCM on calculators limited to 10-digit floats. RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-01-2023 06:18 PM 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 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 Code: B0(m) := { 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. RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-01-2023 06:58 PM (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: 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) := { 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 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-02-2023 05:44 PM (09-01-2023 11:53 AM)John Keith Wrote: The largest row with all numbers < 10^10 is row 11 which should allow accurate computation of B(12). The value of B(12) to 12 digits is -0.253113553114, while the program returns -0.25312, which has only four correct digits. The exact value of B(12) is -691/2730 so we should expect more accurate results. We may be able to correct for tiny errors. B(m) largest denominator is 2*(2^m-1). d = 2*(2^12-1) = 8190 (largest d) d * -0.25312 = -2073.0528 → B(12) = -2073/8190 = -691/2730 Denominator likely smaller, but required work. d = product(p, for (p-1) | even m) This is why d is divisible by 6. 2-1=1, 3-1=2, both divides even m. d = 2*3*5*7*13 = 2730 d * -0.25312 = -691.0176 → B(12) = -691/2730 Amazingly, we get correct B(m = 12), without correction (algorithm based on above link) If m divisible by 4, B(m) is negative. (positive otherwise) lua> m = 12 lua> d = 2*(2^m-1) lua> b = 4*d/(d+2) * fac(m)/pi^m lua> 1+b, d 2073.4899880820685 8190 lua> _ + b*2^-m 2073.995967083065 --> B(12) = -2073/8190 = -691/2730 lua> m = 16 lua> d = 2*(2^m-1) lua> b = 4*d/(d+2) * fac(m)/pi^m lua> 1+b, d 929555.7943024995 131070 lua> _ + b*2^-m 929569.9781830278 lua> _ + b*3^-m 929569.9997771184 --> B(16) = -929569/131070 = -3617/510 Rounding errors might cause off-by-1 error with B(m) numerator. It may be safer to drop numbers by 1/2, equivalent to rounding without 1+ Code: gcd = require'factor'.gcd lua> for m=2, 26, 2 do print(m, ratioB(m)) end 2 1 6 4 -1 30 6 1 42 8 -1 30 10 5 66 12 -691 2730 14 7 6 16 -3617 510 18 43867 798 20 -174611 330 22 854513 138 24 -236364091 2730 26 8553103 6 IEEE double (53 bits mantissa) is able to get upto B(26) exact fractions. For B(28), n reached 54 bits, always even, but numerator should be odd. roughB(big m) still good, just not enough precision to get back fractions. lua> roughB(32) → -15116315767.092175 lua> roughB(36) → -13711655205088.352 RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-03-2023 01:20 AM The reason we sidetrack for Bernoulli numbers is we wanted to use this: With previous post for bernoulli numbers (O(1), speed and storage), we have: Code: function b0(m) Code: function b(m, x) lua> b(3, 5.5) 123.75 lua> b(5, 3.5) 220.9375 lua> b(15, pi) 640433.1973240786 lua> b(16, pi) 1462871.9320025162 Update 1: halved d, and use floor instead of round b0(m) = (n+0.5)/d = (2n+1)/(2d) = odd/even, as expected Update 2: zeta = (1+3^-m)/(1-2^-m) estimate is better than (1+2^-m+3^-m) We could implement this without cost (actually, cheaper!) Also, I remove complex number math for sign of b0(m) Update 3: return bad b0(m) too, before zeta correction. Ratio is zeta(m) lua> good, bad = b0(2) lua> good, bad, good/bad -- = zeta(2) 0.16666666666666666 0.10132118364233778 1.6449340668482262 |