Bernoulli numbers are used in some series expansions of several functions (trigonometric, hyperbolic, gamma, etc.), and are extremely important in number theory and analysis.
Note that there are two definitions of Bernoulli numbers; this task will be using the modern usage (as per The National Institute of Standards and Technology convention).
The nth Bernoulli number is expressed as Bn.
The Akiyama–Tanigawa algorithm for the "second Bernoulli numbers" as taken fromwikipedia is as follows:
formfrom 0by 1tondoA[m] ← 1/(m+1)forjfrommby -1to 1doA[j-1] ←j×(A[j-1] -A[j])returnA[0] (which isBn)
Using a GMP thick binding available athttp://www.codeforge.com/article/422541
WITHGMP.Rationals,GMP.Integers,Ada.Text_IO,Ada.Strings.Fixed,Ada.Strings;USEGMP.Rationals,GMP.Integers,Ada.Text_IO,Ada.Strings.Fixed,Ada.Strings;PROCEDUREMainISFUNCTIONBernoulli_Number(N:Natural)RETURNUnbounded_FractionISFUNCTION"/"(Left,Right:Natural)RETURNUnbounded_FractionIS(To_Unbounded_Integer(Left)/To_Unbounded_Integer(Right));A:ARRAY(0..N)OFUnbounded_Fraction;BEGINFORMIN0..NLOOPA(M):=1/(M+1);FORJINREVERSE1..MLOOPA(J-1):=(J/1)*(A(J-1)-A(J));ENDLOOP;ENDLOOP;RETURNA(0);ENDBernoulli_Number;BEGINFORIIN0..60LOOPIFIMOD2=0ORI=1THENDECLAREB:Unbounded_Fraction:=Bernoulli_Number(I);S:String:=Image(GMP.Rationals.Numerator(B));BEGINPut_Line("B ("&(IFI<10THEN" "ELSE"")&Trim(I'Img,Left)&")="&(44-S'Length)*" "&Image(B));END;ENDIF;ENDLOOP;ENDMain;
B(0)= 1 / 1B(1)= 1 / 2B(2)= 1 / 6B(4)= -1 / 30B(6)= 1 / 42B(8)= -1 / 30B(10)= 5 / 66B(12)= -691 / 2730B(14)= 7 / 6B(16)= -3617 / 510B(18)= 43867 / 798B(20)= -174611 / 330B(22)= 854513 / 138B(24)= -236364091 / 2730B(26)= 8553103 / 6B(28)= -23749461029 / 870B(30)= 8615841276005 / 14322B(32)= -7709321041217 / 510B(34)= 2577687858367 / 6B(36)= -26315271553053477373 / 1919190B(38)= 2929993913841559 / 6B(40)= -261082718496449122051 / 13530B(42)= 1520097643918070802691 / 1806B(44)= -27833269579301024235023 / 690B(46)= 596451111593912163277961 / 282B(48)= -5609403368997817686249127547 / 46410B(50)= 495057205241079648212477525 / 66B(52)= -801165718135489957347924991853 / 1590B(54)= 29149963634884862421418123812691 / 798B(56)= -2479392929313226753685415739663229 / 870B(58)= 84483613348880041862046775994036021 / 354B(60)=-1215233140483755572040304994079820246041491 / 56786730
Uses the LONG LONG INT mode of Algol 68G which allows large precision integers.
BEGIN # Show Bernoulli numbers B0 to B60 as rational numbers # # Uses code from the Arithmetic/Rational task modified to use # # LONG LONG INT to allow for the large number of digits requried # PR precision 100 PR # sets the precision of LONG LONG INT # # Code from the Arithmetic/Rational task # # ============================================================== # MODE FRAC = STRUCT( LONG LONG INT num #erator#, den #ominator#); PROC gcd = (LONG LONG INT a, b) LONG LONG INT: # greatest common divisor # (a = 0 | b |: b = 0 | a |: ABS a > ABS b | gcd(b, a MOD b) | gcd(a, b MOD a)); PROC lcm = (LONG LONG INT a, b)LONG LONG INT: # least common multiple # a OVER gcd(a, b) * b; PRIO // = 9; # higher then the ** operator # OP // = (LONG LONG INT num, den)FRAC: ( # initialise and normalise # LONG LONG INT common = gcd(num, den); IF den < 0 THEN ( -num OVER common, -den OVER common) ELSE ( num OVER common, den OVER common) FI ); OP + = (FRAC a, b)FRAC: ( LONG LONG INT common = lcm(den OF a, den OF b); FRAC result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common ); num OF result//den OF result ); OP - = (FRAC a, b)FRAC: a + -b; OP - = (FRAC frac)FRAC: (-num OF frac, den OF frac); # ============================================================== # # end code from the Arithmetic/Rational task # # Additional FRACrelated operators # OP * = ( INT a, FRAC b )FRAC: ( num OF b * a ) // den OF b; OP // = ( INT a, INT b )FRAC: LONG LONG INT( a ) // LONG LONG INT( b ); # returns the nth Bernoulli number, n must be >= 0 # # Uses the algorithm suggested by the task, so B(1) is +1/2 # PROC bernoulli = ( INT n )FRAC: IF n < 0 THEN # n is out of range # 0 // 1 ELSE # n is valid # [ 0 : n ]FRAC a; FOR i FROM LWB a TO UPB a DO a[ i ] := 0 // 1 OD; FOR m FROM 0 TO n DO a[ m ] := 1 // ( m + 1 ); FOR j FROM m BY -1 TO 1 DO a[ j - 1 ] := j * ( a[ j - 1 ] - a[ j ] ) OD OD; a[ 0 ] FI # bernoulli # ; FOR n FROM 0 TO 60 DO FRAC bn := bernoulli( n ); IF num OF bn /= 0 THEN # have a non-0 Bn # print( ( "B(", whole( n, -2 ), ") ", whole( num OF bn, -50 ) ) ); print( ( " / ", whole( den OF bn, 0 ), newline ) ) FI ODENDB( 0) 1 / 1B( 1) 1 / 2B( 2) 1 / 6B( 4) -1 / 30B( 6) 1 / 42B( 8) -1 / 30B(10) 5 / 66B(12) -691 / 2730B(14) 7 / 6B(16) -3617 / 510B(18) 43867 / 798B(20) -174611 / 330B(22) 854513 / 138B(24) -236364091 / 2730B(26) 8553103 / 6B(28) -23749461029 / 870B(30) 8615841276005 / 14322B(32) -7709321041217 / 510B(34) 2577687858367 / 6B(36) -26315271553053477373 / 1919190B(38) 2929993913841559 / 6B(40) -261082718496449122051 / 13530B(42) 1520097643918070802691 / 1806B(44) -27833269579301024235023 / 690B(46) 596451111593912163277961 / 282B(48) -5609403368997817686249127547 / 46410B(50) 495057205241079648212477525 / 66B(52) -801165718135489957347924991853 / 1590B(54) 29149963634884862421418123812691 / 798B(56) -2479392929313226753685415739663229 / 870B(58) 84483613348880041862046775994036021 / 354B(60) -1215233140483755572040304994079820246041491 / 56786730
To be able to handle numbers up to B(60) and beyond, this script represents the numbers with lists whose items are the digit values — which of course requires custom math routines.
onbernoullis(n)-- Return a list of "numerator / denominator" texts representing Bernoulli numbers B(0) to B(n).setlistMathScripttogetListMathScript(10)-- Script object providing custom list math routines.setoutputto{}-- Akiyama–Tanigawa algorithm for the "second Bernoulli numbers".-- List 'a' will contain {numerator, denominator} lists representing fractions.-- The numerators and denominators will in turn be lists containing integers representing their (decimal) digits.setato{}repeatwithmfrom0ton-- Append the structure for 1 / (m + 1) to the end of a.set{numerator2,denominator2}to{{1},listMathScript'sintToList(m+1)}seta'sendtoresultrepeatwithjfrommto1by-1-- Retrieve the preceding numerator and denominator.set{numerator1,denominator1}toa'sitemjtelllistMathScript-- Get the two fractions' lowest common denominator and adjust the numerators accordingly.setlcdtoitslcm(denominator1,denominator2)setnumerator1toitsmultiply(numerator1,its|div|(lcd,denominator1))setnumerator2toitsmultiply(numerator2,its|div|(lcd,denominator2))-- Subtract numerator2 from numerator1 and multiply the result by j.-- Assign the results to numerator2 and denominator2 for the next iteration.setnumerator2toitsmultiply(itssubtract(numerator1,numerator2),itsintToList(j))setdenominator2tolcdendtell-- Also store them in a's slot j. No need to reduce them here.seta'sitemjto{numerator2,denominator2}endrepeat-- The fraction just stored in a's first slot is Bernoulli(m). Reduce it and append a text representation to the output.telllistMathScriptsetgcdtoitshcf(numerator2,denominator2)setnumerator2toits|div|(numerator2,gcd)setdenominator2toits|div|(denominator2,gcd)setendofoutputtoitslistToText(numerator2)&(" / "&itslistToText(denominator2))endtellendrepeatreturnoutputendbernoullisongetListMathScript(base)scriptonmultiply(lst1,lst2)-- Multiply lst1 by lst2.setlst1Lengthto(countlst1)setlst2Lengthto(countlst2)setproductLengthtolst1Length+lst2Length-1setproductto{}repeatproductLengthtimessetproduct'sendto0endrepeat-- Long multiplication algorithm, updating product digits on the fly instead of summing rows at the end.repeatwithlst2Indexfrom-1to-lst2Lengthby-1setlst2Digittolst2'sitemlst2Indexif(lst2Digitis not0)thensetcarryto0setproductIndextolst2Indexrepeatwithlst1Indexfromlst1'slengthto1by-1telllst2Digit*(lst1'sitemlst1Index)+carry+(product'sitemproductIndex)setproduct'sitemproductIndexto(itmodbase)setcarryto(itdivbase)endtellsetproductIndextoproductIndex-1endrepeatif(carry=0)thenelseif(productIndex<-productLength)thensetproduct'sbeginningtocarryelsesetproduct'sitemproductIndexto(product'sitemproductIndex)+carryendifendifendrepeatreturnproductendmultiplyonsubtract(lst1,lst2)-- Subtract lst2 from lst1.setlst1Lengthto(countlst1)setlst2Lengthto(countlst2)-- Pad copies to equal lengths.copylst1tolst1repeat(lst2Length-lst1Length)timessetlst1'sbeginningto0endrepeatcopylst2tolst2repeat(lst1Length-lst2Length)timessetlst2'sbeginningto0endrepeat-- Is lst2's numeric value greater than lst1's?setpaddedLengthto(countlst1)repeatwithifrom1topaddedLengthsetlst1Digittolst1'sitemisetlst2Digittolst2'sitemisetlst2Greaterto(lst2Digit>lst1Digit)if((lst2Greater)or(lst1Digit>lst2Digit))thenexitrepeatendrepeat-- If so, set up to subtract lst1 from lst2 instead. We'll invert the result's sign at the end.if(lst2Greater)thentelllst2setlst2tolst1setlst1toitendtell-- The subtraction at last!setdifferenceto{}setborrowto0repeatwithifrompaddedLengthto1by-1tell(lst1'sitemi)+base-borrow-(lst2'sitemi)setdifference'sbeginningto(itmodbase)setborrowto1-(itdivbase)endtellendrepeatif(lst2Greater)theninvert(difference)returndifferenceendsubtracton|div|(lst1,lst2)-- List lst1 div lst2.returndivide(lst1,lst2)'squotientend|div|on|mod|(lst1,lst2)-- List lst1 mod lst2.returndivide(lst1,lst2)'sremainderend|mod|ondivide(lst1,lst2)-- Divide lst1 by lst2. Return a record containing separate lists for the quotient and remainder.setdividendtotrim(lst1)setdivisortotrim(lst2)setdividendLengthto(countdividend)setdivisorLengthto(countdivisor)if(divisorLength>dividendLength)thenreturn{quotient:{0},remainder:dividend}-- Note the dividend's and divisor's signs, but use absolute values in the division.setdividendNegativeto(dividend'sbeginning<0)if(dividendNegative)theninvert(dividend)setdivisorNegativeto(divisor'sbeginning<0)if(divisorNegative)theninvert(divisor)-- Long-division algorithm, but quotient digits are subtraction counts.setquotientto{}if(divisorLength>1)thensetremaindertodividend'sitems1thru(divisorLength-1)elsesetremainderto{}endifrepeatwithnextSlotfromdivisorLengthtodividendLengthsetremainder'sendtodividend'sitemnextSlotrepeatwithsubtractionCountfrom0tobase-- Only ever reaches base - 1.setsubtractionResulttotrim(subtract(remainder,divisor))if(subtractionResult'sbeginning<0)thenexitrepeatsetremaindertosubtractionResultendrepeatsetendofquotienttosubtractionCountendrepeat-- The quotient's negative if the input signs are different. Positive otherwise.if(dividendNegative≠divisorNegative)theninvert(quotient)-- The remainder has the same sign as the dividend.if(dividendNegative)theninvert(remainder)return{quotient:quotient,remainder:remainder}enddivideonlcm(lst1,lst2)-- Lowest common multiple of lst1 and lst2.returnmultiply(lst2,|div|(lst1,hcf(lst1,lst2)))endlcmonhcf(lst1,lst2)-- Highest common factor of lst1 and lst2.setlst1totrim(lst1)setlst2totrim(lst2)repeatuntil(lst2={0})setxtolst1setlst1tolst2setlst2totrim(|mod|(x,lst2))endrepeatif(lst1'sbeginning<0)theninvert(lst1)returnlst1endhcfoninvert(lst)-- Invert the sign of all lst's "digits".repeatwiththisDigitinlstsetthisDigit'scontentsto-thisDigitendrepeatendinvertontrim(lst)-- Return a copy of lst with no leading zeros.repeatwithifrom1to(countlst)if(lst'sitemiis not0)thenexitrepeatendrepeatreturnlst'sitemsithruendendtrimonintToList(n)-- Return a list of numbers representing n's digits.setlstto{nmodbase}setntondivbaserepeatuntil(n=0)setbeginningoflsttonmodbaseasintegersetntondivbaseendrepeatreturnlstendintToListonlistToText(lst)-- Return the number represented by the input list as text.-- This lazily assumes 2 <= base <= 10. :)setlsttotrim(lst)if(lst'sbeginning<0)theninvert(lst)setlst'sbeginningto"-"endifreturnjoin(lst,"")endlistToTextendscriptreturnresultendgetListMathScriptonjoin(lst,delim)setastidtoAppleScript'stext item delimiterssetAppleScript'stext item delimiterstodelimsettxttolstastextsetAppleScript'stext item delimiterstoastidreturntxtendjoinontask()setmaxNto60setoutputto{""}setpaddingto" = "setbernoulliNumberstobernoullis(maxN)repeatwithnfrom0tomaxNsetbernietobernoulliNumbers'sitem(n+1)if(berniedoesnotstart with"0")thensetBnto"B("&n&")"setoutput'sendtoBn&¬text((countBn)-3)thru(50-(offsetof"/"inbernie))ofpadding&¬bernieendifendrepeatreturnjoin(output,linefeed)endtasktask()
"B(0) = 1 / 1B(1) = 1 / 2B(2) = 1 / 6B(4) = -1 / 30B(6) = 1 / 42B(8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730"
( BernoulliList = B Bs answer indLn indexLen indexPadding , n numberPadding p solPos solidusPos sp . ( B = m A a j b . -1:?m & :?A & whl ' ( 1+!m:~>!arg:?m & ((!m+1:?j)^-1:?a) map $ ( ( = .(-1+!j:?j)*(!arg+-1*!a):?a ) . !A ) : ?A ) & !A:? @?b & !b ) & -1:?n & :?Bs & whl ' ( 1+!n:~>!arg:?n & B$!n !Bs:?Bs ) & @(!arg:? [?indexLen) & 1+!indexLen:?indexLen & !Bs:%@(?:? "/" [?solidusPos ?) ? & 1+!solidusPos:?solidusPos:?p & :?sp & whl ' (!p+-1:~<0:?p&" " !sp:?sp) & :?answer & whl ' ( !Bs:%?B ?Bs & ( !B:0 | (!B:/|str$(!B "/1"):?B) & @(!B:? "/" [?solPos ?) & @(!arg:? [?indLn) & !sp : ? [(-1*!indexLen+!indLn) ?indexPadding : ? [(-1*!solidusPos+!solPos) ?numberPadding & "B(" !arg ")=" !indexPadding !numberPadding (!B:>0&" "|) !B \n !answer : ?answer ) & -1+!arg:?arg ) & str$!answer )& BernoulliList$60;B(0)= 1/1B(1)= 1/2B(2)= 1/6B(4)= -1/30B(6)= 1/42B(8)= -1/30B(10)= 5/66B(12)= -691/2730B(14)= 7/6B(16)= -3617/510B(18)= 43867/798B(20)= -174611/330B(22)= 854513/138B(24)= -236364091/2730B(26)= 8553103/6B(28)= -23749461029/870B(30)= 8615841276005/14322B(32)= -7709321041217/510B(34)= 2577687858367/6B(36)= -26315271553053477373/1919190B(38)= 2929993913841559/6B(40)= -261082718496449122051/13530B(42)= 1520097643918070802691/1806B(44)= -27833269579301024235023/690B(46)= 596451111593912163277961/282B(48)= -5609403368997817686249127547/46410B(50)= 495057205241079648212477525/66B(52)= -801165718135489957347924991853/1590B(54)= 29149963634884862421418123812691/798B(56)= -2479392929313226753685415739663229/870B(58)= 84483613348880041862046775994036021/354B(60)=-1215233140483755572040304994079820246041491/56786730
#include<stdlib.h>#include<gmp.h>#define mpq_for(buf, op, n)\ do {\ size_t i;\ for (i = 0; i < (n); ++i)\ mpq_##op(buf[i]);\ } while (0)voidbernoulli(mpq_trop,unsignedintn){unsignedintm,j;mpq_t*a=malloc(sizeof(mpq_t)*(n+1));mpq_for(a,init,n+1);for(m=0;m<=n;++m){mpq_set_ui(a[m],1,m+1);for(j=m;j>0;--j){mpq_sub(a[j-1],a[j],a[j-1]);mpq_set_ui(rop,j,1);mpq_mul(a[j-1],a[j-1],rop);}}mpq_set(rop,a[0]);mpq_for(a,clear,n+1);free(a);}intmain(void){mpq_trop;mpz_tn,d;mpq_init(rop);mpz_inits(n,d,NULL);unsignedinti;for(i=0;i<=60;++i){bernoulli(rop,i);if(mpq_cmp_ui(rop,0,1)){mpq_get_num(n,rop);mpq_get_den(d,rop);gmp_printf("B(%-2u) = %44Zd / %Zd\n",i,n,d);}}mpz_clears(n,d,NULL);mpq_clear(rop);return0;}
B(0 ) = 1 / 1B(1 ) = -1 / 2B(2 ) = 1 / 6B(4 ) = -1 / 30B(6 ) = 1 / 42B(8 ) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
Translation of the C implementation
usingMpir.NET;usingSystem;namespaceBernoulli{classProgram{privatestaticvoidbernoulli(mpq_trop,uintn){mpq_t[]a=newmpq_t[n+1];for(uinti=0;i<n+1;i++){a[i]=newmpq_t();}for(uintm=0;m<=n;++m){mpir.mpq_set_ui(a[m],1,m+1);for(uintj=m;j>0;--j){mpir.mpq_sub(a[j-1],a[j],a[j-1]);mpir.mpq_set_ui(rop,j,1);mpir.mpq_mul(a[j-1],a[j-1],rop);}mpir.mpq_set(rop,a[0]);}}staticvoidMain(string[]args){mpq_trop=newmpq_t();mpz_tn=newmpz_t();mpz_td=newmpz_t();for(uinti=0;i<=60;++i){bernoulli(rop,i);if(mpir.mpq_cmp_ui(rop,0,1)!=0){mpir.mpq_get_num(n,rop);mpir.mpq_get_den(d,rop);Console.WriteLine(string.Format("B({0, 2}) = {1, 44} / {2}",i,n,d));}}Console.ReadKey();}}}
B(0 ) = 1 / 1B(1 ) = -1 / 2B(2 ) = 1 / 6B(4 ) = -1 / 30B(6 ) = 1 / 42B(8 ) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
usingSystem;usingSystem.Console;usingSystem.Linq;usingMathNet.Numerics;namespaceRosettacode.Rational.CS{classProgram{privatestaticreadonlyFunc<int,BigRational>ℚ=BigRational.FromInt;privatestaticBigRationalCalculateBernoulli(intn){vara=InitializeArray(n);foreach(varminEnumerable.Range(1,n)){a[m]=ℚ(1)/(ℚ(m)+ℚ(1));for(varj=m;j>=1;j--){a[j-1]=ℚ(j)*(a[j-1]-a[j]);}}returna[0];}privatestaticBigRational[]InitializeArray(intn){vara=newBigRational[n+1];for(varx=0;x<a.Length;x++){a[x]=ℚ(x+1);}returna;}staticvoidMain(){Enumerable.Range(0,61)// the second parameter is the number of range elements, and is not the final item of the range..Select(n=>new{N=n,BernoulliNumber=CalculateBernoulli(n)}).Where(b=>!b.BernoulliNumber.Numerator.IsZero).Select(b=>string.Format("B({0, 2}) = {1, 44} / {2}",b.N,b.BernoulliNumber.Numerator,b.BernoulliNumber.Denominator)).ToList().ForEach(WriteLine);}}}
B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
Algo based on the example provided in the header of this RC page (the one from Wikipedia).
Extra feature - one can override the default of 60 by supplying a suitable number on the command line. The column widths are not hard-coded, but will adapt to the widths of the items listed.
usingSystem;usingSystem.Numerics;usingSystem.Collections.Generic;namespacebern{classProgram{structBerNum{publicintindex;publicBigIntegerNumer,Denomin;};staticintw1=1,w2=1;// widths for formatting outputstaticintmax=60;// default maximum, can override on command line// returns nth Bernoulli numberstaticBerNumCalcBernoulli(intn){BerNumres;BigIntegerf;BigInteger[]nu=newBigInteger[n+1],de=newBigInteger[n+1];for(intm=0;m<=n;m++){nu[m]=1;de[m]=m+1;for(intj=m;j>0;j--)if((f=BigInteger.GreatestCommonDivisor(nu[j-1]=j*(de[j]*nu[j-1]-de[j-1]*nu[j]),de[j-1]*=de[j]))!=BigInteger.One){nu[j-1]/=f;de[j-1]/=f;}}res.index=n;res.Numer=nu[0];res.Denomin=de[0];w1=Math.Max(n.ToString().Length,w1);// ratchet up widthsw2=Math.Max(res.Numer.ToString().Length,w2);if(max>50)Console.Write(".");// progress dots appear for larger valuesreturnres;}staticvoidMain(string[]args){List<BerNum>BNumbList=newList<BerNum>();// defaults to 60 when no (or invalid) command line parameter is presentif(args.Length>0){int.TryParse(args[0],outmax);if(max<1||max>Int16.MaxValue)max=60;if(args[0]=="0")max=0;}for(inti=0;i<=max;i++)// fill list with values{BerNumBNumb=CalcBernoulli(i);if(BNumb.Numer!=BigInteger.Zero)BNumbList.Add(BNumb);}if(max>50)Console.WriteLine();stringstrFmt="B({0, "+w1.ToString()+"}) = {1, "+w2.ToString()+"} / {2}";// display formatted listforeach(BerNumbninBNumbList)Console.WriteLine(strFmt,bn.index,bn.Numer,bn.Denomin);if(System.Diagnostics.Debugger.IsAttached)Console.Read();}}}
Default (nothing entered on command line):
.............................................................B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
Output with "8" entered on command line:
B(0) = 1 / 1B(1) = 1 / 2B(2) = 1 / 6B(4) = -1 / 30B(6) = 1 / 42B(8) = -1 / 30
Output with "126" entered on the command line:
...............................................................................................................................B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B( 10) = 5 / 66B( 12) = -691 / 2730B( 14) = 7 / 6B( 16) = -3617 / 510B( 18) = 43867 / 798B( 20) = -174611 / 330B( 22) = 854513 / 138B( 24) = -236364091 / 2730B( 26) = 8553103 / 6B( 28) = -23749461029 / 870B( 30) = 8615841276005 / 14322B( 32) = -7709321041217 / 510B( 34) = 2577687858367 / 6B( 36) = -26315271553053477373 / 1919190B( 38) = 2929993913841559 / 6B( 40) = -261082718496449122051 / 13530B( 42) = 1520097643918070802691 / 1806B( 44) = -27833269579301024235023 / 690B( 46) = 596451111593912163277961 / 282B( 48) = -5609403368997817686249127547 / 46410B( 50) = 495057205241079648212477525 / 66B( 52) = -801165718135489957347924991853 / 1590B( 54) = 29149963634884862421418123812691 / 798B( 56) = -2479392929313226753685415739663229 / 870B( 58) = 84483613348880041862046775994036021 / 354B( 60) = -1215233140483755572040304994079820246041491 / 56786730B( 62) = 12300585434086858541953039857403386151 / 6B( 64) = -106783830147866529886385444979142647942017 / 510B( 66) = 1472600022126335654051619428551932342241899101 / 64722B( 68) = -78773130858718728141909149208474606244347001 / 30B( 70) = 1505381347333367003803076567377857208511438160235 / 4686B( 72) = -5827954961669944110438277244641067365282488301844260429 / 140100870B( 74) = 34152417289221168014330073731472635186688307783087 / 6B( 76) = -24655088825935372707687196040585199904365267828865801 / 30B( 78) = 414846365575400828295179035549542073492199375372400483487 / 3318B( 80) = -4603784299479457646935574969019046849794257872751288919656867 / 230010B( 82) = 1677014149185145836823154509786269900207736027570253414881613 / 498B( 84) = -2024576195935290360231131160111731009989917391198090877281083932477 / 3404310B( 86) = 660714619417678653573847847426261496277830686653388931761996983 / 6B( 88) = -1311426488674017507995511424019311843345750275572028644296919890574047 / 61410B( 90) = 1179057279021082799884123351249215083775254949669647116231545215727922535 / 272118B( 92) = -1295585948207537527989427828538576749659341483719435143023316326829946247 / 1410B( 94) = 1220813806579744469607301679413201203958508415202696621436215105284649447 / 6B( 96) = -211600449597266513097597728109824233673043954389060234150638733420050668349987259 / 4501770B( 98) = 67908260672905495624051117546403605607342195728504487509073961249992947058239 / 6B(100) = -94598037819122125295227433069493721872702841533066936133385696204311395415197247711 / 33330B(102) = 3204019410860907078243020782116241775491817197152717450679002501086861530836678158791 / 4326B(104) = -319533631363830011287103352796174274671189606078272738327103470162849568365549721224053 / 1590B(106) = 36373903172617414408151820151593427169231298640581690038930816378281879873386202346572901 / 642B(108) = -3469342247847828789552088659323852541399766785760491146870005891371501266319724897592306597338057 / 209191710B(110) = 7645992940484742892248134246724347500528752413412307906683593870759797606269585779977930217515 / 1518B(112) = -2650879602155099713352597214685162014443151499192509896451788427680966756514875515366781203552600109 / 1671270B(114) = 21737832319369163333310761086652991475721156679090831360806110114933605484234593650904188618562649 / 42B(116) = -309553916571842976912513458033841416869004128064329844245504045721008957524571968271388199595754752259 / 1770B(118) = 366963119969713111534947151585585006684606361080699204301059440676414485045806461889371776354517095799 / 6B(120) = -51507486535079109061843996857849983274095170353262675213092869167199297474922985358811329367077682677803282070131 / 2328255930B(122) = 49633666079262581912532637475990757438722790311060139770309311793150683214100431329033113678098037968564431 / 6B(124) = -95876775334247128750774903107542444620578830013297336819553512729358593354435944413631943610268472689094609001 / 30B(126) = 5556330281949274850616324408918951380525567307126747246796782304333594286400508981287241419934529638692081513802696639 / 4357878
/** * Configured with: --prefix=/Library/Developer/CommandLineTools/usr --with-gxx-include-dir=/usr/include/c++/4.2.1 * Apple LLVM version 9.1.0 (clang-902.0.39.1) * Target: x86_64-apple-darwin17.5.0 * Thread model: posix*/#include<boost/multiprecision/cpp_int.hpp> // 1024bit precision#include<boost/rational.hpp> // Rationals#include<iostream> // formatting with std::cout#include<vector> // Containertypedefboost::rational<boost::multiprecision::int1024_t>rational;// reduce boilerplaterationalbernoulli(size_tn){autoout=std::vector<rational>();for(size_tm=0;m<=n;m++){out.emplace_back(1,(m+1));// automatically constructs objectfor(size_tj=m;j>=1;j--){out[j-1]=rational(j)*(out[j-1]-out[j]);}}returnout[0];}intmain(){for(size_tn=0;n<=60;n+=n>=2?2:1){autob=bernoulli(n);std::cout<<"B("<<std::right<<std::setw(2)<<n<<") = ";std::cout<<std::right<<std::setw(44)<<b.numerator();std::cout<<" / "<<b.denominator()<<std::endl;}return0;}
B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
nstest-project-intellij.core(:gen-class))(defna-t[n]" Used Akiyama-Tanigawa algorithm with a single loop rather than double nested loop "" Clojure does fractional arithmetic automatically so that part is easy "(loop[m0jmA(vec(map#(/1%)(range1(+n2))))]; Prefil A(m) with 1/(m+1), for m = 1 to n(cond; Three way conditional allows single loop(>=j1)(recurm(decj)(assocA(decj)(*j(-(nthA(decj))(nthAj))))); A[j-1] ← j×(A[j-1] - A[j]) ;(<mn)(recur(incm)(incm)A); increment m, reset j = m:else(nthA0))))(defnformat-ans[ans]" Formats answer so that '/' is aligned for all answers "(if(=ans1)(format"%50d / %8d"11)(format"%50d / %8d"(numeratorans)(denominatorans))));; Generate a set of results for [0 1 2 4 ... 60](doseq[q(flatten[01(range2622)]):let[ans(a-tq)]](printlnq":"(format-ansans)))
0 : 1 / 11 : 1 / 22 : 1 / 64 : -1 / 306 : 1 / 428 : -1 / 3010 : 5 / 6612 : -691 / 273014 : 7 / 616 : -3617 / 51018 : 43867 / 79820 : -174611 / 33022 : 854513 / 13824 : -236364091 / 273026 : 8553103 / 628 : -23749461029 / 87030 : 8615841276005 / 1432232 : -7709321041217 / 51034 : 2577687858367 / 636 : -26315271553053477373 / 191919038 : 2929993913841559 / 640 : -261082718496449122051 / 1353042 : 1520097643918070802691 / 180644 : -27833269579301024235023 / 69046 : 596451111593912163277961 / 28248 : -5609403368997817686249127547 / 4641050 : 495057205241079648212477525 / 6652 : -801165718135489957347924991853 / 159054 : 29149963634884862421418123812691 / 79856 : -2479392929313226753685415739663229 / 87058 : 84483613348880041862046775994036021 / 35460 : -1215233140483755572040304994079820246041491 / 56786730
An implementation of the simple algorithm.
Be advised that the pseudocode algorithm specifies (j * (a[j-1] - a[j])) in the inner loop; implementing that as-is gives the wrong value (1/2) where n = 1, whereas subtracting a[j]-a[j-1] yields the correct value (B[1]=-1/2). Seethe numerator list.
(defunbernouilli(n)(loopwitha=(make-array(list(1+n)))formfrom0tondo(setf(arefam)(/1(+m1)))(loopforjfrommdownto1do(setf(arefa(-j1))(*j(-(arefaj)(arefa(-j1))))))finally(return(arefa0))));;Print outputs to stdout:(loopfornfrom0to60do(let((b(bernouillin)))(when(not(zeropb))(formatt"~a: ~a~%"nb))));;For the "extra credit" challenge, we need to align the slashes.(let(results);;collect the results(loopfornfrom0to60do(let((b(bernouillin)))(when(not(zeropb))(push(consbn)results))));;parse the numerators into strings; save the greatest length in max-length(let((max-length(apply#'max(mapcar(lambda(r)(length(formatnil"~a"(numeratorr))))(mapcar#'carresults)))));;Print the numbers with using the fixed-width formatter: ~Nd, where N is;;the number of leading spaces. We can't just pass in the width variable;;but we can splice together a formatting string that includes it.;;We also can't use the fixed-width formatter on a ratio, so we have to split;;the ratio and splice it back together like idiots.(loopfornin(mapcar#'cdr(reverseresults))forrin(mapcar#'car(reverseresults))do(formatt(concatenate'string"B(~2d): ~"(formatnil"~a"max-length)"d/~a~%")n(numeratorr)(denominatorr)))))
B( 0): 1/1B( 1): -1/2B( 2): 1/6B( 4): -1/30B( 6): 1/42B( 8): -1/30B(10): 5/66B(12): -691/2730B(14): 7/6B(16): -3617/510B(18): 43867/798B(20): -174611/330B(22): 854513/138B(24): -236364091/2730B(26): 8553103/6B(28): -23749461029/870B(30): 8615841276005/14322B(32): -7709321041217/510B(34): 2577687858367/6B(36): -26315271553053477373/1919190B(38): 2929993913841559/6B(40): -261082718496449122051/13530B(42): 1520097643918070802691/1806B(44): -27833269579301024235023/690B(46): 596451111593912163277961/282B(48): -5609403368997817686249127547/46410B(50): 495057205241079648212477525/66B(52): -801165718135489957347924991853/1590B(54): 29149963634884862421418123812691/798B(56): -2479392929313226753685415739663229/870B(58): 84483613348880041862046775994036021/354B(60): -1215233140483755572040304994079820246041491/56786730
require"big"classBernoulliincludeIterator(Tuple(Int32,BigRational))definitialize@a=[]ofBigRational@m=0enddefnext@a<<BigRational.new(1,@m+1)@m.downto(1){|j|@a[j-1]=j*(@a[j-1]-@a[j])}v=@m.odd?&&@m!=1?BigRational.new(0,1):@a.firstreturn{@m,v}ensure@m+=1endendb=Bernoulli.newbn=b.first(61).to_amax_width=bn.map{|_,v|v.numerator.to_s.size}.maxbn.reject{|i,v|v.zero?}.eachdo|i,v|puts"B(%2i) = %*i/%i"%[i,max_width,v.numerator,v.denominator]end
Version 1: compute each number separately.
require"big"defbernoulli(n)ar=[]ofBigRational(0..n).eachdo|m|ar<<BigRational.new(1,m+1)m.downto(1){|j|ar[j-1]=j*(ar[j-1]-ar[j])}endar[0]# (which is Bn)endb_nums=(0..61).map{|i|bernoulli(i)}width=b_nums.map{|b|b.numerator.to_s.size}.maxb_nums.each_with_index{|b,i|puts"B(%2i) = %*i/%i"%[i,width,b.numerator,b.denominator]unlessb.zero?}
Version 2: create faster generator to compute array of numbers once.
require"big"defbernoulli2(limit)ar=[]ofBigRational(0..limit).eachdo|m|ar<<BigRational.new(1,m+1)m.downto(1){|j|ar[j-1]=j*(ar[j-1]-ar[j])}yieldar[0]# use Bn value in required blockendendb_nums=[]ofBigRationalbernoulli2(61){|b|b_nums<<b}width=b_nums.map{|b|b.numerator.to_s.size}.maxb_nums.each_with_index{|b,i|puts"B(%2i) = %*i/%i"%[i,width,b.numerator,b.denominator]unlessb.zero?}
B( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
This uses the D module from the Arithmetic/Rational task.
importstd.stdio,std.range,std.algorithm,std.conv,arithmetic_rational;autobernoulli(inuintn)purenothrow/*@safe*/{autoA=newRational[n+1];foreach(immutablem;0..n+1){A[m]=Rational(1,m+1);foreach_reverse(immutablej;1..m+1)A[j-1]=j*(A[j-1]-A[j]);}returnA[0];}voidmain(){immutableberns=61.iota.map!bernoulli.enumerate.filter!(t=>t[1]).array;immutablewidth=berns.map!(b=>b[1].numerator.text.length).reduce!max;foreach(immutableb;berns)writefln("B(%2d) = %*d/%d",b[0],width,b[1].tupleof);}
The output is exactly the same as the Python entry.
Thanks Rudy Velthuis for theVelthuis.BigRationals library.
programBernoulli_numbers;{$APPTYPE CONSOLE}usesSystem.SysUtils,Velthuis.BigRationals;functionb(n:Integer):BigRational;beginvara:TArray<BigRational>;SetLength(a,n+1);forvarm:=0toHigh(a)dobegina[m]:=BigRational.Create(1,m+1);forvarj:=mdownto1dobegina[j-1]:=(a[j-1]-a[j])*j;end;end;Result:=a[0];end;beginforvarn:=0to60dobeginvarbb:=b(n);ifbb.Numerator.BitLength>0thenwriteln(format('B(%2d) =%45s/%s',[n,bb.Numerator.ToString,bb.Denominator.ToString]));end;readln;end.
This example isin need of improvement:
|
Only 'small' rationals are supported in EchoLisp, i.e numerator and demominator < 2^31. So, we create a class of 'large' rationals, supported by the bigint library, and then apply the magic formula.
(lib'bigint);; lerge numbers(lib'gloops);; classes(define-classRationalnull((a:initform#0)(b:initform#1)))(define-methodtostring(Rational)(lambda(r)(format"%50d / %d"r.ar.b)))(define-methodnormalize(Rational)(lambda(r);; divide a and b by gcd(let((g(gcdr.ar.b)))(set!r.a(/r.ag))(set!r.b(/r.bg))(when(<r.b0)(set!r.a(-r.a))(set!r.b(-r.b)));; denominator > 0r)))(define-methodinitialize(Rational)(lambda(r)(normalizer)))(define-methodadd(Rational)(lambda(rn);; + Rational any number(normalize(Rational(+(*(+#0n)r.b)r.a)r.b))))(define-methodadd(RationalRational)(lambda(rq);;; + Rational Rational(normalize(Rational(+(*r.aq.b)(*r.bq.a))(*r.bq.b)))))(define-methodsub(RationalRational)(lambda(rq)(normalize(Rational(-(*r.aq.b)(*r.bq.a))(*r.bq.b)))))(define-methodmul(RationalRational)(lambda(rq)(normalize(Rational(*r.aq.a)(*r.bq.b)))))(define-methodmul(Rational)(lambda(rn)(normalize(Rational(*r.a(+#0n))r.b))))(define-methoddiv(RationalRational)(lambda(rq)(normalize(Rational(*r.aq.b)(*r.bq.a)))))
;; Bernoulli numbers;; http://rosettacode.org/wiki/Bernoulli_numbers(defineA(make-vector1000))(define(Bn)(for((m(1+n)));; #1 creates a large integer(vector-set!Am(Rational#1(+#1m)))(for((j(in-rangem0-1)))(vector-set!A(1-j)(mul(sub(vector-refA(1-j))(vector-refAj))j))))(vector-refA0))(for((b(in-range0622)))(writelnb(Bb)))→01/121/64-1/3061/428-1/30105/6612-691/2730147/616-3617/5101843867/79820-174611/33022854513/13824-236364091/2730268553103/628-23749461029/870308615841276005/1432232-7709321041217/510342577687858367/636-26315271553053477373/1919190382929993913841559/640-261082718496449122051/13530421520097643918070802691/180644-27833269579301024235023/69046596451111593912163277961/28248-5609403368997817686249127547/4641050495057205241079648212477525/6652-801165718135489957347924991853/15905429149963634884862421418123812691/79856-2479392929313226753685415739663229/8705884483613348880041862046775994036021/35460-1215233140483755572040304994079820246041491/56786730(B1)→1/2
defmoduleBernoullidodefmoduleRationaldoimportKernel,except:[div:2]defstructnumerator:0,denominator:1defnew(numerator,denominator\\1)dosign=ifnumerator*denominator<0,do:-1,else:1{numerator,denominator}={abs(numerator),abs(denominator)}gcd=gcd(numerator,denominator)%Rational{numerator:sign*Kernel.div(numerator,gcd),denominator:Kernel.div(denominator,gcd)}enddefsub(a,b)donew(a.numerator*b.denominator-b.numerator*a.denominator,a.denominator*b.denominator)enddefmul(a,b)whenis_integer(a)donew(a*b.numerator,b.denominator)enddefpgcd(a,0),do:adefpgcd(a,b),do:gcd(b,rem(a,b))enddefnumbers(n)doStream.transform(0..n,{},fnm,acc->acc=Tuple.append(acc,Rational.new(1,m+1))ifm>0donew=Enum.reduce(m..1,acc,fnj,ar->put_elem(ar,j-1,Rational.mul(j,Rational.sub(elem(ar,j-1),elem(ar,j))))end){[elem(new,0)],new}else{[elem(acc,0)],acc}endend)|>Enum.to_listenddeftask(n\\61)dob_nums=numbers(n)width=Enum.map(b_nums,fnb->b.numerator|>to_string|>String.lengthend)|>Enum.maxformat='B(~2w) = ~#{width}w / ~w~n'Enum.with_index(b_nums)|>Enum.each(fn{b,i}->ifb.numerator!=0,do::io.fwriteformat,[i,b.numerator,b.denominator]end)endendBernoulli.task
B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
openMathNet.NumericsopenSystemopenSystem.Collections.GenericletcalculateBernoullin=letℚ(x)=BigRational.FromIntxletA=Array.init<BigRational>(n+1)(funx->ℚ(x+1))formin[1..n]doA.[m]<-ℚ(1)/(ℚ(m)+ℚ(1))forjin[m..(-1)..1]doA.[j-1]<-ℚ(j)*(A.[j-1]-A.[j])A.[0][<EntryPoint>]letmainargv=fornin[0..60]doletbernoulliNumber=calculateBernoullinmatchbernoulliNumber.Numerator.IsZerowith|false->letformatedString=String.Format("B({0, 2}) = {1, 44} / {2}",n,bernoulliNumber.Numerator,bernoulliNumber.Denominator)printfn"%s"formatedString|true->printf""0
B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
One could use the "bernoulli" word from the math.extras vocabulary as follows:
IN: scratchpad [ 0 1 1 "%2d : %d / %d\n" printf 1 -1 2 "%2d : %d / %d\n" printf 30 iota [ 1 + 2 * dup bernoulli [ numerator ] [ denominator ] bi "%2d : %d / %d\n" printf ] each ] time 0 : 1 / 1 1 : -1 / 2 2 : 1 / 6 4 : -1 / 30 6 : 1 / 42 8 : -1 / 3010 : 5 / 6612 : -691 / 273014 : 7 / 616 : -3617 / 51018 : 43867 / 79820 : -174611 / 33022 : 854513 / 13824 : -236364091 / 273026 : 8553103 / 628 : -23749461029 / 87030 : 8615841276005 / 1432232 : -7709321041217 / 51034 : 2577687858367 / 636 : -26315271553053477373 / 191919038 : 2929993913841559 / 640 : -261082718496449122051 / 1353042 : 1520097643918070802691 / 180644 : -27833269579301024235023 / 69046 : 596451111593912163277961 / 28248 : -5609403368997817686249127547 / 4641050 : 495057205241079648212477525 / 6652 : -801165718135489957347924991853 / 159054 : 29149963634884862421418123812691 / 79856 : -2479392929313226753685415739663229 / 87058 : 84483613348880041862046775994036021 / 35460 : -1215233140483755572040304994079820246041491 / 56786730Running time: 0.00489444 seconds
Alternatively a method described by Brent and Harvey (2011) in "Fast computation of Bernoulli, Tangent and Secant numbers"https://arxiv.org/pdf/1108.0286.pdf is shown.
:: bernoulli-numbers ( n -- ) n 1 + 0 <array> :> tab 1 1 tab set-nth 2 n [a,b] [| k | k 1 - dup tab nth * k tab set-nth ] each 2 n [a,b] [| k | k n [a,b] [| j | j tab nth j k - 2 + * j 1 - tab nth j k - * + j tab set-nth ] each ] each 1 :> s! 1 n [a,b] [| k | k 2 * dup 2^ dup 1 - * k tab nth swap / * s * k tab set-nth s -1 * s! ] each 0 1 1 "%2d : %d / %d\n" printf 1 -1 2 "%2d : %d / %d\n" printf 1 n [a,b] [| k | k 2 * k tab nth [ numerator ] [ denominator ] bi "%2d : %d / %d\n" printf ] each;
It gives the same result as the native implementation, but is slightly faster.
[ 30 bernoulli-numbers ] time...Running time: 0.004331652 seconds
Func Bern(m) = Sigma<k=0,m>[Sigma<v=0,k>[(-1)^v*Bin(k,v)*(v+1)^m/(k+1)]].;for i=0, 60 do b:=Bern(i); if b<>0 then !!(i,b) fi od;
0 1
1 1 / 22 1 / 64 -1 / 306 1 / 428 -1 / 3010 5 / 6612 -691 / 273014 7 / 616 -3617 / 51018 43867 / 79820 -174611 / 33022 854513 / 13824 -236364091 / 273026 8553103 / 628 -23749461029 / 87030 8615841276005 / 1432232 -7709321041217 / 51034 2577687858367 / 636 -26315271553053477373 / 191919038 2929993913841559 / 640 -261082718496449122051 / 1353042 1520097643918070802691 / 180644 -27833269579301024235023 / 69046 596451111593912163277961 / 28248 -5609403368997817686249127547 / 4641050 495057205241079648212477525 / 6652 -801165718135489957347924991853 / 159054 29149963634884862421418123812691 / 79856 -2479392929313226753685415739663229 / 87058 84483613348880041862046775994036021 / 35460 -1215233140483755572040304994079820246041491 / 56786730
' version 08-10-2016' compile with: fbc -s console' uses gmp#IncludeOnce"gmp.bi"#Definemax60DimAsLongnDimAsZStringPtrgmp_str:gmp_str=Allocate(1000)' 1000 charDimSharedAsMpq_ptrtmp,big_jtmp=Allocate(Len(__mpq_struct)):Mpq_init(tmp)big_j=Allocate(Len(__mpq_struct)):Mpq_init(big_j)DimSharedAsMpq_ptra(max),b(max)Forn=0TomaxA(n)=Allocate(Len(__mpq_struct)):Mpq_init(A(n))B(n)=Allocate(Len(__mpq_struct)):Mpq_init(B(n))NextFunctionBernoulli(nAsInteger)AsMpq_ptrDimAsLongm,jForm=0TonMpq_set_ui(A(m),1,m+1)Forj=mTo1Step-1Mpq_sub(tmp,A(j-1),A(j))Mpq_set_ui(big_j,j,1)'big_j = jMpq_mul(A(j-1),big_j,tmp)NextNextReturnA(0)EndFunction' ------=< MAIN >=------Forn=0TomaxMpq_set(B(n),Bernoulli(n))Mpq_get_str(gmp_str,10,B(n))If*gmp_str<>"0"ThenIf*gmp_str="1"Then*gmp_str="1/1"PrintUsing"B(##) = ";n;PrintSpace(45-InStr(*gmp_str,"/"));*gmp_strEndIfNext' empty keyboard bufferWhileInkey<>"":WendPrint:Print"hit any key to end program"SleepEnd
B( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
BernoulliNumber[n] :={ a = new array for m = 0 to n { a@m = 1/(m+1) for j = m to 1 step -1 a@(j-1) = j * (a@(j-1) - a@j) } return a@0}result = new arrayfor n=0 to 60{ b = BernoulliNumber[n] if b != 0 { [num,den] = numeratorDenominator[b] result.push[[n, num, "/", den]] }}println[formatTable[result, "right"]]0 1 / 1 1 1 / 2 2 1 / 6 4 -1 / 30 6 1 / 42 8 -1 / 3010 5 / 6612 -691 / 273014 7 / 616 -3617 / 51018 43867 / 79820 -174611 / 33022 854513 / 13824 -236364091 / 273026 8553103 / 628 -23749461029 / 87030 8615841276005 / 1432232 -7709321041217 / 51034 2577687858367 / 636 -26315271553053477373 / 191919038 2929993913841559 / 640 -261082718496449122051 / 1353042 1520097643918070802691 / 180644 -27833269579301024235023 / 69046 596451111593912163277961 / 28248 -5609403368997817686249127547 / 4641050 495057205241079648212477525 / 6652 -801165718135489957347924991853 / 159054 29149963634884862421418123812691 / 79856 -2479392929313226753685415739663229 / 87058 84483613348880041862046775994036021 / 35460 -1215233140483755572040304994079820246041491 / 56786730
FunL has pre-defined functionB in moduleintegers, which is defined as:
import integers.choosedef B( n ) = sum( 1/(k + 1)*sum((if 2|r then 1 else -1)*choose(k, r)*(r^n) | r <- 0..k) | k <- 0..n )for i <- 0..60 if i == 1 or 2|i printf( "B(%2d) = %s\n", i, B(i) )
B( 0) = 1B( 1) = -1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in itswebsite.
Inthis page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.
Solution. The following function reduces to the n-th Bernoulli number. It is a replica of the Akiyama–Tanigawa algorithm.
Test case. Showing the Bernoulli numbers B0 to B60
Using FutureBasic 7.0.35
FB GMP/MPFR frameworks and headers (external link)
//// Bernoulli Number// Using FutureBasic 7.0.35//// September 2025, R.W//include "gmp.incl"_kMaxN = 200 // max Bernoulli termsmpq_t tmpdim a( _kMaxN ) as mpq_t // workspace a(0.._kMaxN)// Bernoulli:// B_m = -(1 / (m+1)) * sum_{k=0..m-1} binomial(m+1, k) * B_klocal fn Bernoulli( rop as mpq_t, n as long ) long m, j if n < 0 then n = 0 if n > _kMaxN then n = _kMaxN for m = 1 to n + 1 mpq_set_si( a(m-1), 1, m ) // j = m-1..1 for j = m - 1 to 1 step -1 // a[j] = a[j+1] - a[j] with 0-based mapping mpq_sub( a(j-1), a(j), a(j-1) ) // multiply by j mpq_set_si( tmp, j, 1 ) mpq_mul( a(j-1), a(j-1), tmp ) next next mpq_set( rop, a(0) )end fn//-------------- Main ---------Window 1, @"Bernoulli Numbers",(0,0,500,500)CFStringRef cfnum, cfden, pLine, filler,sp,sp2Int termslong iterms = 60mpq_init( tmp )for i = 0 to terms mpq_init( a(i) )nextfiller = @" "mpq_t ropmpz_t num, denmpq_init( rop )mpz_init( num )mpz_init( den )pLine = @""for i = 0 to 60 step 2 if i = 2 then i-- if i = 3 then i-- fn Bernoulli( rop, i ) mpq_get_num( num, rop ) if i == 1 mpz_neg( num, num ) end if if fn mpz_cmp_si( num, 0 ) // nonzero? mpq_get_den( den, rop ) cfnum = fn mpz_cf(10, num) cfden = fn mpz_cf(10, den) if i < 10 sp = @" " else sp = @"" end if sp2 = left(filler,44 - len(cfnum)) pLine = concat(sp, @"B(",mid(str(i),1),@") = ",sp2 ,cfnum, @" / ", cfden) print@ pLine end ifnext// clean upfor i = 0 to terms mpq_clear( a(i) )nextmpq_clear( tmp )mpz_clear( den )mpz_clear( num )mpq_clear( rop )HandleEvents//B(0) = 1 / 1 B(1) = 1 / 2 B(2) = 1 / 6 B(4) = -1 / 30 B(6) = 1 / 42 B(8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
forainFiltered(List([0..60],n->[n,Bernoulli(n)]),x->x[2]<>0)doPrint(a,"\n");od;[0,1][1,-1/2][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,-23749461029/870][30,8615841276005/14322][32,-7709321041217/510][34,2577687858367/6][36,-26315271553053477373/1919190][38,2929993913841559/6][40,-261082718496449122051/13530][42,1520097643918070802691/1806][44,-27833269579301024235023/690][46,596451111593912163277961/282][48,-5609403368997817686249127547/46410][50,495057205241079648212477525/66][52,-801165718135489957347924991853/1590][54,29149963634884862421418123812691/798][56,-2479392929313226753685415739663229/870][58,84483613348880041862046775994036021/354][60,-1215233140483755572040304994079820246041491/56786730]
packagemainimport("fmt""math/big")funcb(nint)*big.Rat{varfbig.Rata:=make([]big.Rat,n+1)form:=rangea{a[m].SetFrac64(1,int64(m+1))forj:=m;j>=1;j--{d:=&a[j-1]d.Mul(f.SetInt64(int64(j)),d.Sub(d,&a[j]))}}returnf.Set(&a[0])}funcmain(){forn:=0;n<=60;n++{ifb:=b(n);b.Num().BitLen()>0{fmt.Printf("B(%2d) =%45s/%s\n",n,b.Num(),b.Denom())}}}
B( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
This program works as a command line utility, that reads from stdin the number of elements to compute (default 60) and prints them in stdout. The implementation of the algorithm is in the function bernoullis. The rest is for printing the results.
importData.RatioimportSystem.Environmentmain=getArgs>>=printM.defaultArgwheredefaultArgas=ifnullasthen60elseread(headas)printMm=mapM_(putStrLn.printP).takeWhile((<=m).fst).filter(\(_,b)->b/=0%1).zip[0..]$bernoullisprintP(i,r)="B("++showi++") = "++show(numeratorr)++"/"++show(denominatorr)bernoullis=maphead.iterate(ulli1).mapberno$enumFrom0wherebernoi=1%(i+1)ulli_[_]=[]ullii(x:y:xs)=(i%1)*(x-y):ulli(i+1)(y:xs)
B(0) = 1/1B(1) = 1/2B(2) = 1/6B(4) = -1/30B(6) = 1/42B(8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
importData.Bool(bool)importData.Ratio(Ratio,denominator,numerator,(%))-------------------- BERNOULLI NUMBERS -------------------bernouillis::Integer->[Rational]bernouillis=fmaphead.tail.scanlfaulhaber[].enumFromTo0faulhaber::[RatioInteger]->Integer->[RatioInteger]faulhaberrsn=(:)=<<(-)1.sum$zipWith((*).(n%))[2..]rs--------------------------- TEST -------------------------main::IO()main=doletxs=bernouillis60w=length$show(numerator(lastxs))putStrLn$fTable"Bernouillis from Faulhaber triangle:\n"(show.fst)(showRatiow.snd)id(filter((0/=).snd)$zip[0..]xs)------------------------ FORMATTING ----------------------fTable::String->(a->String)->(b->String)->(a->b)->[a]->StringfTablesxShowfxShowfxs=letw=maximum(length.xShow<$>xs)inunlines$s:fmap(((<>).rjustw' '.xShow)<*>((" -> "<>).fxShow.f))xsshowRatio::Int->Rational->StringshowRatiowr=letd=denominatorrinrjustw' '$show(numeratorr)<>bool[](" / "<>showd)(1/=d)rjust::Int->a->[a]->[a]rjustnc=drop.length<*>(replicatenc<>)
Bernouillis from Faulhaber triangle: 0 -> 1 1 -> 1 / 2 2 -> 1 / 6 4 -> -1 / 30 6 -> 1 / 42 8 -> -1 / 3010 -> 5 / 6612 -> -691 / 273014 -> 7 / 616 -> -3617 / 51018 -> 43867 / 79820 -> -174611 / 33022 -> 854513 / 13824 -> -236364091 / 273026 -> 8553103 / 628 -> -23749461029 / 87030 -> 8615841276005 / 1432232 -> -7709321041217 / 51034 -> 2577687858367 / 636 -> -26315271553053477373 / 191919038 -> 2929993913841559 / 640 -> -261082718496449122051 / 1353042 -> 1520097643918070802691 / 180644 -> -27833269579301024235023 / 69046 -> 596451111593912163277961 / 28248 -> -5609403368997817686249127547 / 4641050 -> 495057205241079648212477525 / 6652 -> -801165718135489957347924991853 / 159054 -> 29149963634884862421418123812691 / 79856 -> -2479392929313226753685415739663229 / 87058 -> 84483613348880041862046775994036021 / 35460 -> -1215233140483755572040304994079820246041491 / 56786730
Note that only the few invocationsprimcorec/definition in the last subsection are the implementation; the rest is about the correctness proof.
Adapted from the corresponding entry in theArchive of Formal Proofs.
Defines Bernoulli numbers as the sequence of numbers with the exponential generating function and then proves that the Akiyama–Tanigawa transform takes an ordinary generating function to the exponential generating function, from which it follows that applying it to the sequence with ordinary generating function gives exactly the Bernoulli numbers.
theoryMini_BernoulliimportsComplex_Main"HOL-Computational_Algebra.Computational_Algebra""HOL-Combinatorics.Stirling""HOL-Library.Stream""Coinductive.Coinductive_List""HOL-Library.Code_Target_Numeral""HOL-Library.Code_Lazy"beginsubsection‹Definition of Bernoulli numbers›(* We define Bernoulli numbers in the standard fashion as those numbers B_n whose exponential generating function function is X / (1 - e^(-X)): *)definitionbernoulli_fps::"rat fps"where"bernoulli_fps = fps_X / (1 - fps_exp (-1))"definitionbernoulli::"nat ⇒ rat"where"bernoulli n = fps_nth bernoulli_fps n * fact n"subsection‹Stirling numbers of the 2nd kind›(* We use the Stirling numbers from the Isabelle library and prove a few additional things about them. *)lemmaStirling_n_0:"Stirling n 0 = (if n = 0 then 1 else 0)"by(casesn)simp_alllemmasum_Stirling_binomial:"Stirling (Suc n) (Suc m) = (∑i=0..n. Stirling i m * (n choose i))"proof-have"real (Stirling (Suc n) (Suc m)) = real (∑i = 0..n. Stirling i m * (n choose i))"proof(inductionnarbitrary:m)case(Sucnm)have"real (∑i = 0..Suc n. Stirling i m * (Suc n choose i)) = real (∑i = 0..n. Stirling (Suc i) m * (Suc n choose Suc i)) + real (Stirling 0 m)"by(substsum.atLeast0_atMost_Suc_shift)simp_allalsohave"real (∑i = 0..n. Stirling (Suc i) m * (Suc n choose Suc i)) = real (∑i = 0..n. (n choose i) * Stirling (Suc i) m) + real (∑i = 0..n. (n choose Suc i) * Stirling (Suc i) m)"by(simpadd:algebra_simpssum.distrib)alsohave"(∑i = 0..n. (n choose Suc i) * Stirling (Suc i) m) = (∑i = Suc 0..Suc n. (n choose i) * Stirling i m)"by(substsum.shift_bounds_cl_Suc_ivl)simp_allalsohave"… = (∑i = Suc 0..n. (n choose i) * Stirling i m)"by(introsum.mono_neutral_right)autoalsohave"… = real (∑i = 0..n. Stirling i m * (n choose i)) - real (Stirling 0 m)"by(simpadd:sum.atLeast_Suc_atMostmult_ac)alsohave"real (∑i = 0..n. Stirling i m * (n choose i)) = real (Stirling (Suc n) (Suc m))"by(ruleSuc.IH[symmetric])alsohave"real (∑i = 0..n. (n choose i) * Stirling (Suc i) m) = real m * real (Stirling (Suc n) (Suc m)) + real (Stirling (Suc n) m)"by(casesm;(simponly:Suc.IH,simpadd:algebra_simpssum.distribsum_distrib_leftsum_distrib_right))alsohave"… + (real (Stirling (Suc n) (Suc m)) - real (Stirling 0 m)) + real (Stirling 0 m) = real (Suc m * Stirling (Suc n) (Suc m) + Stirling (Suc n) m)"by(simpadd:algebra_simpsdel:Stirling.simps)alsohave"Suc m * Stirling (Suc n) (Suc m) + Stirling (Suc n) m = Stirling (Suc (Suc n)) (Suc m)"by(ruleStirling.simps(4)[symmetric])finallyshow?case..qedsimp_allthus?thesisby(subst(asm)of_nat_eq_iff)qed(* The exponential generating function of the Stirling numbers (of the 2nd kind) with fixed second argument $m$ is S_m(X) = (e^X - 1)^m / m!. We use this fact to derive a closed form for them. *)lemmaStirling_fps_aux:"(fps_exp 1 - 1) ^ m $ n * fact n = (fact m * of_nat (Stirling n m) :: 'a :: field_char_0)"proof(inductionmarbitrary:n)case0thus?caseby(simpadd:Stirling_n_0)nextcase(Sucmn)show?caseproof(casesn)case0thus?thesisbysimpnextcase(Sucn')hence"(fps_exp 1 - 1 :: 'a fps) ^ Suc m $ n * fact n = fps_deriv ((fps_exp 1 - 1) ^ Suc m) $ n' * fact n'"by(simp_alladd:algebra_simpsdel:power_Suc)alsohave"fps_deriv ((fps_exp 1 - 1 :: 'a fps) ^ Suc m) = fps_const (of_nat (Suc m)) * ((fps_exp 1 - 1) ^ m * fps_exp 1)"by(substfps_deriv_power)simp_allalsohave"… $ n' * fact n' = of_nat (Suc m) * ((∑i = 0..n'. (fps_exp 1 - 1) ^ m $ i / fact (n' - i)) * fact n')"unfoldingfps_mult_left_const_nthby(simpadd:fps_mult_nthSuc.IHsum_distrib_rightdel:of_nat_Suc)alsohave"(∑i = 0..n'. (fps_exp 1 - 1 :: 'a fps) ^ m $ i / fact (n' - i)) * fact n' = (∑i = 0..n'. (fps_exp 1 - 1) ^ m $ i * fact n' / fact (n' - i))"by(substsum_distrib_right,rulesum.cong)(simp_alladd:divide_simps)alsohave"… = (∑i = 0..n'. (fps_exp 1 - 1) ^ m $ i * fact i * of_nat (n' choose i))"by(introsum.congrefl)(simp_alladd:binomial_fact)alsohave"… = (∑i = 0..n'. fact m * of_nat (Stirling i m) * of_nat (n' choose i))"by(simponly:Suc.IH)alsohave"of_nat (Suc m) * … = (fact (Suc m) :: 'a) * (∑i = 0..n'. of_nat (Stirling i m) * of_nat (n' choose i))"(is"_ = _ * ?S")by(simpadd:sum_distrib_leftsum_distrib_rightmult_acdel:of_nat_Suc)alsohave"?S = of_nat (Stirling (Suc n') (Suc m))"by(substsum_Stirling_binomial)simpalsohave"Suc n' = n"by(simpadd:Suc)finallyshow?thesis.qedqedtheoremStirling_closed_form:"(of_nat (Stirling n k) :: 'a :: field_char_0) = (∑j≤k. (-1)^(k - j) * of_nat (k choose j) * of_nat j ^ n) / fact k"proof-have"(fps_exp 1 - 1 :: 'a fps) = (fps_exp 1 + (-1))"bysimpalsohave"… ^ k = (∑j≤k. of_nat (k choose j) * fps_exp 1 ^ j * (- 1) ^ (k - j))"unfoldingbinomial_ring..alsohave"… = (∑j≤k. fps_const ((-1) ^ (k - j) * of_nat (k choose j)) * fps_exp (of_nat j))"by(simpadd:fps_const_mult[symmetric]fps_const_power[symmetric]fps_const_neg[symmetric]mult_acfps_of_natfps_exp_power_multdel:fps_const_multfps_const_powerfps_const_neg)alsohave"fps_nth … n = (∑j≤k. (- 1) ^ (k - j) * of_nat (k choose j) * of_nat j ^ n) / fact n"by(simpadd:fps_sum_nthsum_divide_distrib)alsohave"… * fact n = (∑j≤k. (- 1) ^ (k - j) * of_nat (k choose j) * of_nat j ^ n)"bysimpalsonoteStirling_fps_aux[ofkn]finallyshow?thesisby(simpadd:atLeast0AtMostfield_simps)qed(* We now define a somewhat ad-hoc operator formal power series: XD' maps a formal power series A(X) to the series (X - 1) A'(X). The relevance of this operator to us is that the n-fold iteration of this operator is related to Stirling numbers.*)definitionfps_XD'::"'a :: field_char_0 fps ⇒ 'a fps"where"fps_XD' = (λb. (fps_X - 1) * fps_deriv b)"lemmafps_XD'_0[simp]:"fps_XD' 0 = 0"andfps_XD'_1[simp]:"fps_XD' 1 = 0"andfps_XD'_add[simp]:"fps_XD' (b + c) = fps_XD' b + fps_XD' c"andfps_XD'_mult:"fps_XD' (b * c) = fps_XD' b * c + b * fps_XD' c"by(simp_alladd:fps_XD'_defalgebra_simps)lemmafps_XD'_power:"fps_XD' (b ^ n) = of_nat n * b ^ (n - 1) * fps_XD' b"by(inductionn)(simp_alladd:algebra_simpsfps_XD'_multpower_eq_if)lemmafps_XD'_sum:"fps_XD' (sum f A) = sum (λx. fps_XD' (f x)) A"by(inductionArule:infinite_finite_induct)simp_alllemmafps_XD'_funpow:defines"S ≡ λn i. fps_const (of_nat (Stirling n i))"shows"(fps_XD' ^^ n) H = (∑m≤n. S n m * (fps_X - 1) ^ m * (fps_deriv ^^ m) H)"proof(inductionnarbitrary:H)case0thus?caseby(simpadd:S_def)nextcase(SucnH)defineGwhere"G = (fps_X - 1 :: 'a fps)"have"(∑m≤Suc n. S (Suc n) m * G ^ m * (fps_deriv ^^ m) H) = (∑i≤n. of_nat (Suc i) * S n (Suc i) * G ^ Suc i * (fps_deriv ^^ Suc i) H) + (∑i≤n. S n i * G ^ Suc i * (fps_deriv ^^ Suc i) H)"(is"_ = sum (λi. ?f (Suc i)) … + ?S2")by(substsum.atMost_Suc_shift)(simp_alladd:sum.distribalgebra_simpsfps_of_natS_deffps_const_add[symmetric]fps_const_mult[symmetric]del:fps_const_addfps_const_mult)alsohave"sum (λi. ?f (Suc i)) {..n} = sum (λi. ?f (Suc i)) {..<n}"by(introsum.mono_neutral_right)(autosimp:S_def)alsohave"… = ?f 0 + …"bysimpalsohave"… = sum ?f {..n}"by(substsum.atMost_shift[symmetric])simp_allalsohave"… + ?S2 = (∑x≤n. fps_XD' (S n x * G ^ x * (fps_deriv ^^ x) H))"unfoldingsum.distrib[symmetric]proof(rulesum.cong,goal_cases)case(2i)thus?caseunfoldingfps_XD'_multfps_XD'_powerby(casesi)(autosimp:fps_XD'_multalgebra_simpsof_nat_diffS_deffps_XD'_defG_def)qedsimp_allalsohave"… = (fps_XD' ^^ Suc n) H"by(simpadd:Suc.IHfps_XD'_sumG_def)finallyshow?caseby(simpadd:G_def)qedsubsection‹The Akiyama–Tanigawa transform›(* a single step in the Akiyama–Tanigawa matrix, i.e. how to get the next line from the current one *)definitionAT_step::"(nat ⇒ 'a :: field_char_0) ⇒ nat ⇒ 'a"where"AT_step f = (λk. of_nat (k+1) * (f k - f (k+1)))"definitionAT_matrix::"(nat ⇒ 'a :: field_char_0) ⇒ nat ⇒ nat ⇒ 'a"where"AT_matrix f n = (AT_step ^^ n) f"(* The following describes the ordinary generating function of the n-th row in the AT matrix. *)definitionAT_fps::"(nat ⇒ 'a :: field_char_0) ⇒ nat ⇒ 'a fps"where"AT_fps f n = (fps_X - 1) * Abs_fps (AT_matrix f n)"lemmaAT_fps_Suc:"AT_fps f (Suc n) = (fps_X - 1) * fps_deriv (AT_fps f n)"proof(rulefps_ext)fixm::natshow"AT_fps f (Suc n) $ m = ((fps_X - 1) * fps_deriv (AT_fps f n)) $ m"by(casesm)(simp_alladd:AT_fps_defAT_matrix_defAT_step_deffps_deriv_defalgebra_simps)qedlemmaAT_fps_altdef:"AT_fps f n = (∑m≤n. fps_const (of_nat (Stirling n m)) * (fps_X - 1)^m * (fps_deriv ^^ m) (AT_fps f 0))"proof-have"AT_fps f n = (fps_XD' ^^ n) (AT_fps f 0)"by(inductionn)(simp_alladd:AT_fps_Sucfps_XD'_def)alsohave"… = (∑m≤n. fps_const (of_nat (Stirling n m)) * (fps_X - 1) ^ m * (fps_deriv ^^ m) (AT_fps f 0))"by(rulefps_XD'_funpow)finallyshow?thesis.qedlemmaAT_fps_0_nth:"AT_fps f 0 $ n = (if n = 0 then -f 0 else f (n - 1) - f n)"by(simpadd:AT_fps_defAT_matrix_defalgebra_simps)(* The following gives a closed form for the first column of the AT matrix, i.e. the result of the transform. *)lemmaAT_matrix_firstcol:"AT_matrix f n 0 = (∑k≤n. (-1) ^ k * fact k * of_nat (Stirling (Suc n) (Suc k)) * f k)"proof(cases"n = 0")caseFalsehave"AT_matrix f n 0 = -(AT_fps f n $ 0)"by(simpadd:AT_fps_def)alsohave"AT_fps f n $ 0 = (∑k≤n. of_nat (Stirling n k) * (- 1) ^ k * (fact k * AT_fps f 0 $ k))"by(substAT_fps_altdef)(simpadd:fps_sum_nthfps_nth_power_0fps_0th_higher_deriv)alsohave"… = (∑k≤n. of_nat (Stirling n k) * (- 1) ^ k * (fact k * (f (k - 1) - f k)))"usingFalseby(introsum.congrefl)(autosimp:Stirling_n_0AT_fps_0_nth)alsohave"… = (∑k≤n. fact k * (of_nat (Stirling n k) * (- 1) ^ k) * f (k - 1)) - (∑k≤n. fact k * (of_nat (Stirling n k) * (- 1) ^ k) * f k)"(is"_ = sum ?f _ - ?S2")by(simpadd:sum_subtractfalgebra_simps)alsofromFalsehave"sum ?f {..n} = sum ?f {0<..n}"by(introsum.mono_neutral_right)(autosimp:Stirling_n_0)alsohave"… = sum ?f {0<..Suc n}"by(introsum.mono_neutral_left)autoalsohave"{0<..Suc n} = {Suc 0..Suc n}"byautoalsohave"sum ?f … = sum (λn. ?f (Suc n)) {0..n}"by(substsum.atLeast_Suc_atMost_Suc_shift)simp_allalsohave"{0..n} = {..n}"byautoalsohave"sum (λn. ?f (Suc n)) … - ?S2 = (∑k≤n. -((-1)^k * fact k * of_nat (Stirling (Suc n) (Suc k)) * f k))"by(substsum_subtractf[symmetric],introsum.cong)(simp_alladd:algebra_simps)alsohave"-… = (∑k≤n. ((-1)^k * fact k * of_nat (Stirling (Suc n) (Suc k)) * f k))"by(simpadd:sum_negf)finallyshow?thesis.qed(simp_alladd:AT_matrix_def)(* The following theorem relates the exponential generating function B(X) of the transformed sequence to the ordinary generating function A(X) of the original sequence. Namely: B(X) = e^X A(1 - e^X). *)theoremAT_fps:"Abs_fps (λn. AT_matrix f n 0 / fact n) = fps_exp 1 * fps_compose (Abs_fps f) (1 - fps_exp 1)"proof(rulefps_ext)fixn::nathave"(fps_const (fact n) * (fps_compose (Abs_fps (λn. AT_matrix f 0 n)) (1 - fps_exp 1) * fps_exp 1)) $ n = (∑m≤n. ∑k≤m. (1 - fps_exp 1) ^ k $ m * fact n / fact (n - m) * f k)"unfoldingfps_mult_left_const_nthby(simpadd:fps_times_deffps_compose_defAT_matrix_firstcolsum_Stirling_binomialfield_simpssum_distrib_leftsum_distrib_rightatLeast0AtMostAT_matrix_defdel:Stirling.simpsof_nat_Suc)alsohave"… = (∑m≤n. ∑k≤m. (-1)^k * fact k * of_nat (Stirling m k * (n choose m)) * f k)"proof(introsum.congrefl,goal_cases)case(1mk)have"(1 - fps_exp 1 :: 'a fps) ^ k = (-fps_exp 1 + 1 :: 'a fps) ^ k"bysimpalsohave"… = (∑i≤k. of_nat (k choose i) * (-1) ^ i * fps_exp (of_nat i))"by(substbinomial_ring)(simpadd:atLeast0AtMostpower_minus'fps_exp_power_multmult.assoc)alsohave"… = (∑i≤k. fps_const (of_nat (k choose i) * (-1) ^ i) * fps_exp (of_nat i))"by(simpadd:fps_const_mult[symmetric]fps_of_natfps_const_power[symmetric]fps_const_neg[symmetric]del:fps_const_multfps_const_powerfps_const_neg)alsohave"… $ m = (∑i≤k. of_nat (k choose i) * (- 1) ^ i * of_nat i ^ m) / fact m"(is"_ = ?S / _")by(simpadd:fps_sum_nthsum_divide_distrib[symmetric])alsohave"?S = (-1) ^ k * (∑i≤k. (-1) ^ (k - i) * of_nat (k choose i) * of_nat i ^ m)"by(substsum_distrib_left,introsum.congrefl)(autosimp:minus_one_power_iff)alsohave"(∑i≤k. (-1) ^ (k - i) * of_nat (k choose i) * of_nat i ^ m) = of_nat (Stirling m k) * (fact k :: 'a)"by(substStirling_closed_form)(simp_alladd:field_simps)finallyhave*:"(1 - fps_exp 1 :: 'a fps) ^ k $ m * fact n / fact (n - m) = (- 1) ^ k * fact k * of_nat (Stirling m k) * of_nat (n choose m)"using1by(simpadd:binomial_factdel:of_nat_Suc)show?caseusing1by(subst*)simpqedalsohave"… = (∑m≤n. ∑k≤n. (- 1) ^ k * fact k * of_nat (Stirling m k * (n choose m)) * f k)"by(rulesum.cong[OFrefl],rulesum.mono_neutral_left)autoalsohave"… = (∑k≤n. ∑m≤n. (- 1) ^ k * fact k * of_nat (Stirling m k * (n choose m)) * f k)"by(rulesum.swap)alsohave"… = AT_matrix f n 0"by(simpadd:AT_matrix_firstcolsum_Stirling_binomialsum_distrib_leftsum_distrib_rightmult.assocatLeast0AtMostdel:Stirling.simps)finallyshow"Abs_fps (λn. AT_matrix f n 0 / fact n) $ n = (fps_exp 1 * (Abs_fps f oo 1 - fps_exp 1)) $ n"by(subst(asm)fps_mult_left_const_nth)(simpadd:field_simpsAT_matrix_defdel:of_nat_Suc)qed(* If we specialise this to the input sequence a(k) = 1 / (k+1), this has the ordinary generating function -ln(1 - X) / X, so the exponential generating function of the AT transform is e^X (-ln (1 - (1 - e^X)) / (1 - e^X)) = X / (1 - e^(-X)), which is exactly the exponential generating function of the Bernoulli numbers. *)corollarybernoulli_conv_AT:"bernoulli n = AT_matrix (λk. 1 / of_nat (k+1)) n 0"proof-definef::"nat ⇒ rat"where"f = (λn. 1 / of_nat (Suc n))"noteAT_fps[off]also{have"fps_ln 1 = fps_X * Abs_fps (λn. (-1)^n / of_nat (Suc n) :: rat)"by(introfps_ext)(simpdel:of_nat_Sucadd:fps_ln_def)hence"fps_ln 1 / fps_X = Abs_fps (λn. (-1)^n / of_nat (Suc n) :: rat)"by(metisfps_X_neq_zerononzero_mult_div_cancel_left)alsohave"fps_compose … (-fps_X) = Abs_fps f"by(simpadd:fps_compose_uminus'fps_eq_ifff_def)finallyhave"Abs_fps f = fps_compose (fps_ln 1 / fps_X) (-fps_X)"..alsohave"fps_ln 1 / fps_X oo - fps_X oo 1 - fps_exp (1::rat) = fps_ln 1 / fps_X oo fps_exp 1 - 1"by(substfps_compose_assoc[symmetric])(simp_alladd:fps_compose_uminus)alsohave"… = (fps_ln 1 oo fps_exp 1 - 1) / (fps_exp 1 - 1)"by(substfps_compose_divide_distrib)autoalsohave"… = fps_X / (fps_exp 1 - 1)"by(simpadd:fps_ln_fps_exp_invfps_inv_fps_exp_compose)finallyhave"Abs_fps f oo 1 - fps_exp 1 = fps_X / (fps_exp 1 - 1)".}alsohave"fps_exp (1::rat) - 1 = (1 - fps_exp (-1)) * fps_exp 1"by(simpadd:algebra_simpsfps_exp_add_mult[symmetric])alsohave"fps_exp 1 * (fps_X / …) = bernoulli_fps"unfoldingbernoulli_fps_defby(substdvd_div_mult2_eq)(autosimp:fps_dvd_iffintro!:subdegree_leI)finallyhave"Abs_fps (λn. AT_matrix f n 0 / fact n) = bernoulli_fps".hence"fps_nth (Abs_fps (λn. AT_matrix f n 0 / fact n)) n = fps_nth bernoulli_fps n"by(rulearg_cong)thus?thesisby(simpadd:fps_eq_ifff_defbernoulli_deffield_simps)qedsubsection‹Efficient code›(* Next, we implement the AT transform using infinite streams. We define the infinite AT matrix as a stream of streams of numbers and then define its leftmost column as the result of the transform. *)primcorecAT_impl_step::"nat ⇒ 'a :: field_char_0 stream ⇒ 'a stream"where"AT_impl_step m xs = of_nat m * (xs !! 0 - xs !! 1) ## AT_impl_step (Suc m) (stl xs)"primcorecAT_matrix_impl::"'a :: field_char_0 stream ⇒ 'a stream stream"where"AT_matrix_impl xs = xs ## AT_matrix_impl (AT_impl_step 1 xs)"definitionAT_impl::"'a :: field_char_0 stream ⇒ 'a :: field_char_0 stream"where"AT_impl xs = smap shd (AT_matrix_impl xs)"lemmasnth_AT_impl_step[simp]:"AT_impl_step m xs !! n = of_nat (m + n) * (xs !! n - xs !! (n+1))"by(inductionnarbitrary:mxs;substAT_impl_step.code)autolemmasnth_AT_matrix_impl[simp]:"AT_matrix_impl xs !! n !! k = AT_matrix (snth xs) n k"by(inductionnarbitrary:xsk;substAT_matrix_impl.code)(autosimp:snth_AT_impl_step[abs_def]AT_matrix_deffunpow_Suc_rightAT_step_defsimpdel:funpow.simpsof_nat_Suc)lemmasnth_AT_impl[simp]:"AT_impl xs !! n = AT_matrix (snth xs) n 0"by(simpadd:AT_impl_defflip:snth.simps(1))definitionbernoulli_stream::"rat stream"where"bernoulli_stream = AT_impl (smap (λi. 1 / of_nat i) (fromN 1))"theorembernoulli_stream_correct:"bernoulli_stream !! n = bernoulli n"by(simpadd:bernoulli_stream_defbernoulli_conv_ATsnth_smap[abs_def])end
code_lazy_typestreamcode_lazy_typellistprimcorecllist_of_stream::"'a stream ⇒ 'a llist"where"llist_of_stream xs = LCons (shd xs) (llist_of_stream (stl xs))"value"list_of (ltakeWhile (λ(i,_). i ≤ 60) (llist_of_stream (sfilter (λ(_, x). x ≠ 0) (szip (fromN 0) bernoulli_stream))))"(* Output"[(0, 1), (1, 1 / 2), (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, - (23749461029 / 870)), (30, 8615841276005 / 14322), (32, - (7709321041217 / 510)), (34, 2577687858367 / 6), (36, - (26315271553053477373 / 1919190)), (38, 2929993913841559 / 6), (40, - (261082718496449122051 / 13530)), (42, 1520097643918070802691 / 1806), (44, - (27833269579301024235023 / 690)), (46, 596451111593912163277961 / 282), (48, - (5609403368997817686249127547 / 46410)), (50, 495057205241079648212477525 / 66), (52, - (801165718135489957347924991853 / 1590)), (54, 29149963634884862421418123812691 / 798), (56, - (2479392929313226753685415739663229 / 870)), (58, 84483613348880041862046775994036021 / 354), (60, - (1215233140483755572040304994079820246041491 / 56786730))]" :: "(nat × rat) list"*)
The following works in both languages:
link"rational"proceduremain(args)limit:=integer(!args)|60everyb:=bernoulli(i:=0tolimit)doifb.numer>0thenwrite(right(i,3),": ",align(rat2str(b),60))endprocedurebernoulli(n)(A:=table(0))[0]:=rational(1,1,1)everym:=1tondo{A[m]:=rational(1,m+1,1)everyj:=mto1by-1doA[j-1]:=mpyrat(rational(j,1,1),subrat(A[j-1],A[j]))}returnA[0]endprocedurealign(r,n)returnrepl(" ",n-find("/",r))||rend
Sample run:
->bernoulli 60 0: (1/1) 1: (1/2) 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: (-23749461029/870) 30: (8615841276005/14322) 32: (-7709321041217/510) 34: (2577687858367/6) 36: (-26315271553053477373/1919190) 38: (2929993913841559/6) 40: (-261082718496449122051/13530) 42: (1520097643918070802691/1806) 44: (-27833269579301024235023/690) 46: (596451111593912163277961/282) 48: (-5609403368997817686249127547/46410) 50: (495057205241079648212477525/66) 52: (-801165718135489957347924991853/1590) 54: (29149963634884862421418123812691/798) 56: (-2479392929313226753685415739663229/870) 58: (84483613348880041862046775994036021/354) 60: (-1215233140483755572040304994079820246041491/56786730)->
Implementation:
SeeBernoulli Numbers Essay on the J wiki.
B=:{.&1%.(i.!])@>:@i.@x:
Task:
'B',.rplc&'r/_-'"1":(#~0~:{:"1)(i.,.B)61B01B1-1/2B21/6B4-1/30B61/42B8-1/30B105/66B12-691/2730B147/6B16-3617/510B1843867/798B20-174611/330B22854513/138B24-236364091/2730B268553103/6B28-23749461029/870B308615841276005/14322B32-7709321041217/510B342577687858367/6B36-26315271553053477373/1919190B382929993913841559/6B40-261082718496449122051/13530B421520097643918070802691/1806B44-27833269579301024235023/690B46596451111593912163277961/282B48-5609403368997817686249127547/46410B50495057205241079648212477525/66B52-801165718135489957347924991853/1590B5429149963634884862421418123812691/798B56-2479392929313226753685415739663229/870B5884483613348880041862046775994036021/354B60-1215233140483755572040304994079820246041491/56786730
importorg.apache.commons.math3.fraction.BigFraction;publicclassBernoulliNumbers{publicstaticvoidmain(String[]args){for(intn=0;n<=60;n++){BigFractionb=bernouilli(n);if(!b.equals(BigFraction.ZERO))System.out.printf("B(%-2d) = %-1s%n",n,b);}}staticBigFractionbernouilli(intn){BigFraction[]A=newBigFraction[n+1];for(intm=0;m<=n;m++){A[m]=newBigFraction(1,(m+1));for(intj=m;j>=1;j--)A[j-1]=(A[j-1].subtract(A[j])).multiply(newBigFraction(j));}returnA[0];}}
B(0 ) = 1B(1 ) = 1 / 2B(2 ) = 1 / 6B(4 ) = -1 / 30B(6 ) = 1 / 42B(8 ) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
classFraction{constructor(numerator,denominator=1){if(denominator===0)thrownewError("Denominator cannot be zero");// Convert to BigIntthis.numerator=BigInt(numerator);this.denominator=BigInt(denominator);// Handle negative denominatorif(this.denominator<0n){this.numerator=-this.numerator;this.denominator=-this.denominator;}// SimplifyconstgcdVal=this.gcd(this.abs(this.numerator),this.abs(this.denominator));this.numerator=this.numerator/gcdVal;this.denominator=this.denominator/gcdVal;}abs(n){returnn<0n?-n:n;}gcd(a,b){while(b!==0n){consttemp=b;b=a%b;a=temp;}returna;}subtract(other){constnum=this.numerator*other.denominator-other.numerator*this.denominator;constden=this.denominator*other.denominator;returnnewFraction(num,den);}multiply(other){constnum=this.numerator*other.numerator;constden=this.denominator*other.denominator;returnnewFraction(num,den);}equals(other){returnthis.numerator===other.numerator&&this.denominator===other.denominator;}toString(){if(this.denominator===1n){returnthis.numerator.toString();}return`${this.numerator} /${this.denominator}`;}}functionbernoulli(n){constA=newArray(n+1);for(letm=0;m<=n;m++){A[m]=newFraction(1,m+1);for(letj=m;j>=1;j--){A[j-1]=A[j-1].subtract(A[j]).multiply(newFraction(j));}}returnA[0];}// Main executionconstZERO=newFraction(0);for(letn=0;n<=60;n++){constb=bernoulli(n);if(!b.equals(ZERO)){console.log(`B(${n.toString().padEnd(2)}) =${b}`);}}
B(0 ) = 1B(1 ) = 1 / 2B(2 ) = 1 / 6B(4 ) = -1 / 30B(6 ) = 1 / 42B(8 ) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
This section uses the Akiyama–Tanigawa algorithm for the second Bernoulli numbers, Bn. Therefore, the sign of B(1) differs from the modern definition.
The implementation presented here is intended for use with a "BigInt" library that uses string representations of decimal integers.Such a library is atBigInt.jq.To make the code in this section self-contained, stubs for the "BigInt" operations are provided in the first subsection.
BigInt Stubs:
# def negate:# def lessOrEqual(x; y): # x <= y# def long_add(x;y): # x+y# def long_minus(x;y): # x-y# def long_multiply(x;y) # x*y# def long_divide(x;y): # x/y => [q,r]# def long_div(x;y) # integer division# def long_mod(x;y) # %# In all cases, x and y must be stringsdef negate: (- tonumber) | tostring;def lessOrEqual(num1; num2): (num1|tonumber) <= (num2|tonumber);def long_add(num1; num2): ((num1|tonumber) + (num2|tonumber)) | tostring;def long_minus(x;y): ((num1|tonumber) - (num2|tonumber)) | tostring;# multiply two decimal strings, which may be signed (+ or -)def long_multiply(num1; num2): ((num1|tonumber) * (num2|tonumber)) | tostring;# return [quotient, remainder] # 0/0 = 1; n/0 => errordef long_divide(xx;yy): # x/y => [q,r] imples x == (y * q) + r def ld(x;y): def abs: if . < 0 then -. else . end; (x|abs) as $x | (y|abs) as $y | (if (x >= 0 and y > 0) or (x < 0 and y < 0) then 1 else -1 end) as $sign | (if x >= 0 then 1 else -1 end) as $sx | [$sign * ($x / $y | floor), $sx * ($x % $y)]; ld( xx|tonumber; yy|tonumber) | map(tostring);def long_div(x;y): long_divide(x;y) | .[0];def long_mod(x;y): ((x|tonumber) % (y|tonumber)) | tostring;
Fractions:
# A fraction is represented by [numerator, denominator] in reduced form, with the sign on top# a and b should be BigInt; return a BigIntdef gcd(a; b): def long_abs: . as $in | if lessOrEqual("0"; $in) then $in else negate end; # subfunction rgcd expects [a,b] as input # i.e. a ~ .[0] and b ~ .[1] def rgcd: .[0] as $a | .[1] as $b | if $b == "0" then $a else [$b, long_mod($a ; $b ) ] | rgcd end; a as $a | b as $b | [$a,$b] | rgcd | long_abs ;def normalize: .[0] as $p | .[1] as $q | if $p == "0" then ["0", "1"] elif lessOrEqual($q ; "0") then [ ($p|negate), ($q|negate)] | normalize else gcd($p; $q) as $g | [ long_div($p;$g), long_div($q;$g) ] end ;# a and b should be fractions expressed in the form [p, q]def add(a; b): a as $a | b as $b | if $a[1] == "1" and $b[1] == "1" then [ long_add($a[0]; $b[0]) , "1"] elif $a[1] == $b[1] then [ long_add( $a[0]; $b[0]), $a[1] ] | normalize elif $a[0] == "0" then $b elif $b[0] == "0" then $a else [ long_add( long_multiply($a[0]; $b[1]) ; long_multiply($b[0]; $a[1])), long_multiply($a[1]; $b[1]) ] | normalize end ;# a and/or b may be BigInts, or [p,q] fractionsdef multiply(a; b): a as $a | b as $b | if ($a|type) == "string" and ($b|type) == "string" then [ long_multiply($a; $b), "1"] else if $a|type == "string" then [ long_multiply( $a; $b[0]), $b[1] ] elif $b|type == "string" then [ long_multiply( $b; $a[0]), $a[1] ] else [ long_multiply( $a[0]; $b[0]), long_multiply($a[1]; $b[1]) ] end | normalize end ;def minus(a; b): a as $a | b as $b | if $a == $b then ["0", "1"] else add($a; [ ($b[0]|negate), $b[1] ] ) end ;Bernoulli Numbers:
# Using the algorithm in the task description:def bernoulli(n): reduce range(0; n+1) as $m ( []; .[$m] = ["1", long_add($m|tostring; "1")] # i.e. 1 / ($m+1) | reduce ($m - range(0 ; $m)) as $j (.; .[$j-1] = multiply( [($j|tostring), "1"]; minus( .[$j-1] ; .[$j]) ) )) | .[0] # (which is Bn) ;
The task:
range(0;61)| if . % 2 == 0 or . == 1 then "\(.): \(bernoulli(.) )" else empty end
The following output was obtained using the previously mentioned BigInt library.
$jq-n-r-fBernoulli.jq0:["1","1"]1:["1","2"]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:["-23749461029","870"]30:["8615841276005","14322"]32:["-7709321041217","510"]34:["2577687858367","6"]36:["-26315271553053477373","1919190"]38:["2929993913841559","6"]40:["-261082718496449122051","13530"]42:["1520097643918070802691","1806"]44:["-27833269579301024235023","690"]46:["596451111593912163277961","282"]48:["-5609403368997817686249127547","46410"]50:["495057205241079648212477525","66"]52:["-801165718135489957347924991853","1590"]54:["29149963634884862421418123812691","798"]56:["-2479392929313226753685415739663229","870"]58:["84483613348880041862046775994036021","354"]60:["-1215233140483755572040304994079820246041491","56786730"]
functionbernoulli(n)A=Vector{Rational{BigInt}}(undef,n+1)form=0:nA[m+1]=1//(m+1)forj=m:-1:1A[j]=j*(A[j]-A[j+1])endendreturnA[1]endfunctiondisplay(n)B=map(bernoulli,0:n)pad=mapreduce(x->ndigits(numerator(x))+Int(x<0),max,B)argdigits=ndigits(n)fori=0:nifnumerator(B[i+1])&1==1println("B(",lpad(i,argdigits),") = ",lpad(numerator(B[i+1]),pad)," / ",denominator(B[i+1]))endendenddisplay(60)# Alternative: Following the comment in the Perl section it is much more efficient# to compute the list of numbers instead of one number after the other.functionBernoulliList(len)A=Vector{Rational{BigInt}}(undef,len+1)B=similar(A)fornin0:lenA[n+1]=1//(n+1)forj=n:-1:1A[j]=j*(A[j]-A[j+1])endB[n+1]=A[1]endreturnBendfor(n,b)inenumerate(BernoulliList(60))isodd(numerator(b))&&println("B($(n-1)) =$b")end
Produces virtually the same output as the Python version.
importorg.apache.commons.math3.fraction.BigFractionobjectBernoulli{operatorfuninvoke(n:Int):BigFraction{valA=Array(n+1,init)for(min0..n)for(jinmdownTo1)A[j-1]=A[j-1].subtract(A[j]).multiply(integers[j])returnA.first()}valmax=60privatevalinit={m:Int->BigFraction(1,m+1)}privatevalintegers=Array(max+1,{m:Int->BigFraction(m)})}funmain(args:Array<String>){for(nin0..Bernoulli.max)if(n%2==0||n==1)System.out.printf("B(%-2d) = %-1s%n",n,Bernoulli(n))}
Produces virtually the same output as the Java version.
LuaJIT version with FFI and GMP library
#!/usr/bin/env luajitlocalgmp=require'gmp'('libgmp')localffi=require'ffi'localmpz,mpq=gmp.types.z,gmp.types.qlocalfunctionmpq_for(buf,op,n)fori=0,n-1doop(buf[i])endendlocalfunctionbernoulli(rop,n)locala=ffi.new("mpq_t[?]",n+1)mpq_for(a,gmp.q_init,n+1)form=0,ndogmp.q_set_ui(a[m],1,m+1)forj=m,1,-1dogmp.q_sub(a[j-1],a[j],a[j-1])gmp.q_set_ui(rop,j,1)gmp.q_mul(a[j-1],a[j-1],rop)endendgmp.q_set(rop,a[0])mpq_for(a,gmp.q_clear,n+1)enddo--MAINlocalrop=mpq()localn,d=mpz(),mpz()gmp.q_init(rop)gmp.z_inits(n,d)localto=arg[1]andtonumber(arg[1])or60localfrom=arg[2]andtonumber(arg[2])or0iffrom~=0thento,from=from,toendfori=from,todobernoulli(rop,i)ifgmp.q_cmp_ui(rop,0,1)~=0thengmp.q_get_num(n,rop)gmp.q_get_den(d,rop)gmp.printf("B(%-2g) = %44Zd / %Zd\n",i,n,d)endendgmp.z_clears(n,d)gmp.q_clear(rop)end
> time ./bernoulli_gmp.lua B(0 ) = 1 / 1B(1 ) = -1 / 2B(2 ) = 1 / 6B(4 ) = -1 / 30B(6 ) = 1 / 42B(8 ) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730./bernoulli_gmp.lua 0,02s user 0,00s system 97% cpu 0,022 total
Time compare: Python 0.591 sec, C 0.023 sec, Lua 0.022-0.025
Using a subset of Class RationalMaximum B(28) (if we use double and not decimal then we get lower maximum: B(20))
moduleBernoulli_Numbers{ClassRational{numeratorasdecimal,denominatorasdecimalgcd=lambda->0lcm=lambda->0operator"+"{Readldenom=.lcm(l.denominator,.denominator).numerator<=denom/l.denominator*l.numerator+denom/.denominator*.numeratorif.numerator==0thendenom=1.denominator<=denom}OperatorUnary{.numerator-!}Operator"-"{ReadlCallOperator"+",-l}GroupReal{value{linkparentnumerator,denominatorton,d=n/d}}GroupToString${value{linkparentnumerator,denominatorton,d=Str$(n)+"/"+Str$(d,"")}}class:ModuleRational(.numerator,.denominator){if.denominator=0then.denominator<=1whilefrac(.numerator)<>0{.numerator*=10@.denominator*=10@}sgn=Sgn(.numerator)*Sgn(.denominator).denominator<=abs(.denominator).numerator<=abs(.numerator)*sgngcd1=lambda(aasdecimal,basdecimal)->{ifa<bthenswapa,bg=amodbwhileg{a=b:b=g:g=amodb}=abs(b)}gdcval=gcd1(abs(.numerator),.denominator)ifgdcval<.denominatorandgdcval<>0then.denominator/=gdcval.numerator/=gdcvalendif.gcd<=gcd1.lcm<=lambdagcd=gcd1(aasdecimal,basdecimal)->{=a/gcd(a,b)*b}}}Open"out.txt"foroutputas#fr1=Rational(1,1)b=0nextBern()b=1nextBern()forb=2to28step2nextBern()nextclose#fsubnextBern()localm=@bernoulli(b)Print#F,format$("B({0::-2})={1}",b,m.tostring$)endsubfunctionbernoulli(nasdecimal)Locala(1ton+1)=Rational(1,1),zlocaldecimalm,jform=0tona(m+1)=Rational(1,m+1)ifm<1thencontinueforj=mto1z=a(J)-a(j+1)+r1a(j)=znn=z.gcd(j,z.denominator)a(j).numerator*=j/nna(j).denominator/=nnnextfora(1){this=rational(.numerator,.denominator)}next=a(1)endfunction}Bernoulli_Numbers
B( 0)= 1/1B( 1)= 3/2B( 2)= 1/6B( 4)=-1/30B( 6)= 1/42B( 8)=-1/30B(10)= 5/66B(12)=-691/2730B(14)= 7/6B(16)=-3617/510B(18)= 43867/798B(20)=-174611/330B(22)= 854513/138B(24)=-236364091/2730B(26)= 8553103/6B(28)=-23749461029/870
print(select(n->n[2]<>0,[seq([n,bernoulli(n,1)],n=0..60)]));
[[0, 1], [1, 1/2], [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, -23749461029/870], [30, 8615841276005/14322], [32, -7709321041217/510], [34, 2577687858367/6], [36, -26315271553053477373/1919190], [38, 2929993913841559/6], [40, -261082718496449122051/13530], [42, 1520097643918070802691/1806], [44, -27833269579301024235023/690], [46, 596451111593912163277961/282], [48, -5609403368997817686249127547/46410], [50, 495057205241079648212477525/66], [52, -801165718135489957347924991853/1590], [54, 29149963634884862421418123812691/798], [56, -2479392929313226753685415739663229/870], [58, 84483613348880041862046775994036021/354], [60, -1215233140483755572040304994079820246041491/56786730]]
Mathematica has no native way for starting an array at index 0. I therefore had to build the array from 1 to n+1 instead of from 0 to n, adjusting the formula accordingly.
bernoulli[n_]:=Module[{a=ConstantArray[0,n+2]},Do[a[[m]]=1/m;If[m==1&&a[[1]]!=0,Print[{m-1,a[[1]]}]];Do[a[[j-1]]=(j-1)*(a[[j-1]]-a[[j]]);If[j==2&&a[[1]]!=0,Print[{m-1,a[[1]]}]];,{j,m,2,-1}];,{m,1,n+1}];]bernoulli[60]
{0,1}{1,1/2}{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,-(23749461029/870)}{30,8615841276005/14322}{32,-(7709321041217/510)}{34,2577687858367/6}{36,-(26315271553053477373/1919190)}{38,2929993913841559/6}{40,-(261082718496449122051/13530)}{42,1520097643918070802691/1806}{44,-(27833269579301024235023/690)}{46,596451111593912163277961/282}{48,-(5609403368997817686249127547/46410)}{50,495057205241079648212477525/66}{52,-(801165718135489957347924991853/1590)}{54,29149963634884862421418123812691/798}{56,-(2479392929313226753685415739663229/870)}{58,84483613348880041862046775994036021/354}{60,-(1215233140483755572040304994079820246041491/56786730)}Or, it's permissible to use the native Bernoulli number function instead of being forced to use the specified algorithm, we very simply have:
(Note from task's author: nobody is forced to use any specific algorithm, the one shown is just a suggestion.)
Table[{i,BernoulliB[i]},{i,0,60}];Select[%,#[[2]]!=0&]//TableForm
011-(1/2)21/64-(1/30)61/428-(1/30)105/6612-(691/2730)147/616-(3617/510)1843867/79820-(174611/330)22854513/13824-(236364091/2730)268553103/628-(23749461029/870)308615841276005/1432232-(7709321041217/510)342577687858367/636-(26315271553053477373/1919190)382929993913841559/640-(261082718496449122051/13530)421520097643918070802691/180644-(27833269579301024235023/690)46596451111593912163277961/28248-(5609403368997817686249127547/46410)50495057205241079648212477525/6652-(801165718135489957347924991853/1590)5429149963634884862421418123812691/79856-(2479392929313226753685415739663229/870)5884483613348880041862046775994036021/35460-(1215233140483755572040304994079820246041491/56786730)
Using built-in function bern
block(makelist([sconcat("B","(",i,")","="),bern(i)],i,0,60),sublist(%%,lambda([x],x[2]#0)),table_form(%%))
matrix(["B(0)=",1],["B(1)=",-1/2],["B(2)=",1/6],["B(4)=",-1/30],["B(6)=",1/42],["B(8)=",-1/30],["B(10)=",5/66],["B(12)=",-691/2730],["B(14)=",7/6],["B(16)=",-3617/510],["B(18)=",43867/798],["B(20)=",-174611/330],["B(22)=",854513/138],["B(24)=",-236364091/2730],["B(26)=",8553103/6],["B(28)=",-23749461029/870],["B(30)=",8615841276005/14322],["B(32)=",-7709321041217/510],["B(34)=",2577687858367/6],["B(36)=",-26315271553053477373/1919190],["B(38)=",2929993913841559/6],["B(40)=",-261082718496449122051/13530],["B(42)=",1520097643918070802691/1806],["B(44)=",-27833269579301024235023/690],["B(46)=",596451111593912163277961/282],["B(48)=",-5609403368997817686249127547/46410],["B(50)=",495057205241079648212477525/66],["B(52)=",-801165718135489957347924991853/1590],["B(54)=",29149963634884862421418123812691/798],["B(56)=",-2479392929313226753685415739663229/870],["B(58)=",84483613348880041862046775994036021/354],["B(60)=",-1215233140483755572040304994079820246041491/56786730])
importbignumimportstrformatconstLim=60#---------------------------------------------------------------------------------------------------procbernoulli(n:Natural):Rat=## Compute a Bernoulli number using Akiyama–Tanigawa algorithm.vara=newSeq[Rat](n+1)formin0..n:a[m]=newRat(1,m+1)forjincountdown(m,1):a[j-1]=j*(a[j]-a[j-1])result=a[0]#———————————————————————————————————————————————————————————————————————————————————————————————————typeInfo=tuplen:int# Number index in Bernoulli sequence.val:Rat# Bernoulli number.varvalues:seq[Info]# List of values as Info tuples.varmaxLen=-1# Maximum length.# First step: compute the values and prepare for display.fornin0..Lim:# Compute value.ifn!=1and(nand1)==1:continue# Ignore odd "n" except 1.letb=bernoulli(n)# Check numerator length.letlen=($b.num).leniflen>maxLen:maxLen=len# Store information for next step.values.add((n,b))# Second step: display the values with '/' aligned.for(n,b)invalues:lets=fmt"{($b.num).alignString(maxLen, '>')} / {b.denom}"echofmt"{n:2}: {s}"
0: 1 / 1 1: -1 / 2 2: 1 / 6 4: -1 / 30 6: 1 / 42 8: -1 / 3010: 5 / 6612: -691 / 273014: 7 / 616: -3617 / 51018: 43867 / 79820: -174611 / 33022: 854513 / 13824: -236364091 / 273026: 8553103 / 628: -23749461029 / 87030: 8615841276005 / 1432232: -7709321041217 / 51034: 2577687858367 / 636: -26315271553053477373 / 191919038: 2929993913841559 / 640: -261082718496449122051 / 1353042: 1520097643918070802691 / 180644: -27833269579301024235023 / 69046: 596451111593912163277961 / 28248: -5609403368997817686249127547 / 4641050: 495057205241079648212477525 / 6652: -801165718135489957347924991853 / 159054: 29149963634884862421418123812691 / 79856: -2479392929313226753685415739663229 / 87058: 84483613348880041862046775994036021 / 35460: -1215233140483755572040304994079820246041491 / 56786730
for(n=0,60,t=bernfrac(n);if(t,print(n" "t)))
0 11 -1/22 1/64 -1/306 1/428 -1/3010 5/6612 -691/273014 7/616 -3617/51018 43867/79820 -174611/33022 854513/13824 -236364091/273026 8553103/628 -23749461029/87030 8615841276005/1432232 -7709321041217/51034 2577687858367/636 -26315271553053477373/191919038 2929993913841559/640 -261082718496449122051/1353042 1520097643918070802691/180644 -27833269579301024235023/69046 596451111593912163277961/28248 -5609403368997817686249127547/4641050 495057205241079648212477525/6652 -801165718135489957347924991853/159054 29149963634884862421418123812691/79856 -2479392929313226753685415739663229/87058 84483613348880041862046775994036021/35460 -1215233140483755572040304994079820246041491/56786730
Tested with fpc 3.0.4
(* Taken from the 'Ada 99' project, https://marquisdegeek.com/code_ada99 *)programBernoulliForAda99;usesBigDecimalMath;{library for arbitary high precision BCD numbers}typeFraction=objectprivatenumerator,denominator:BigDecimal;publicprocedureassign(n,d:Int64);proceduresubtract(rhs:Fraction);proceduremultiply(value:Int64);procedurereduce();procedurewriteOutput();end;functiongcd(a,b:BigDecimal):BigDecimal;beginif(b=0)thenbegingcd:=a;endelsebegingcd:=gcd(b,amodb);end;end;procedureFraction.writeOutput();varsign:char;beginsign:=' ';if(numerator<0)thensign:='-';if(denominator<0)thensign:='-';write(sign+BigDecimalToStr(abs(numerator)):45);write(' / ');write(BigDecimalToStr(abs(denominator)));end;procedureFraction.assign(n,d:Int64);beginnumerator:=n;denominator:=d;end;procedureFraction.subtract(rhs:Fraction);beginnumerator:=numerator*rhs.denominator;numerator:=numerator-(rhs.numerator*denominator);denominator:=denominator*rhs.denominator;end;procedureFraction.multiply(value:Int64);vartemp:BigDecimal;begintemp:=value;numerator:=numerator*temp;end;procedureFraction.reduce();vargcdResult:BigDecimal;begingcdResult:=gcd(numerator,denominator);beginnumerator:=numeratordivgcdResult;(* div is Int64 division *)denominator:=denominatordivgcdResult;(* could also use round(d/r) *)end;end;functioncalculateBernoulli(n:Int64):Fraction;varm,j:Int64;results:arrayofFraction;beginsetlength(results,60);{largest value 60}form:=0tondobeginresults[m].assign(1,m+1);forj:=mdownto1dobeginresults[j-1].subtract(results[j]);results[j-1].multiply(j);results[j-1].reduce();end;end;calculateBernoulli:=results[0];end;(* Main program starts here *)varb:Int64;result:Fraction;beginwriteln('Calculating Bernoulli numbers...');writeln('B( 0) : 1 / 1');forb:=1to60dobeginif(b<3)or((bmod2)=0)thenbeginresult:=calculateBernoulli(b);write('B(',b:2,')');write(' : ');result.writeOutput();writeln;end;end;end.
Calculating Bernoulli numbers...B( 0) : 1 / 1B( 1) : 1 / 2B( 2) : 1 / 6B( 4) : -1 / 30B( 6) : 1 / 42B( 8) : -1 / 30B(10) : 5 / 66B(12) : -691 / 2730B(14) : -7 / 6B(16) : -3617 / 510B(18) : 43867 / 798B(20) : -174611 / 330B(22) : 854513 / 138B(24) : -236364091 / 2730B(26) : 8553103 / 6B(28) : -23749461029 / 870B(30) : 8615841276005 / 14322B(32) : -7709321041217 / 510B(34) : 2577687858367 / 6B(36) : -26315271553053477373 / 1919190B(38) : 2929993913841559 / 6B(40) : -261082718496449122051 / 13530B(42) : 1520097643918070802691 / 1806B(44) : -27833269579301024235023 / 690B(46) : -596451111593912163277961 / 282B(48) : -5609403368997817686249127547 / 46410B(50) : 495057205241079648212477525 / 66B(52) : -801165718135489957347924991853 / 1590B(54) : 29149963634884862421418123812691 / 798B(56) : -2479392929313226753685415739663229 / 870B(58) : 84483613348880041862046775994036021 / 354B(60) : -1215233140483755572040304994079820246041491 / 56786730
The only thing in the suggested algorithm which depends on N is the number of times through the inner block. This means that all but the last iteration through the loop produce the exact same values of A.
Instead of doing the same calculations over and over again, I retain the A array until the final Bernoulli number is produced.
#!perlusestrict;usewarnings;useList::Utilqw(max);useMath::BigRat;my$one=Math::BigRat->new(1);subbernoulli_print{my@a;formy$m(0..60){push@a,$one/($m+1);formy$j(reverse1..$m){# This line:($a[$j-1]-=$a[$j])*=$j;# is a faster version of the following line:# $a[$j-1] = $j * ($a[$j-1] - $a[$j]);# since it avoids unnecessary object creation.}nextunless$a[0];printf"B(%2d) = %44s/%s\n",$m,$a[0]->parts;}}bernoulli_print();
The output is exactly the same as the Python entry.
We can also use modules for faster results. E.g.
usentheoryqw/bernfrac/;formy$n(0..60){my($num,$den)=bernfrac($n);printf"B(%2d) = %44s/%s\n",$n,$num,$denif$num!=0;}
with identical output. Or:
useMath::Pariqw/bernfrac/;formy$n(0..60){my($num,$den)=split"/",bernfrac($n);printf("B(%2d) = %44s/%s\n",$n,$num,$den||1)if$num!=0;}
with the difference being that Pari chooses = -½.
withjavascript_semanticsincludebuiltins/mpfr.eprocedurebernoulli(mpqrop,integern)sequencea=mpq_inits(n+1)form=1ton+1dompq_set_si(a[m],1,m)forj=m-1to1by-1dompq_sub(a[j],a[j+1],a[j])mpq_set_si(rop,j,1)mpq_mul(a[j],a[j],rop)endforendformpq_set(rop,a[1])a=mpq_free(a)endprocedurempqrop=mpq_init()mpzn=mpz_init(),d=mpz_init()fori=0to60dobernoulli(rop,i)ifmpq_cmp_si(rop,0,1)thenmpq_get_num(n,rop)mpq_get_den(d,rop)stringns=mpz_get_str(n),ds=mpz_get_str(d)printf(1,"B(%2d) = %44s / %s\n",{i,ns,ds})endifendfor{n,d}=mpz_free({n,d})rop=mpq_free(rop)
B( 0) = 1 / 1B( 1) = -1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
Brute force and method by Srinivasa Ramanujan.
(load "@lib/frac.l")(de fact (N) (cache '(NIL) N (if (=0 N) 1 (apply * (range 1 N))) ) )(de binomial (N K) (frac (/ (fact N) (* (fact (- N K)) (fact K)) ) 1 ) ) (de A (N M) (let Sum (0 . 1) (for X M (setq Sum (f+ Sum (f* (binomial (+ N 3) (- N (* X 6))) (berno (- N (* X 6)) ) ) ) ) ) Sum ) )(de berno (N) (cache '(NIL) N (cond ((=0 N) (1 . 1)) ((= 1 N) (-1 . 2)) ((bit? 1 N) (0 . 1)) (T (case (% N 6) (0 (f/ (f- (frac (+ N 3) 3) (A N (/ N 6)) ) (binomial (+ N 3) N) ) ) (2 (f/ (f- (frac (+ N 3) 3) (A N (/ (- N 2) 6)) ) (binomial (+ N 3) N) ) ) (4 (f/ (f- (f* (-1 . 1) (frac (+ N 3) 6)) (A N (/ (- N 4) 6)) ) (binomial (+ N 3) N) ) ) ) ) ) ) ) (de berno-brute (N) (cache '(NIL) N (let Sum (0 . 1) (cond ((=0 N) (1 . 1)) ((= 1 N) (-1 . 2)) ((bit? 1 N) (0 . 1)) (T (for (X 0 (> N X) (inc X)) (setq Sum (f+ Sum (f* (binomial (inc N) X) (berno-brute X)) ) ) ) (f/ (f* (-1 . 1) Sum) (binomial (inc N) N)) ) ) ) ) )(for (N 0 (> 62 N) (inc N)) (if (or (= N 1) (not (bit? 1 N))) (tab (2 4 -60) N " => " (sym (berno N))) ) )(for (N 0 (> 400 N) (inc N)) (test (berno N) (berno-brute N)) )(bye)
Bern: procedure options (main); /* 4 July 2014 */ declare i fixed binary; declare B complex fixed (31);Bernoulli: procedure (n) returns (complex fixed (31)); declare n fixed binary; declare anum(0:n) fixed (31), aden(0:n) fixed (31); declare (j, m) fixed; declare F fixed (31); do m = 0 to n; anum(m) = 1; aden(m) = m+1; do j = m to 1 by -1; anum(j-1) = j*( aden(j)*anum(j-1) - aden(j-1)*anum(j) ); aden(j-1) = ( aden(j-1) * aden(j) ); F = gcd(abs(anum(j-1)), abs(aden(j-1)) ); if F ^= 1 then do; anum(j-1) = anum(j-1) / F; aden(j-1) = aden(j-1) / F; end; end; end; return ( complex(anum(0), aden(0)) );end Bernoulli; do i = 0, 1, 2 to 36 by 2; /* 36 is upper limit imposed by hardware. */ B = Bernoulli(i); put skip edit ('B(' , trim(i) , ')=' , real(B) , '/' , trim(imag(B)) ) (3 A, column(10), F(32), 2 A); end;end Bern;The above uses GCD (see Rosetta Code) extended for 31-digit working.
Results obtained by this program are limited to the entries shown below due to the restrictions imposed by storing numbers in fixed decimal (31 digits).
B(0)= 1/1B(1)= 1/2B(2)= 1/6B(4)= -1/30B(6)= 1/42B(8)= -1/30B(10)= 5/66B(12)= -691/2730B(14)= 7/6B(16)= -3617/510B(18)= 43867/798B(20)= -174611/330B(22)= 854513/138B(24)= -236364091/2730B(26)= 8553103/6B(28)= -23749461029/870B(30)= 8615841276005/14322B(32)= -7709321041217/510B(34)= 2577687858367/6B(36)= -26315271553053477373/1919190
require"bignum"localfmt=require"fmt"localfunctionbernoulli(n)assert(n>=0,"Argument must be non-negative")locala=table.create(n+1)form=0,ndoa[m+1]=bigrat.recip(m+1)localj=mwhilej>=1doa[j]=(a[j]-a[j+1])*bigrat.int(j)j-=1endendreturnn!=1?a[1]:-a[1]--'first' Bernoulli numberendlocalzero=bigrat.int(0)forn=0,60dolocalb=bernoulli(n)ifb!=zerothenfmt.print("B(%2d) = %44s / %s",n,b:getNum(),b:getDen())endend
B( 0) = 1 / 1B( 1) = -1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
fromfractionsimportFractionasFrdefbernoulli(n):A=[0]*(n+1)forminrange(n+1):A[m]=Fr(1,m+1)forjinrange(m,0,-1):A[j-1]=j*(A[j-1]-A[j])returnA[0]# (which is Bn)bn=[(i,bernoulli(i))foriinrange(61)]bn=[(i,b)fori,binbnifb]width=max(len(str(b.numerator))fori,binbn)fori,binbn:print('B(%2i) =%*i/%i'%(i,width,b.numerator,b.denominator))
B( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
Using the optimization mentioned in the Perl entry to reduce intermediate calculations we create and use the generator bernoulli2():
defbernoulli2():A,m=[],0whileTrue:A.append(Fr(1,m+1))forjinrange(m,0,-1):A[j-1]=j*(A[j-1]-A[j])yieldA[0]# (which is Bm)m+=1bn2=[ixforixinzip(range(61),bernoulli2())]bn2=[(i,b)fori,binbn2ifb]width=max(len(str(b.numerator))fori,binbn2)fori,binbn2:print('B(%2i) =%*i/%i'%(i,width,b.numerator,b.denominator))
Output is exactly the same as before.
$ "bigrat.qky" loadfile [ 1+ ' [ [] ] over of swap times [ i^ 1+ n->v 1/v join swap i^ poke i^ times [ dup i 1+ peek do dip over swap i peek do v- i 1+ n->v v* join swap i poke ] ] 1 split drop do ] is bernoulli ( n --> n/d ) 61 times [ i^ bernoulli 2dup v0= iff 2drop else [ i^ 10 < if sp i^ echo sp vulgar$ char / over find 44 swap - times sp echo$ cr ] ]
0 1/1 1 -1/2 2 1/6 4 -1/30 6 1/42 8 -1/3010 5/6612 -691/273014 7/616 -3617/51018 43867/79820 -174611/33022 854513/13824 -236364091/273026 8553103/628 -23749461029/87030 8615841276005/1432232 -7709321041217/51034 2577687858367/636 -26315271553053477373/191919038 2929993913841559/640 -261082718496449122051/1353042 1520097643918070802691/180644 -27833269579301024235023/69046 596451111593912163277961/28248 -5609403368997817686249127547/4641050 495057205241079648212477525/6652 -801165718135489957347924991853/159054 29149963634884862421418123812691/79856 -2479392929313226753685415739663229/87058 84483613348880041862046775994036021/35460 -1215233140483755572040304994079820246041491/56786730
library(gmp)indices<-c(0,1,2*(1:30))bnums<-sapply(indices,BernoulliQ)names(bnums)<-paste0("B(",indices,")")print(bnums,initLine=FALSE)
$`B(0)`[1] 1$`B(1)`[1] 1/2$`B(2)`[1] 1/6$`B(4)`[1] -1/30$`B(6)`[1] 1/42$`B(8)`[1] -1/30$`B(10)`[1] 5/66$`B(12)`[1] -691/2730$`B(14)`[1] 7/6$`B(16)`[1] -3617/510$`B(18)`[1] 43867/798$`B(20)`[1] -174611/330$`B(22)`[1] 854513/138$`B(24)`[1] -236364091/2730$`B(26)`[1] 8553103/6$`B(28)`[1] -23749461029/870$`B(30)`[1] 8615841276005/14322$`B(32)`[1] -7709321041217/510$`B(34)`[1] 2577687858367/6$`B(36)`[1] -26315271553053477373/1919190$`B(38)`[1] 2929993913841559/6$`B(40)`[1] -261082718496449122051/13530$`B(42)`[1] 1520097643918070802691/1806$`B(44)`[1] -27833269579301024235023/690$`B(46)`[1] 596451111593912163277961/282$`B(48)`[1] -5609403368997817686249127547/46410$`B(50)`[1] 495057205241079648212477525/66$`B(52)`[1] -801165718135489957347924991853/1590$`B(54)`[1] 29149963634884862421418123812691/798$`B(56)`[1] -2479392929313226753685415739663229/870$`B(58)`[1] 84483613348880041862046775994036021/354$`B(60)`[1] -1215233140483755572040304994079820246041491/56786730
This implements, firstly, the algorithm specified with the task... then the better performingbernoulli.3,which uses the "double sum formula" listed under REXX. The number generators all (there is also abernoulli.2)use the same emmitter... it's just a matter of how long to wait for the emission.
#lang racket;; For: http://rosettacode.org/wiki/Bernoulli_numbers;; As described in task...(define (bernoulli.1 n) (define A (make-vector (add1 n))) (for ((m (in-range 0 (add1 n)))) (vector-set! A m (/ (add1 m))) (for ((j (in-range m (sub1 1) -1))) (define new-A_j-1 (* j (- (vector-ref A (sub1 j)) (vector-ref A j)))) (vector-set! A (sub1 j) new-A_j-1))) (vector-ref A 0))(define (non-zero-bernoulli-indices s) (sequence-filter (λ (n) (or (even? n) (= n 1))) s))(define (bernoulli_0..n B N) (for/list ((n (non-zero-bernoulli-indices (in-range (add1 N))))) (B n)));; From REXX description / http://mathworld.wolfram.com/BernoulliNumber.html #33;; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~;; bernoulli.2 is for illustrative purposes, binomial is very costly if there is no memoisation;; (which math/number-theory doesn't do)(require (only-in math/number-theory binomial))(define (bernoulli.2 n) (for/sum ((k (in-range 0 (add1 n)))) (* (/ (add1 k)) (for/sum ((r (in-range 0 (add1 k)))) (* (expt -1 r) (binomial k r) (expt r n))))));; Three things to do:;; 1. (expt -1 r): is 1 for even r, -1 for odd r... split the sum between those two.;; 2. splitting the sum might has arithmetic advantages, too. We're using rationals, so the smaller;; summations should require less normalisation of intermediate, fractional results;; 3. a memoised binomial... although the one from math/number-theory is fast, it is (and its;; factorials are) computed every time which is redundant(define kCr-memo (make-hasheq))(define !-memo (make-vector 1000 #f))(vector-set! !-memo 0 1) ;; seed the memo(define (! k) (cond [(vector-ref !-memo k) => values] [else (define k! (* k (! (- k 1)))) (vector-set! !-memo k k!) k!]))(define (kCr k r) ; If we want (kCr ... r>1000000) we'll have to reconsider this. However, until then... (define hash-key (+ (* 1000000 k) r)) (hash-ref! kCr-memo hash-key (λ () (/ (! k) (! r) (! (- k r))))))(define (bernoulli.3 n) (for/sum ((k (in-range 0 (add1 n)))) (define k+1 (add1 k)) (* (/ k+1) (- (for/sum ((r (in-range 0 k+1 2))) (* (kCr k r) (expt r n))) (for/sum ((r (in-range 1 k+1 2))) (* (kCr k r) (expt r n)))))))(define (display/align-fractions caption/idx-fmt Bs) ;; widths are one more than the order of magnitude (define oom+1 (compose add1 order-of-magnitude)) (define-values (I-width N-width D-width) (for/fold ((I 0) (N 0) (D 0)) ((b Bs) (n (non-zero-bernoulli-indices (in-naturals)))) (define +b (abs b)) (values (max I (oom+1 (max n 1))) (max N (+ (oom+1 (numerator +b)) (if (negative? b) 1 0))) (max D (oom+1 (denominator +b)))))) (define (~a/w/a n w a) (~a n #:width w #:align a)) (for ((n (non-zero-bernoulli-indices (in-naturals))) (b Bs)) (printf "~a ~a/~a~%" (format caption/idx-fmt (~a/w/a n I-width 'right)) (~a/w/a (numerator b) N-width 'right) (~a/w/a (denominator b) D-width 'left))))(module+ main (display/align-fractions "B(~a) =" (bernoulli_0..n bernoulli.3 60)))(module+ test (require rackunit) ; correctness and timing tests (check-match (time (bernoulli_0..n bernoulli.1 60)) (list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...)) (check-match (time (bernoulli_0..n bernoulli.2 60)) (list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...)) (check-match (time (bernoulli_0..n bernoulli.3 60)) (list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...)) ; timing only ... (void (time (bernoulli_0..n bernoulli.3 100))))
B( 0) = 1/1 B( 1) = -1/2 B( 2) = 1/6 B( 4) = -1/30 B( 6) = 1/42 B( 8) = -1/30 B(10) = 5/66 B(12) = -691/2730 B(14) = 7/6 B(16) = -3617/510 B(18) = 43867/798 B(20) = -174611/330 B(22) = 854513/138 B(24) = -236364091/2730 B(26) = 8553103/6 B(28) = -23749461029/870 B(30) = 8615841276005/14322 B(32) = -7709321041217/510 B(34) = 2577687858367/6 B(36) = -26315271553053477373/1919190 B(38) = 2929993913841559/6 B(40) = -261082718496449122051/13530 B(42) = 1520097643918070802691/1806 B(44) = -27833269579301024235023/690 B(46) = 596451111593912163277961/282 B(48) = -5609403368997817686249127547/46410 B(50) = 495057205241079648212477525/66 B(52) = -801165718135489957347924991853/1590 B(54) = 29149963634884862421418123812691/798 B(56) = -2479392929313226753685415739663229/870 B(58) = 84483613348880041862046775994036021/354 B(60) = -1215233140483755572040304994079820246041491/56786730
(formerly Perl 6)
First, a straighforward implementation of the naïve algorithm in the task description.
subbernoulli($n) {my@a;for0..$n ->$m {@a[$m] =FatRat.new(1,$m +1);forreverse1..$m ->$j {@a[$j -1] =$j * (@a[$j -1] -@a[$j]); } }return@a[0];}constant@bpairs =grep *.value.so, ($_ =>bernoulli($_)for0..60);my$width =max@bpairs.map: *.value.numerator.chars;my$form ="B(%2d) = \%{$width}d/%d\n";printf$form, .key, .value.nudefor@bpairs;
B( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
Here is a much faster way, following the Perl solution that avoids recalculating previous values each time through the function. We do this in Raku by not defining it as a function at all, but by defining it as an infinite sequence that we can read however many values we like from (52, in this case, to get up to B(100)). In this solution we've also avoided subscripting operations; rather we use a sequence operator (...) iterated over the list of the previous solution to find the next solution. We reverse the array in this case to make reference to the previous value in the list more natural, which means we take the last value of the list rather than the first value, and do so conditionally to avoid 0 values.
constantbernoulli =gather {my@a;for0..* ->$m {@a =FatRat.new(1,$m +1), ->$prev {my$j =@a.elems;$j * (@a.shift -$prev); } ... {not@a.elems }take$m =>@a[*-1]if@a[*-1]; }}constant@bpairs =bernoulli[^52];my$width =max@bpairs.map: *.value.numerator.chars;my$form ="B(%d)\t= \%{$width}d/%d\n";printf$form, .key, .value.nudefor@bpairs;
B(0)= 1/1B(1)= 1/2B(2)= 1/6B(4)= -1/30B(6)= 1/42B(8)= -1/30B(10)= 5/66B(12)= -691/2730B(14)= 7/6B(16)= -3617/510B(18)= 43867/798B(20)= -174611/330B(22)= 854513/138B(24)= -236364091/2730B(26)= 8553103/6B(28)= -23749461029/870B(30)= 8615841276005/14322B(32)= -7709321041217/510B(34)= 2577687858367/6B(36)= -26315271553053477373/1919190B(38)= 2929993913841559/6B(40)= -261082718496449122051/13530B(42)= 1520097643918070802691/1806B(44)= -27833269579301024235023/690B(46)= 596451111593912163277961/282B(48)= -5609403368997817686249127547/46410B(50)= 495057205241079648212477525/66B(52)= -801165718135489957347924991853/1590B(54)= 29149963634884862421418123812691/798B(56)= -2479392929313226753685415739663229/870B(58)= 84483613348880041862046775994036021/354B(60)= -1215233140483755572040304994079820246041491/56786730B(62)= 12300585434086858541953039857403386151/6B(64)= -106783830147866529886385444979142647942017/510B(66)= 1472600022126335654051619428551932342241899101/64722B(68)= -78773130858718728141909149208474606244347001/30B(70)= 1505381347333367003803076567377857208511438160235/4686B(72)= -5827954961669944110438277244641067365282488301844260429/140100870B(74)= 34152417289221168014330073731472635186688307783087/6B(76)= -24655088825935372707687196040585199904365267828865801/30B(78)= 414846365575400828295179035549542073492199375372400483487/3318B(80)= -4603784299479457646935574969019046849794257872751288919656867/230010B(82)= 1677014149185145836823154509786269900207736027570253414881613/498B(84)= -2024576195935290360231131160111731009989917391198090877281083932477/3404310B(86)= 660714619417678653573847847426261496277830686653388931761996983/6B(88)= -1311426488674017507995511424019311843345750275572028644296919890574047/61410B(90)= 1179057279021082799884123351249215083775254949669647116231545215727922535/272118B(92)= -1295585948207537527989427828538576749659341483719435143023316326829946247/1410B(94)= 1220813806579744469607301679413201203958508415202696621436215105284649447/6B(96)= -211600449597266513097597728109824233673043954389060234150638733420050668349987259/4501770B(98)= 67908260672905495624051117546403605607342195728504487509073961249992947058239/6B(100)= -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330
And if you're a pure enough FP programmer to dislike destroying and reconstructing the array each time, here's the same algorithm without side effects. We use zip with the pair constructor=> to keep values associated with their indices. This provides sufficient local information that we can define our own binary operator "bop" to reduce between each two terms, using the "triangle" form (called "scan" in Haskell) to return the intermediate results that will be important to compute the next Bernoulli number.
subinfix:<bop>(\prev, \this) {this.key =>this.key * (this.value -prev.value)}subnext-bernoulli ( (:key($pm), :value(@pa)) ) {$pm +1 => [map *.value, [\bop] ($pm +2 ...1)Z=>FatRat.new(1,$pm +2), |@pa ]}constantbernoulli =grep *.value,map { .key => .value[*-1] }, (0 => [FatRat.new(1,1)],&next-bernoulli ... *);constant@bpairs =bernoulli[^52];my$width =max@bpairs.map: *.value.numerator.chars;my$form ="B(%d)\t= \%{$width}d/%d\n";printf$form, .key, .value.nudefor@bpairs;
Same output as memoization example
IncludeHow to use
IncludeSource code
We may choose from several algorithms: Double sum formula, Generating function and Akiyama-Tanigawa.
Of these, the summation derived from the generating function looks most simple. Any of these may be implemented using Fractional or Floating calculations. If Bernoulli numbers are used in series expansions or formulas, the numerator/denominator representation is not very useful.
In below program, the first 2 routines use fractions (integer arithmetic), as per the task requirements. The last 2 routines perform decimal calculations, and do not try to find the numerators and denominators.
Rational1 uses functions from module Rational, making it a quite straightforward translation of the generating function summation. Rational2, using the double summation formula, is highly optimized. Decimal1 uses again the generating function formula, and Decimal2 implements Akiyama-Tanigawa.
--23Aug2025includeSettingnumericdigits16argxxifxx=''thenxx=60ss=(xx>0);xx=Abs(xx)say'BERNOULLI NUMBERS'sayversionsaycallRational1xxifssthencallShowFractionsxxcallTimer'R'callRational2xxifssthencallShowFractionsxxcallTimer'R'callDecimal1xxifssthencallShowDecimalsxxcallTimer'R'callDecimal2xxifssthencallShowDecimalsxxcallTimer'R'exitRational1:procedureexposeMemo.Bern.argxxsay'Rational arithmetic'say'cf Generating function'numericdigitsMax(Digits(),3*xx)Bern.=0;Bern.0='1 1';Bern.1='-1 2'don=2by2toxxs=0dok=0ton-1s=Addq(s,Scaleq(Bern.k,Comb(n+1,k)))endkBern.n=Mulq(Negq(s),Invq(n+1))endnreturnRational2:procedureexposeMemo.Bern.argxxsay'Numerator and denominator calculations'say'cf Double sum formula optimized'numericdigitsMax(Digits(),3*xx)Bern.=0;Bern.0='1 1';Bern.1='-1 2'doj=2by2toxxjp=j+1;sn=1-j;sd=2dok=2by2toj-1bn=bn.k;ad=ad.k;an=Comb(jp,k)*bn;tl=Lcm(sd,ad)sn=tl%sd*sn;sd=tl;an=tl%ad*an;sn=sn+anendksn=-sn;sd=sd*jp;bn.j=sn;ad.j=sd;g=Gcd(Abs(sn),sd)Bern.j=sn/gsd/gendjreturnShowFractions:procedureexposeBern.argxxw=Length(Word(Bern.xx,1))+2say' Bn'Right('Num',w)'/'Left('Den',10)'Decimal'doi=0toxxifBern.i<>0thendoparsevarBern.inumdensayRight(i,3)Right(num,w)'/'Left(den,10)Digit(num/den,16)endendireturnDecimal1:procedureexposeBern.Memo.argxxsay'Floating point calculations'say'cf Generating function'numericdigitsMax(Digits(),2*xx)Bern.=0;Bern.0=1;Bern.1='-0.5'don=2by2toxxs=0dok=0ton-1s=s+Bern.k*Comb(n+1,k)endkBern.n=-s/(n+1)endnreturnxxDecimal2:procedureexposeBern.Memo.argxxsay'Floating point calculations'say'cf Akiyama-Tanigawa'numericdigitsMax(Digits(),3*xx)Bern.=0;Bern.0=1;Bern.1=-0.5doi=2by2toxxdom=0toia.m=1/(m+1)doj=mby-1to1j1=j-1;a.j1=j*(a.j1-a.j)endendmBern.i=a.0endireturnxxShowDecimals:procedureexposeBern.argxxsay' Bn''Decimal'doi=0toxxifBern.i<>0thensayRight(i,3)Digit(Bern.i,16)endireturnincludeMath
Running without parameters (thus up to B60).
BERNOULLI NUMBERSREXX-Regina_3.9.7(MT) 5.00 18 Mar 2025Rational arithmeticcf Generating function Bn Num / Den Decimal 0 1 / 1 1 1 -1 / 2 -0.5 2 1 / 6 0.1666666666666667 4 -1 / 30 -0.03333333333333333 6 1 / 42 0.02380952380952381 8 -1 / 30 -0.03333333333333333 10 5 / 66 0.07575757575757576 12 -691 / 2730 -0.2531135531135531 14 7 / 6 1.166666666666667 16 -3617 / 510 -7.092156862745098 18 43867 / 798 54.97117794486216 20 -174611 / 330 -529.1242424242424 22 854513 / 138 6192.123188405797 24 -236364091 / 2730 -86580.25311355311 26 8553103 / 6 1425517.166666667 28 -23749461029 / 870 -27298231.06781609 30 8615841276005 / 14322 601580873.9006424 32 -7709321041217 / 510 -15116315767.09216 34 2577687858367 / 6 429614643061.1667 36 -26315271553053477373 / 1919190 -13711655205088.33 38 2929993913841559 / 6 488332318973593.2 40 -261082718496449122051 / 13530 -1.929657934194007E+16 42 1520097643918070802691 / 1806 8.416930475736826E+17 44 -27833269579301024235023 / 690 -4.033807185405946E+19 46 596451111593912163277961 / 282 2.115074863808199E+21 48 -5609403368997817686249127547 / 46410 -1.208662652229653E+23 50 495057205241079648212477525 / 66 7.500866746076964E+24 52 -801165718135489957347924991853 / 1590 -5.038778101481069E+26 54 29149963634884862421418123812691 / 798 3.652877648481812E+28 56 -2479392929313226753685415739663229 / 870 -2.849876930245088E+30 58 84483613348880041862046775994036021 / 354 2.386542749968363E+32 60 -1215233140483755572040304994079820246041491 / 56786730 -2.139994925722533E+340.026 secondsNumerator and denominator calculationscf Double sum formula optimized Bn Num / Den Decimal 0 1 / 1 1 1 -1 / 2 -0.5 2 1 / 6 0.1666666666666667 4 -1 / 30 -0.03333333333333333 6 1 / 42 0.02380952380952381 8 -1 / 30 -0.03333333333333333 10 5 / 66 0.07575757575757576 12 -691 / 2730 -0.2531135531135531 14 7 / 6 1.166666666666667 16 -3617 / 510 -7.092156862745098 18 43867 / 798 54.97117794486216 20 -174611 / 330 -529.1242424242424 22 854513 / 138 6192.123188405797 24 -236364091 / 2730 -86580.25311355311 26 8553103 / 6 1425517.166666667 28 -23749461029 / 870 -27298231.06781609 30 8615841276005 / 14322 601580873.9006424 32 -7709321041217 / 510 -15116315767.09216 34 2577687858367 / 6 429614643061.1667 36 -26315271553053477373 / 1919190 -13711655205088.33 38 2929993913841559 / 6 488332318973593.2 40 -261082718496449122051 / 13530 -1.929657934194007E+16 42 1520097643918070802691 / 1806 8.416930475736826E+17 44 -27833269579301024235023 / 690 -4.033807185405946E+19 46 596451111593912163277961 / 282 2.115074863808199E+21 48 -5609403368997817686249127547 / 46410 -1.208662652229653E+23 50 495057205241079648212477525 / 66 7.500866746076964E+24 52 -801165718135489957347924991853 / 1590 -5.038778101481069E+26 54 29149963634884862421418123812691 / 798 3.652877648481812E+28 56 -2479392929313226753685415739663229 / 870 -2.849876930245088E+30 58 84483613348880041862046775994036021 / 354 2.386542749968363E+32 60 -1215233140483755572040304994079820246041491 / 56786730 -2.139994925722533E+340.005 secondsFloating point calculationscf Generating function Bn Decimal 0 1 1 -0.5 2 0.1666666666666667 4 -0.03333333333333333 6 0.02380952380952381 8 -0.03333333333333333 10 0.07575757575757576 12 -0.2531135531135531 14 1.166666666666667 16 -7.092156862745098 18 54.97117794486216 20 -529.1242424242424 22 6192.123188405797 24 -86580.25311355311 26 1425517.166666667 28 -27298231.06781609 30 601580873.9006424 32 -15116315767.09216 34 429614643061.1667 36 -13711655205088.33 38 488332318973593.2 40 -1.929657934194007E+16 42 8.416930475736826E+17 44 -4.033807185405946E+19 46 2.115074863808199E+21 48 -1.208662652229653E+23 50 7.500866746076964E+24 52 -5.038778101481069E+26 54 3.652877648481812E+28 56 -2.849876930245088E+30 58 2.386542749968363E+32 60 -2.139994925722533E+340.005 secondsFloating point calculationscf Akiyama-Tanigawa Bn Decimal 0 1 1 -0.5 2 0.1666666666666667 4 -0.03333333333333333 6 0.02380952380952381 8 -0.03333333333333333 10 0.07575757575757576 12 -0.2531135531135531 14 1.166666666666667 16 -7.092156862745098 18 54.97117794486216 20 -529.1242424242424 22 6192.123188405797 24 -86580.25311355311 26 1425517.166666667 28 -27298231.06781609 30 601580873.9006424 32 -15116315767.09216 34 429614643061.1667 36 -13711655205088.33 38 488332318973593.2 40 -1.929657934194007E+16 42 8.416930475736826E+17 44 -4.033807185405946E+19 46 2.115074863808199E+21 48 -1.208662652229653E+23 50 7.500866746076964E+24 52 -5.038778101481069E+26 54 3.652877648481812E+28 56 -2.849876930245088E+30 58 2.386542749968363E+32 60 -2.139994925722533E+340.048 seconds
This is very fast. But which algorithm is fastest? Let's run with parameter -500 (up to B500 having about 750 decimal digits, but only timings shown).
BERNOULLI NUMBERSREXX-ooRexx_5.1.0(MT)_64-bit 6.05 2 May 2025Rational arithmeticcf Generating function35.874 secondsNumerator and denominator calculationscf Double sum formula optimized23.285 secondsFloating point calculationscf Generating function24.895 secondsFloating point calculationscf Akiyama-Tanigawa165.937 seconds
Fractions such as a/b are here handled through the complex number data structure(a,b).2 local words support the algorithm suggested by the task:FracSub substracts 2 fractions andFracSimp make a fraction irreducibleUnfortunately, floating-point precision prevents from going beyond B22.
| RPL code | Comment |
|---|---|
≪ DUP2 IM * ROT IM ROT RE * -≫ 'FracSub' STO≪ DUP C→R ABS SWAP ABS DUP2 < ≪ SWAP ≫ IFTWHILE DUPREPEAT SWAP OVER MODEND DROP /≫ 'FracSimp' STO≪ { } 1 ROT 1 +FOR m 1 m R→C +IF m 2 ≥THEN m 2FOR j DUP j 1 - GET OVER j GETFracSub C→R SWAP j 1 - * SWAP R→CFracSimp j 1 - SWAP PUT -1STEP ENDNEXT 1 GET DUP RE SWAP 0IFTE≫ 'BRNOU' STO | ( (a,b) (c,d) -- (e,f) ) with e/f = a/b - c/d( (a,b) -- (c,d) ) with c/d simplified fraction of a/bget GCD of a and bdivide (a,b) by GCD( n -- (a,b) ) with a/b = B(n)For m from 1 to n+1 do A[m] ← 1/m for j from m to 2 do (A[j-1] - A[j]) . * (j-1) A[j-1] ←return A[1] as a fraction or zero |
5BRNOU22BRNOU
2: 01: (854513,138)
Latest RPL implementations can natively handle long fractions and generate Bernoulli numbers.
≪ { } 0 ROTFOR nIF n 2 > LASTARG MOD AND NOTTHEN n IBERNOULLI +ENDNEXT≫ 'TASK' STO60TASK1: {1 -1/2 1/6 -1/30 1/42 -1/30 5/66 -691/2730 7/6 -3617/510 43867/798 -174611/330 854513/138 -236364091/2730 8553103/6 -23749461029/870 8615841276005/14322 -7709321041217/510 2577687858367/6 -26315271553053477373/1919190 2929993913841559/6 -261082718496449122051/13530 1520097643918070802691/1806 -27833269579301024235023/690 596451111593912163277961/282 -5609403368997817686249127547/46410 495057205241079648212477525/66 -801165718135489957347924991853/1590 9149963634884862421418123812691/798 -2479392929313226753685415739663229/870 84483613348880041862046775994036021/354 -1215233140483755572040304994079820246041491/56786730}Runs in 3 minutes 40 on a HP-50g, against 1 hour and 30 minutes if calculating Bernoulli numbers with the above function.
bernoulli=Enumerator.newdo|y|ar=[]0.stepdo|m|ar<<Rational(1,m+1)m.downto(1){|j|ar[j-1]=j*(ar[j-1]-ar[j])}y<<ar.first# yieldendendb_nums=bernoulli.take(61)width=b_nums.map{|b|b.numerator.to_s.size}.maxb_nums.each_with_index{|b,i|puts"B(%2i) = %*i/%i"%[i,width,b.numerator,b.denominator]unlessb.zero?}
B( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
// 2.5 implementations presented here: naive, optimized, and an iterator using// the optimized function. The speeds vary significantly: relative// speeds of optimized:iterator:naive implementations is 625:25:1.#![feature(test)]externcratenum;externcratetest;usenum::bigint::{BigInt,ToBigInt};usenum::rational::{BigRational};usestd::cmp::max;usestd::env;usestd::ops::{Mul,Sub};usestd::process;structBn{value:BigRational,index:i32}structContext{bigone_const:BigInt,a:Vec<BigRational>,index:i32// Counter for iterator implementation}implContext{pubfnnew()->Context{letbigone=1.to_bigint().unwrap();leta_vec:Vec<BigRational>=vec![];Context{bigone_const:bigone,a:a_vec,index:-1}}}implIteratorforContext{typeItem=Bn;fnnext(&mutself)->Option<Bn>{self.index+=1;Some(Bn{value:bernoulli(self.indexasusize,self),index:self.index})}}fnhelp(){println!("Usage: bernoulli_numbers <up_to>");}fnmain(){letargs:Vec<String>=env::args().collect();letmutup_to:usize=60;matchargs.len(){1=>{},2=>{up_to=args[1].parse::<usize>().unwrap();},_=>{help();process::exit(0);}}letcontext=Context::new();// Collect the solutions by using the Context iterator// (this is not as fast as calling the optimized function directly).letres=context.take(up_to+1).collect::<Vec<_>>();letwidth=res.iter().fold(0,|a,r|max(a,r.value.numer().to_string().len()));forrinres.iter().filter(|r|*r.value.numer()!=ToBigInt::to_bigint(&0).unwrap()){println!("B({:>2}) = {:>2$} / {denom}",r.index,r.value.numer(),width,denom=r.value.denom());}}// Implementation with no reused calculations.fn_bernoulli_naive(n:usize,c:&mutContext)->BigRational{formin0..n+1{c.a.push(BigRational::new(c.bigone_const.clone(),(m+1).to_bigint().unwrap()));forjin(1..m+1).rev(){c.a[j-1]=(c.a[j-1].clone().sub(c.a[j].clone())).mul(BigRational::new(j.to_bigint().unwrap(),c.bigone_const.clone()));}}c.a[0].reduced()}// Implementation with reused calculations (does not require sequential calls).fnbernoulli(n:usize,c:&mutContext)->BigRational{foriin0..n+1{ifi>=c.a.len(){c.a.push(BigRational::new(c.bigone_const.clone(),(i+1).to_bigint().unwrap()));forjin(1..i+1).rev(){c.a[j-1]=(c.a[j-1].clone().sub(c.a[j].clone())).mul(BigRational::new(j.to_bigint().unwrap(),c.bigone_const.clone()));}}}c.a[0].reduced()}#[cfg(test)]modtests{usesuper::{Bn,Context,bernoulli,_bernoulli_naive};usenum::rational::{BigRational};usestd::str::FromStr;usetest::Bencher;// [tests elided]#[bench]fnbench_bernoulli_naive(b:&mutBencher){letmutcontext=Context::new();b.iter(||{letmutres:Vec<Bn>=vec![];fornin0..30+1{letb=_bernoulli_naive(n,&mutcontext);res.push(Bn{value:b.clone(),index:nasi32});}});}#[bench]fnbench_bernoulli(b:&mutBencher){letmutcontext=Context::new();b.iter(||{letmutres:Vec<Bn>=vec![];fornin0..30+1{letb=bernoulli(n,&mutcontext);res.push(Bn{value:b.clone(),index:nasi32});}});}#[bench]fnbench_bernoulli_iter(b:&mutBencher){b.iter(||{letcontext=Context::new();let_res=context.take(30+1).collect::<Vec<_>>();});}}
B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
With Custom Rational Number Class
(code will run in Scala REPL with a cut-and-paste without need for a third-party library)
/** Roll our own pared-down BigFraction class just for these Bernoulli Numbers */caseclassBFraction(numerator:BigInt,denominator:BigInt){require(denominator!=BigInt(0),"Denominator cannot be zero")valgcd=numerator.gcd(denominator)valnum=numerator/gcdvalden=denominator/gcddefunary_-=BFraction(-num,den)def-(that:BFraction)=thatmatch{casefiff.num==BigInt(0)=>thiscasefiff.den==this.den=>BFraction(this.num-f.num,this.den)casef=>BFraction(((this.num*f.den)-(f.num*this.den)),this.den*f.den)}def*(that:Int)=BFraction(num*that,den)overridedeftoString=num+" / "+den}defbernoulliB(n:Int):BFraction={valaa:Array[BFraction]=Array.ofDim(n+1)for(m<-0ton){aa(m)=BFraction(1,(m+1))for(n<-mto1by-1){aa(n-1)=(aa(n-1)-aa(n))*n}}aa(0)}assert({valb12=bernoulliB(12);b12.num==-691&&b12.den==2730})valr=for(n<-0to60;b=bernoulliB(n)ifb.num!=0)yield(n,b)valnumeratorSize=r.map(_._2.num.toString.length).max// Print the resultsrforeach{case(i,b)=>{vallabel=f"b($i)"valnum=(" "*(numeratorSize-b.num.toString.length))+b.numprintln(f"$label%-6s$num /${b.den}")}}
b(0) 1 / 1b(1) 1 / 2b(2) 1 / 6b(4) -1 / 30b(6) 1 / 42b(8) -1 / 30b(10) 5 / 66b(12) -691 / 2730b(14) 7 / 6b(16) -3617 / 510b(18) 43867 / 798b(20) -174611 / 330b(22) 854513 / 138b(24) -236364091 / 2730b(26) 8553103 / 6b(28) -23749461029 / 870b(30) 8615841276005 / 14322b(32) -7709321041217 / 510b(34) 2577687858367 / 6b(36) -26315271553053477373 / 1919190b(38) 2929993913841559 / 6b(40) -261082718496449122051 / 13530b(42) 1520097643918070802691 / 1806b(44) -27833269579301024235023 / 690b(46) 596451111593912163277961 / 282b(48) -5609403368997817686249127547 / 46410b(50) 495057205241079648212477525 / 66b(52) -801165718135489957347924991853 / 1590b(54) 29149963634884862421418123812691 / 798b(56) -2479392929313226753685415739663229 / 870b(58) 84483613348880041862046775994036021 / 354b(60) -1215233140483755572040304994079820246041491 / 56786730
; Return the n'th Bernoulli number.(definebernoulli(lambda(n)(let((a(make-vector(1+n))))(do((m0(1+m)))((>mn))(vector-set!am(/1(1+m)))(do((jm(1-j)))((<j1))(vector-set!a(1-j)(*j(-(vector-refa(1-j))(vector-refaj))))))(vector-refa0)))); Convert a rational to a string. If an integer, ends with "/1".(definerational->string(lambda(rational)(format"~a/~a"(numeratorrational)(denominatorrational)))); Returns the string length of the numerator of a rational.(definerational-numerator-length(lambda(rational)(string-length(format"~a"(numeratorrational))))); Formats a rational with left-padding such that total length to the slash is as given.(definerational-padded(lambda(rationaltotal-length-to-slash)(let*((length-padding(-total-length-to-slash(rational-numerator-lengthrational)))(padding-string(make-stringlength-padding#\)))(string-appendpadding-string(rational->stringrational))))); Return the Bernoulli numbers 0 through n in a list.(definemake-bernoulli-list(lambda(n)(if(=n0)(list(bernoullin))(append(make-bernoulli-list(1-n))(list(bernoullin)))))); Print the non-zero Bernoulli numbers 0 through 60 aligning the slashes.(let*((bernoullis-list(make-bernoulli-list60))(numerator-lengths(maprational-numerator-lengthbernoullis-list))(max-numerator-length(applymaxnumerator-lengths)))(letprint-bernoulli((index0)(numbersbernoullis-list))(cond((null?numbers))((=0(carnumbers))(print-bernoulli(1+index)(cdrnumbers)))(else(printf"B(~2@a) = ~a~%"index(rational-padded(carnumbers)max-numerator-length))(print-bernoulli(1+index)(cdrnumbers))))))
$ scheme --script bernoulli.scmB( 0) = 1/1B( 1) = 1/2B( 2) = 1/6B( 4) = -1/30B( 6) = 1/42B( 8) = -1/30B(10) = 5/66B(12) = -691/2730B(14) = 7/6B(16) = -3617/510B(18) = 43867/798B(20) = -174611/330B(22) = 854513/138B(24) = -236364091/2730B(26) = 8553103/6B(28) = -23749461029/870B(30) = 8615841276005/14322B(32) = -7709321041217/510B(34) = 2577687858367/6B(36) = -26315271553053477373/1919190B(38) = 2929993913841559/6B(40) = -261082718496449122051/13530B(42) = 1520097643918070802691/1806B(44) = -27833269579301024235023/690B(46) = 596451111593912163277961/282B(48) = -5609403368997817686249127547/46410B(50) = 495057205241079648212477525/66B(52) = -801165718135489957347924991853/1590B(54) = 29149963634884862421418123812691/798B(56) = -2479392929313226753685415739663229/870B(58) = 84483613348880041862046775994036021/354B(60) = -1215233140483755572040304994079820246041491/56786730
The program below usesbigRationalnumbers. The Bernoulli numbers are written as fraction and as decimal number, with possible repeating decimals.The conversion of a bigRational number tostring is donewith the functionstr. Thisfunction automatically writes repeating decimals in parentheses, when necessary.
$ include "seed7_05.s7i"; include "bigrat.s7i";const func bigRational: bernoulli (in integer: n) is func result var bigRational: bernoulli is bigRational.value; local var integer: m is 0; var integer: j is 0; var array bigRational: a is 0 times bigRational.value; begin a := [0 .. n] times bigRational.value; for m range 0 to n do a[m] := 1_ / bigInteger(succ(m)); for j range m downto 1 do a[pred(j)] := bigRational(j) * (a[j] - a[pred(j)]); end for; end for; bernoulli := a[0]; end func;const proc: main is func local var bigRational: bernoulli is bigRational.value; var integer: i is 0; begin for i range 0 to 60 do bernoulli := bernoulli(i); if bernoulli <> bigRational.value then writeln("B(" <& i lpad 2 <& ") = " <& bernoulli.numerator lpad 44 <& " / " <& bernoulli.denominator rpad 8 <& " " <& bernoulli); end if; end for; end func;B( 0) = 1 / 1 1.0B( 1) = -1 / 2 -0.5B( 2) = 1 / 6 0.1(6)B( 4) = -1 / 30 -0.0(3)B( 6) = 1 / 42 0.0(238095)B( 8) = -1 / 30 -0.0(3)B(10) = 5 / 66 0.0(75)B(12) = -691 / 2730 -0.2(531135)B(14) = 7 / 6 1.1(6)B(16) = -3617 / 510 -7.0(9215686274509803)B(18) = 43867 / 798 54.9(711779448621553884)B(20) = -174611 / 330 -529.1(24)B(22) = 854513 / 138 6192.1(2318840579710144927536)B(24) = -236364091 / 2730 -86580.2(531135)B(26) = 8553103 / 6 1425517.1(6)B(28) = -23749461029 / 870 -27298231.0(6781609195402298850574712643)B(30) = 8615841276005 / 14322 601580873.9(006423683843038681748359167714)B(32) = -7709321041217 / 510 -15116315767.0(9215686274509803)B(34) = 2577687858367 / 6 429614643061.1(6)B(36) = -26315271553053477373 / 1919190 -13711655205088.3(327721590879485616)B(38) = 2929993913841559 / 6 488332318973593.1(6)B(40) = -261082718496449122051 / 13530 -19296579341940068.1(4863266814)B(42) = 1520097643918070802691 / 1806 841693047573682615.0(005537098560354374307862679955703211517165)B(44) = -27833269579301024235023 / 690 -40338071854059455413.0(7681159420289855072463)B(46) = 596451111593912163277961 / 282 2115074863808199160560.1(4539007092198581560283687943262411347517730496)B(48) = -5609403368997817686249127547 / 46410 -120866265222965259346027.3(119370825253178194354664942900237017884076707606)B(50) = 495057205241079648212477525 / 66 7500866746076964366855720.0(75)B(52) = -801165718135489957347924991853 / 1590 -503877810148106891413789303.0(5220125786163)B(54) = 29149963634884862421418123812691 / 798 36528776484818123335110430842.9(711779448621553884)B(56) = -2479392929313226753685415739663229 / 870 -2849876930245088222626914643291.0(6781609195402298850574712643)B(58) = 84483613348880041862046775994036021 / 354 238654274996836276446459819192192.1(4971751412429378531073446327683615819209039548022598870056)B(60) = -1215233140483755572040304994079820246041491 / 56786730 -21399949257225333665810744765191097.3(926741511617238745742183076926598872659158222352299560126106)
Built-in:
saybernoulli(42).as_frac#=> 1520097643918070802691/1806
Recursive solution (with auto-memoization):
funcbernoulli_number(n)iscached{n.is_one&&return1/2n.is_odd&&return01-sum(^n,{|k|binomial(n,k)*__FUNC__(k)/(n-k+1)})}fornin(0..60){varBn=bernoulli_number(n)||nextprintf("B(%2d) = %44s / %s\n",n,Bn.nude)}
Using Ramanujan's congruences (pretty fast):
funcramanujan_bernoulli_number(n)iscached{return1/2ifn.is_onereturn0ifn.is_odd((n%6==4?-1/2:1)*(n+3)/3-sum(1..(n-n%6)/6,{|k|binomial(n+3,n-6*k)*__FUNC__(n-6*k)}))/binomial(n+3,n)}
Using Euler's product formula for the Riemann zeta function and the Von Staudt–Clausen theorem (very fast):
funcbernoulli_number_from_zeta(n){n.is_zero&&return1n.is_one&&return1/2n.is_odd&&return0varlog2B=(log(4*Num.tau*n)/2+n*log(n)-n*log(Num.tau)-n)/log(2)localNum!PREC=*(int(n+log2B)+(n<=90?18:0))varK=2*(n!/ Num.tau**n) var d = n.divisors.grep {|k| is_prime(k+1) }.prod {|k| k+1 } var z = ceil((K*d).root(n-1)).primes.prod {|p| 1 - p.float**(-n) } (-1)**(n/2+1)*int(ceil(d*K/z))/d}
The Akiyama–Tanigawa algorithm:
funcbernoulli_print{vara=[]formin(0..60){a<<1/(m+1)forjin(1..m->flip){(a[j-1]-=a[j])*=j}a[0]||nextprintf("B(%2d) = %44s / %s\n",m,a[0].nude)}} bernoulli_print()
B( 0) = 1 / 1B( 1) = 1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
for n in 0..60 | (b:=bernoulli(n)$INTHEORY; b~=0) repeat print [n,b]
Package:IntegerNumberTheoryFunctions
===============Format: [n,B_n] =============== [0,1] 1 [1,- -] 2 1 [2,-] 6 1 [4,- --] 30 1 [6,--] 42 1 [8,- --] 30 5 [10,--] 66 691 [12,- ----] 2730 7 [14,-] 6 3617 [16,- ----] 510 43867 [18,-----] 798 174611 [20,- ------] 330 854513 [22,------] 138 236364091 [24,- ---------] 2730 8553103 [26,-------] 6 23749461029 [28,- -----------] 870 8615841276005 [30,-------------] 14322 7709321041217 [32,- -------------] 510 2577687858367 [34,-------------] 6 26315271553053477373 [36,- --------------------] 1919190 2929993913841559 [38,----------------] 6 261082718496449122051 [40,- ---------------------] 13530 1520097643918070802691 [42,----------------------] 1806 27833269579301024235023 [44,- -----------------------] 690 596451111593912163277961 [46,------------------------] 282 5609403368997817686249127547 [48,- ----------------------------] 46410 495057205241079648212477525 [50,---------------------------] 66 801165718135489957347924991853 [52,- ------------------------------] 1590 29149963634884862421418123812691 [54,--------------------------------] 798 2479392929313226753685415739663229 [56,- ----------------------------------] 870 84483613348880041862046775994036021 [58,-----------------------------------] 354 1215233140483755572040304994079820246041491 [60,- -------------------------------------------] 56786730 Type: Void
Uses the Frac type defined in theRational task.
importBigIntpublicfuncbernoulli<T:BinaryInteger&SignedNumeric>(n:Int)->Frac<T>{guardn!=0else{return1}vararr=[Frac<T>]()formin0...n{arr.append(Frac(numerator:1,denominator:T(m)+1))forjinstride(from:m,through:1,by:-1){arr[j-1]=(arr[j-1]-arr[j])*Frac(numerator:T(j),denominator:1)}}returnarr[0]}fornin0...60{letb=bernoulli(n:n)asFrac<BigInt>guardb!=0else{continue}print("B(\(n)) =\(b)")}
B(0) = Frac(1 / 1)B(1) = Frac(1 / 2)B(2) = Frac(1 / 6)B(4) = Frac(-1 / 30)B(6) = Frac(1 / 42)B(8) = Frac(-1 / 30)B(10) = Frac(5 / 66)B(12) = Frac(-691 / 2730)B(14) = Frac(7 / 6)B(16) = Frac(-3617 / 510)B(18) = Frac(43867 / 798)B(20) = Frac(-174611 / 330)B(22) = Frac(854513 / 138)B(24) = Frac(-236364091 / 2730)B(26) = Frac(8553103 / 6)B(28) = Frac(-23749461029 / 870)B(30) = Frac(8615841276005 / 14322)B(32) = Frac(-7709321041217 / 510)B(34) = Frac(2577687858367 / 6)B(36) = Frac(-26315271553053477373 / 1919190)B(38) = Frac(2929993913841559 / 6)B(40) = Frac(-261082718496449122051 / 13530)B(42) = Frac(1520097643918070802691 / 1806)B(44) = Frac(-27833269579301024235023 / 690)B(46) = Frac(596451111593912163277961 / 282)B(48) = Frac(-5609403368997817686249127547 / 46410)B(50) = Frac(495057205241079648212477525 / 66)B(52) = Frac(-801165718135489957347924991853 / 1590)B(54) = Frac(29149963634884862421418123812691 / 798)B(56) = Frac(-2479392929313226753685415739663229 / 870)B(58) = Frac(84483613348880041862046775994036021 / 354)B(60) = Frac(-1215233140483755572040304994079820246041491 / 56786730)
bernoulli (n): a =: new row for m =: from 0 upto na[m] =: 1 / (m + 1)for j =: from m downto 1 a[j-1] =: j * (a[j-1] - a[j]) return a[0]main (p):+ max =: string p.1 as integer else 60 \ max number as parameter for i =: from 2 upto max step 2b =: bernoulli i\ to align the slash, print numerator and denominator separatelyprint format "B(%2;) = %50; / %;" i, b.num, b.den
B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
Iterative version as used by Ada Lovelace in her Note G (old numbering)Seehttps://rclab.de/analyticalengine/inhalt
Significantly quicker in execution
main (parms):+ n =: string parms[1] as integer else 30 print n bernoullis\ ierative version as used by Ada Lovelace print (max) bernoullis: b =: [] b[1] =: 1/6 print 'B[1]=' _ b[1] ?# n =: from 2 upto maxk =: 2*n - 1n2 =: 2*na0 =: - (n2 - 1) / (n2 + 1) / 2a1 =: n2 / 2x =: a0 + (b[1] * a1)aj =: a1?# j =: from 3 upto k-1 step 2 aj =* (n2 - (j-2)) / j aj =* (n2 - (j-1)) / (j + 1) x =+ b[j] * ajb[k] =: -xprint format "B(%2;) = %50; / %;" k, b[k].num, b[k].den
Assigning common subexpression avoids very long and complex expressions that are difficult to read and verify. And its closer to the original.
Same as above version
procbernoulli{n}{for{setm0}{$m<=$n}{incrm}{lappendA[list1[expr{$m+1}]]for{setj$m}{[seti$j]>=1}{}{lassign[lindex$A[incrj-1]]a1b1lassign[lindex$A$i]a2b2setx[setp[expr{$i*($a1*$b2-$a2*$b1)}]]sety[setq[expr{$b1*$b2}]]while{$q}{setq[expr{$p%[setp$q]}]}lsetA$j[list[expr{$x/$p}][expr{$y/$p}]]}}return[lindex$A0]}setlen0for{setn0}{$n<=60}{incrn}{setb[bernoulli$n]if{[lindex$b0]}{lappendresult$n{*}$bsetlen[expr{max($len,[stringlength[lindex$b0]])}]}}foreach{nnumdenom}$result{puts[format{B_%-2d=%*lld/%lld}$n$len$num$denom]}
B_0 = 1/1B_1 = 1/2B_2 = 1/6B_4 = -1/30B_6 = 1/42B_8 = -1/30B_10 = 5/66B_12 = -691/2730B_14 = 7/6B_16 = -3617/510B_18 = 43867/798B_20 = -174611/330B_22 = 854513/138B_24 = -236364091/2730B_26 = 8553103/6B_28 = -23749461029/870B_30 = 8615841276005/14322B_32 = -7709321041217/510B_34 = 2577687858367/6B_36 = -26315271553053477373/1919190B_38 = 2929993913841559/6B_40 = -261082718496449122051/13530B_42 = 1520097643918070802691/1806B_44 = -27833269579301024235023/690B_46 = 596451111593912163277961/282B_48 = -5609403368997817686249127547/46410B_50 = 495057205241079648212477525/66B_52 = -801165718135489957347924991853/1590B_54 = 29149963634884862421418123812691/798B_56 = -2479392929313226753685415739663229/870B_58 = 84483613348880041862046775994036021/354B_60 = -1215233140483755572040304994079820246041491/56786730
' Bernoulli numbers - vb.net - 06/03/2017ImportsSystem.Numerics'BigIntegerModuleBernoulli_numbersFunctiongcd_BigInt(ByValxAsBigInteger,ByValyAsBigInteger)AsBigIntegerDimy2AsBigIntegerx=BigInteger.Abs(x)Doy2=BigInteger.Remainder(x,y)x=yy=y2LoopUntily=0ReturnxEndFunction'gcd_BigIntSubbernoul_BigInt(nAsInteger,ByRefbnumAsBigInteger,ByRefbdenAsBigInteger)Dimj,mAsIntegerDimfAsBigIntegerDimanum(),aden()AsBigIntegerReDimanum(n+1),aden(n+1)Form=0Tonanum(m+1)=1aden(m+1)=m+1Forj=mTo1Step-1anum(j)=j*(aden(j+1)*anum(j)-aden(j)*anum(j+1))aden(j)=aden(j)*aden(j+1)f=gcd_BigInt(BigInteger.Abs(anum(j)),BigInteger.Abs(aden(j)))Iff<>1Thenanum(j)=anum(j)/faden(j)=aden(j)/fEndIfNextNextbnum=anum(1):bden=aden(1)EndSub'bernoul_BigIntSubbernoulli_BigInt()DimiAsIntegerDimbnum,bdenAsBigIntegerbnum=0:bden=0Fori=0To60bernoul_BigInt(i,bnum,bden)Ifbnum<>0ThenConsole.WriteLine("B("&i&")="&bnum.ToString("D")&"/"&bden.ToString("D"))EndIfNextiEndSub'bernoulli_BigIntEndModule'Bernoulli_numbers
B(0)=1/1B(1)=1/2B(2)=1/6B(4)=-1/30B(6)=1/42B(8)=-1/30B(10)=5/66B(12)=-691/2730B(14)=7/6B(16)=-3617/510B(18)=43867/798B(20)=-174611/330B(22)=854513/138B(24)=-236364091/2730B(26)=8553103/6B(28)=-23749461029/870B(30)=8615841276005/14322B(32)=-7709321041217/510B(34)=2577687858367/6B(36)=-26315271553053477373/1919190B(38)=2929993913841559/6B(40)=-261082718496449122051/13530B(42)=1520097643918070802691/1806B(44)=-27833269579301024235023/690B(46)=596451111593912163277961/282B(48)=-5609403368997817686249127547/46410B(50)=495057205241079648212477525/66B(52)=-801165718135489957347924991853/1590B(54)=29149963634884862421418123812691/798B(56)=-2479392929313226753685415739663229/870B(58)=84483613348880041862046775994036021/354B(60)=-1215233140483755572040304994079820246041491/56786730
import"./fmt"forFmtimport"./big"forBigRatvarbernoulli=Fn.new{|n|if(n<0)Fiber.abort("Argument must be non-negative")vara=List.filled(n+1,null)for(min0..n){a[m]=BigRat.new(1,m+1)varj=mwhile(j>=1){a[j-1]=(a[j-1]-a[j])*BigRat.new(j,1)j=j-1}}return(n!=1)?a[0]:-a[0]// 'first' Bernoulli number}for(nin0..60){varb=bernoulli.call(n)if(b!=BigRat.zero)Fmt.print("B($2d) = $44i / $i",n,b.num,b.den)}
B( 0) = 1 / 1B( 1) = -1 / 2B( 2) = 1 / 6B( 4) = -1 / 30B( 6) = 1 / 42B( 8) = -1 / 30B(10) = 5 / 66B(12) = -691 / 2730B(14) = 7 / 6B(16) = -3617 / 510B(18) = 43867 / 798B(20) = -174611 / 330B(22) = 854513 / 138B(24) = -236364091 / 2730B(26) = 8553103 / 6B(28) = -23749461029 / 870B(30) = 8615841276005 / 14322B(32) = -7709321041217 / 510B(34) = 2577687858367 / 6B(36) = -26315271553053477373 / 1919190B(38) = 2929993913841559 / 6B(40) = -261082718496449122051 / 13530B(42) = 1520097643918070802691 / 1806B(44) = -27833269579301024235023 / 690B(46) = 596451111593912163277961 / 282B(48) = -5609403368997817686249127547 / 46410B(50) = 495057205241079648212477525 / 66B(52) = -801165718135489957347924991853 / 1590B(54) = 29149963634884862421418123812691 / 798B(56) = -2479392929313226753685415739663229 / 870B(58) = 84483613348880041862046775994036021 / 354B(60) = -1215233140483755572040304994079820246041491 / 56786730
Uses lib GMP (GNU MP Bignum Library).
class Rational{ // Weenie Rational class, can handle BigInts fcn init(_a,_b){ var a=_a, b=_b; normalize(); } fcn toString{ "%50d / %d".fmt(a,b) } fcn normalize{ // divide a and b by gcd g:= a.gcd(b); a/=g; b/=g; if(b<0){ a=-a; b=-b; } // denominator > 0 self } fcn __opAdd(n){ if(Rational.isChildOf(n)) self(a*n.b + b*n.a, b*n.b); // Rat + Rat else self(b*n + a, b); // Rat + Int } fcn __opSub(n){ self(a*n.b - b*n.a, b*n.b) } // Rat - Rat fcn __opMul(n){ if(Rational.isChildOf(n)) self(a*n.a, b*n.b); // Rat * Rat else self(a*n, b); // Rat * Int } fcn __opDiv(n){ self(a*n.b,b*n.a) } // Rat / Rat}var [const] BN=Import.lib("zklBigNum");// libGMP (GNU MP Bignum Library)fcn B(N){// calculate Bernoulli(n) var A=List.createLong(100,0); // aka static aka not thread safe foreach m in (N+1){ A[m]=Rational(BN(1),BN(m+1)); foreach j in ([m..1, -1]){ A[j-1]= (A[j-1] - A[j])*j; } } A[0]}foreach b in ([0..1].chain([2..60,2])){ println("B(%2d)%s".fmt(b,B(b))) }B( 0) 1 / 1B( 1) 1 / 2B( 2) 1 / 6B( 4) -1 / 30B( 6) 1 / 42B( 8) -1 / 30B(10) 5 / 66B(12) -691 / 2730B(14) 7 / 6B(16) -3617 / 510B(18) 43867 / 798B(20) -174611 / 330B(22) 854513 / 138B(24) -236364091 / 2730B(26) 8553103 / 6B(28) -23749461029 / 870B(30) 8615841276005 / 14322B(32) -7709321041217 / 510B(34) 2577687858367 / 6B(36) -26315271553053477373 / 1919190B(38) 2929993913841559 / 6B(40) -261082718496449122051 / 13530B(42) 1520097643918070802691 / 1806B(44) -27833269579301024235023 / 690B(46) 596451111593912163277961 / 282B(48) -5609403368997817686249127547 / 46410B(50) 495057205241079648212477525 / 66B(52) -801165718135489957347924991853 / 1590B(54) 29149963634884862421418123812691 / 798B(56) -2479392929313226753685415739663229 / 870B(58) 84483613348880041862046775994036021 / 354B(60) -1215233140483755572040304994079820246041491 / 56786730