(12C+) Bernoulli Number - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (12C+) Bernoulli Number (/thread-13348.html) (12C+) Bernoulli Number - Gamo - 07-27-2019 06:41 AM In need of the Bernoulli Number using HP-12C ? Here is an attempt to generate a Bernoulli Number constant using 12C Without a Pi function this program use 355/113 which give out about 4 to 5 digits precision. ------------------------------------------------- To run: If you need to know B10 divide it by 2 is 5 5 [R/S] display 0.07576 [R/S] 5 [X<>Y] 66 Answer: B10 is 0.07576 or in fraction is 5/66 -------------------------------------------------- B12 12 ÷ 2 = 6 6 [R/S] display 0.25311 [R/S] 61 [X<>Y] 241 Answer: B12 since 12 is divisible by 4 answer is Negative -0.25311 in fraction is -61/241 -------------------------------------------------- Remark: To find B(n) divide it by 2 and calculate. This program do not give answer of the alternate negative value such as B2 = 1/6 where B4 = -1/30 For B(n) that divisible by 4 answer is "Negative" -------------------------------------------------- Program: Code:  2 x STO 2 355 ENTER 113 ÷ STO 3 1 STO 0 STO 1 ----------------- RCL 0  // line 16 2 x 1 + RCL 2 CHS Y^X RCL 1 + RCL 1 X<>Y X≤Y GTO 34 STO 1 1 STO+0 GTO 16 RCL 2  // Line 34 n! 2 x RCL 1 x 2 RCL 2 Y^X 1 - RCL 2 RCL 3 X<>Y Y^X x ÷  // B(n) constant end here ---------------- R/S  // Line 51 STO 0  // Decimal to Fraction start here STO 1 0 STO 2 1 RCL 0 INTG  // Line 58 RCL 2 x + STO 2 RCL 0 x . 5 // decimal and five // Line 66 + INTG STO 3 RCL 2 ÷ RND RCL 0 RND -  // Subtract sign X=0 GTO 86 Rv  // Roll Down RCL 2 X<>Y RCL 1 FRAC 1/x STO 1 GTO 58 RCL 2  // Line 86 RCL 3 GTO 00  // Line 88 Formula use to calculate Bernoulli Number B(n) = [2(2n)! ÷ ((2^2n) - 1)(Pi^2n)] [1 + (1/3^2n) + (1/5^2n) + ...] Gamo RE: (12C+) Bernoulli Number - Albert Chan - 07-27-2019 12:41 PM (07-27-2019 06:41 AM)Gamo Wrote:  Formula use to calculate Bernoulli Number B(n) = [2(2n)! ÷ ((2^2n) - 1)(Pi^2n)] [1 + (1/3^2n) + (1/5^2n) + ...] Hi, Gamo I think there is a typo: B(n) should be abs(B(2n)) Can you provide source ? From https://wstein.org/edu/fall05/168/projects/kevin_mcgown/bernproj.pdf zeta(x) = sum(k^-x, where k = 1 to ∞) = 1 / product(1 - p^-x, where p is prime) B(2n) = (-1)^(n+1) * 2*(2n)!/(2 pi)^(2n) * zeta(2n) RE: (12C+) Bernoulli Number - Gamo - 07-27-2019 01:40 PM Hello Albert Chan Thanks for the review. Here is the formula: [attachment=7559] Look like this formula is only good for Small B(n) Gamo RE: (12C+) Bernoulli Number - John Keith - 07-27-2019 07:49 PM It looks like an exact formula to me but it is an infinite series, and I have no idea how fast it converges. The Wikipedia page on Bernoulli numbers has lots of useful information including some exact algorithms. The fastest methods do seem to require a lot of memory so they may not be suitable for a small calculator like the 12C. RE: (12C+) Bernoulli Number - Albert Chan - 07-28-2019 12:02 AM I lookup "A Source Book in Mathematics", chapter "On the Bernoulli numbers": B(n) is just sum of k^n formula linear term coefficient. Example, this is how B(6) is calculated, by doing k^6 forward difference (see thread: https://www.hpmuseum.org/forum/thread-12248-post-110972.html#pid110972) 1 64 729 4096 15625 46656 117649 // value of 1^6 to 7^6 63 665 3367 11529 31031 70993 ﻿ ﻿ ﻿ ﻿ // forward differences 602 2702 8162 19502 39962 2100 5460 11340 20460 3360 5880 9120 2520 3240 720 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}$$ B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42 Correction: B(n) is sum of k^n formula, Σ(i^k, i = 0 to n-1) linear term coefficient B(1) = linear coefficient of n(n+1)/2 - n = 1/2 - 1 = -1/2 (see http://www.mikeraugh.org/Talks/BernoulliSummation-LACC.pdf, slide 31 to 34) The Python code (next post), work for n=1 because it flipped sign for all odd n. The "bug" actually fix the sign RE: (12C+) Bernoulli Number - Albert Chan - 07-28-2019 01:08 AM Python code for Bernoulli Number (also return the forward difference table) Code: from fractions import Fraction def B_fdiff(n):     lst = [1,1] + [i**n for i in range(2, n+2)]     for i in range(n):         for j in range(n,i,-1): lst[j+1] -= lst[j]  # forward diff         lst[0] = Fraction(lst[i+2],i+2) - lst[0]    # B(n) value     return lst >>> B_fdiff(6) [Fraction(1, 42), 1, 63, 602, 2100, 3360, 2520, 720] >>> for i in range(13): print i, B_fdiff(i)[0] ... 0 1 1 -1/2 2 1/6 3 0 4 -1/30 5 0 6 1/42 7 0 8 -1/30 9 0 10 5/66 11 0 12 -691/2730 RE: (12C+) Bernoulli Number - Gamo - 07-28-2019 02:29 AM I use this same program on RPN-67 app the HP-67 emulator. This RPN-67 got the printer feature so I did the loop print out. The result give good approximate answer from B2 to B112 B114 just give out incorrect answer to a whole number 107407655 the correct answer should be 5.17567436175..... Correct answer continue from B116 to B124 Then B126 and over all result are 1.000000000 *** Is this have to do with the limited value of the Factorial function? Clip: https://youtu.be/LH7yDHan08M Gamo RE: (12C+) Bernoulli Number - John Keith - 07-28-2019 11:21 AM (07-28-2019 12:02 AM)Albert Chan Wrote:  I lookup "A Source Book in Mathematics", chapter "On the Bernoulli numbers": B(n) is just sum of k^n formula linear term coefficient. Example, this is how B(6) is calculated, by doing k^6 forward difference (see thread: https://www.hpmuseum.org/forum/thread-12248-post-110972.html#pid110972) 1 64 729 4096 15625 46656 117649 // value of 1^6 to 7^6 63 665 3367 11529 31031 70993 ﻿ ﻿ ﻿ ﻿ // forward differences 602 2702 8162 19502 39962 2100 5460 11340 20460 3360 5880 9120 2520 3240 720 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}$$ B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42 That is a very neat method, I was not aware of that one. However, it seems that all similar exact methods require n+(n-1) storage registers to calculate B(n) since one needs to keep the (n-1)th row of the difference table in memory while calculating the nth row. EDIT: I tried your program as well as the Akiyama-Tanigawa method as used in the third program here on the HP-48G and both methods fail due to catastrophic rounding error for n>10. These methods may only be practical for languages that use exact rational arithmetic. RE: (12C+) Bernoulli Number - Albert Chan - 07-31-2019 05:14 PM Knowing the pattern of sk(n) = (Σi^k, k=0 to n-1) = n^(k+1)/(k+1) - n^k/2 + k/12 * n^(k-1) + 0 * x^(k-2) + ..., we can get Bernoulli constants another way. Example: to get upto B(6), s6(n) = n^7/7 - n^6/2 + n^5/2 + a*n^3 + b*n, where a, b are unknown → a*n^2 + b = s6(n)/n - n^4*(n^2/7 - (n-1)/2) For n=1, eqn1 = a + b = 0/1 - 1*(1/7 - 0/2) = -1/7 For n=2, eqn2 = 4a + b = 1/2 - 16*(4/7 - 1/2) = -9/14 eqn2 - eqn1: 3a = -9/14 + 2/14 = -1/2 → a = -1/6 B(6) = b = -1/7 - a = 1/42 B(4) = a * 3 / $$\binom{6}{4}$$ = -1/30 B(2) = (1/2) * 5 / $$\binom{6}{2}$$ = 1/6 B(1) = (-1/2) * 6 / $$\binom{6}{1}$$ = -1/2 B(0) = (1/7) * 7 / $$\binom{6}{0}$$ = 1 RE: (12C+) Bernoulli Number - Albert Chan - 08-30-2023 09:46 PM (07-28-2019 11:21 AM)John Keith Wrote:  EDIT: I tried your program as well as the Akiyama-Tanigawa method as used in the third program here on the HP-48G and both methods fail due to catastrophic rounding error for n>10. If we use scaled numbers, we can push it a bit more. Example, this is Akiyama-Tanigawa method, scaled by denominator LCM. Note: XCas 1.9.0 default float now is 48 bits, truncated rounding. Note: round(L[0]) does not change L[0] value, just its type, real → integer. XCas> m := 22; XCas> h := lcm(range(2, m+2))         → 5354228880 XCas> L := float(h ./ range(2, m+2)) XCas> for(k:=2; k<=m; k++) {                 L := -deltalist(L) .* range(1,len(L));                 if (even(k)) print(k, round(L[0])/h);           }; 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 48 bits float mangaged to get upto B(22) IEEE double can reach B(24), all exact fractions. RE: (12C+) Bernoulli Number - Albert Chan - 09-11-2023 03:48 PM (07-28-2019 12:02 AM)Albert Chan Wrote:  1 64 729 4096 15625 46656 117649 // value of 1^6 to 7^6 63 665 3367 11529 31031 70993 ﻿ ﻿ ﻿ ﻿ // forward differences 602 2702 8162 19502 39962 2100 5460 11340 20460 3360 5880 9120 2520 3240 720 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}$$ B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42 Note: above Σk^6 formula is from 0 to n = B(7,n)/7 + n^6 We may use symmetry, and do 0^6 to 3^6, then flip it, to get values (-3)^6 .. (+3)^6 This keep table numbers smaller, thus smaller forward difference numbers: $$x^6 = 729\binom{x+3}{0} -665\binom{x+3}{1} +602\binom{x+3}{2} -540\binom{x+3}{3} +480\binom{x+3}{4} -360\binom{x+3}{5} +720\binom{x+3}{6}$$ Integrate (actually sum, but look like integration) $$\large \frac{B^7}{7} \normalsize = 729\binom{x+3}{1} -665\binom{x+3}{2} +602\binom{x+3}{3} -540\binom{x+3}{4} +480\binom{x+3}{5} -360\binom{x+3}{6} +720\binom{x+3}{7}$$ + C Differentiate, evaluated at x=0, we get B(6) = RHS constant term. First half of terms is not as easy to evaluate, because there was no factor x (with a factor x, derivative at x=0 make the term goes away!) But, it may still worth it, to keep forward difference table entries small. Note: for B(m), n = floor(m/2), k = 1 .. n $$\displaystyle \binom{x+n}{k} = \left(\frac{x+n}{k}\right) \binom{x+n-1}{k-1}$$ $$\displaystyle \binom{x+n}{k}' \bigg\rvert_{x=0} = \left(\frac{1}{k}\right) \binom{n-1}{k-1} + \left(\frac{n}{k}\right) \binom{x+n-1}{k-1}'\bigg\rvert_{x=0}$$ I had managed to undo recursion, into iterative sums (written prove welcome!) $$\displaystyle \binom{x+n}{k}' \bigg\rvert_{x=0} \;=\; \sum_{j=1}^{k} \frac{(-1)^{j+1}}{j}\; \binom{n}{k-j} \;=\; \binom{n}{k}\;\sum_{j=n-(k-1)}^{n} \frac{1}{j}$$ Note that coefficient formula good for half terms (the other half is easy, see code) Code: from gmpy2 import mpq, mpz, xmpz from itertools import izip, count def bernoulli(m):     if m == 0: return 1     lst = [0, 1] + [ mpz(i)**m for i in range(2, (m+3)//2) ]     lst = [(-x if m&1 else x) for x in lst[-1-(m&1):0:-1]] + lst     lst = map(xmpz, lst)        # speedup fdiff_inplace     fdiff_inplace(lst)          # forward diff table     return sum(c*d for (c,d) in izip(bernoulli_coef(m//2), lst)) def fdiff_inplace(f):     n = len(f) - 1     for i in xrange(n):         # n forward diff(s)         for j in xrange(n,i,-1):             f[j] -= f[j-1] def bernoulli_coef(n, mpq=mpq):     'diff(comb(x+n,k)), k=1,2,3,..., evaluate at x=0'     c, s, N = mpz(1), 0, n+1     for j in xrange(n, 0, -1):  # tricky half         c = c * j // (N-j)         s += mpq(1, j)         yield c * s     s = mpq(1, N)     for c in count(1):          # easy half         yield s         s *= mpq(-c, N+c) I did not optimize away B(odd m), just to see if it matches ... it does! B(1) = -1/2, B(odd>1) = 0 >>> for m in range(2,31,2): print m, bernoulli(m) ... 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 28     -23749461029L/870 30     8615841276005L/14322 RE: (12C+) Bernoulli Number - Albert Chan - 09-12-2023 05:59 PM Bernoulli number related links ... Compare with forward difference table, Stirling numbers 2nd kind numbers are smaller see thread, (HP15C)(HP67)(HP41C) Bernoulli Polynomials Euler zigzag numbers (A) smaller still. B(m) = A(m-1) * (-1)^(m/2-1) * m / (2^m*(2^m-1)) B(6) = 16 * (-1)^2 * 6 / (2^6*(2^6-1)) = 1/42 see thread, (28 48 49 50) Bernoulli Numbers BCMATH programs for some number theory functions CalcBn V3.0 from The Bernoulli Number Page This use zeta function correction, same as this thread, but *much* faster! On my laptop, it get B(100,000) exact fraction in 12 second!