プログラミング初心者でもベルヌーイ数を計算したい - 音楽室と化学室と美術室とPC室の融合部屋　所謂自由室

$\frac{x}{e^x - 1}$
の0を中心としたテイラー展開
$\frac{x}{e^x - 1} = \sum_{n=0}^{\infty} \frac{B_n}{n!}x^n$
について項ごとに表れる $B_n$ はベルヌーイ数と言われる．
このベルヌーイ数は例えばリーマンゼータ関数
$\zeta (s) = \sum_{n=1}^{\infty}\frac{1}{n^s}$
についてsが正の偶数のときに
$\zeta(2n) = (-1)^{n+1} \frac{B_{2n}\cdot (2π)^{2n}}{2\cdot(2n)!}$
と，いった形で現れたりする．
今回はこのベルヌーイ数を計算させようと思ったが，素直にテイラー展開を微分から求めようとすることは非常に困難なので，漸化式から求めることにする．

この漸化式の導出は

$\frac{x}{e^x - 1} \cdot \frac{e^x - 1}{x} = 1$

という恒等式を作り

$\frac{x}{e^x - 1} = \sum_{n=0}^{\infty} \frac{B_n}{n!}x^n$

及び

$\frac{e^x - 1}{x} = \frac{\sum_{n=1}^{\infty} \frac{x^n}{n!}}{x}$

より

$\frac{e^x - 1}{x} = \sum_{n=0}^{\infty} \frac{x^n}{(n+1)!}$

という式変形ができるので

$\frac{x}{e^x - 1} \cdot \frac{e^x - 1}{x} = \left(\sum_{n=0}^{\infty} \frac{B_n}{n!}x^n\right)\left(\sum_{m=0}^{\infty} \frac{x^m}{(m+1)!}\right) = 1$

となり，ピネ・コーシーの恒等式より

$\sum_{n=0}^{\infty} \sum_{m=0}^{\infty} \frac{B_n}{n!(m+1)!}x^{n+m}= 1$

と変形ができる．

$x^{n+m}$ の指数を1つにまとめるために $u = n+m$ ， $v = n$ とすると， $n+m = u$ とおいているため $u\geq v$ となるので

$\sum_{u=0}^{\infty}\left(\sum_{v=0}^{u}\frac{B_v}{v!(u-v+1)!}\right)x^u = 1$

になる．

これが恒等式であるならば $x^u$ について， $u=0$ のときは括弧の中は $1$ ， $u\neq 0$ のときは括弧の中は $0$ にならないといけない．

$u=0$ のときは

$\sum_{v=0}^{0}\frac{B_v}{v!(-v+1)!} = 1 = \frac{B_0}{0!\cdot 1!} = B_0 = 1$

となる一方 $u\gt 0$ のときは

$\sum_{v=0}^u\frac{B_v}{v!(u-v+1)!} = \sum_{v=0}^{u-1}\frac{B_v}{v!(u-v+1)!} + \frac{B_u}{u!1!} =0$

とすると

$B_u = -\sum_{v=0}^{u-1}\frac{u! \cdot B_v}{v! ( (u+1)-v)!}$

$B_u = -\sum_{v=0}^{u-1}\frac{(u+1)!}{v! ( (u+1)-v)!}\cdot \frac{B_v}{u+1}$

になるので

$B_u = -\frac{1}{u+1}\sum_{v=0}^{u-1}\left(_{n+1}C_v\cdot B_v\right)$

が得られる．

この漸化式をC++で組み計算を行った．

#include <iostream>

#include <boost/multiprecision/cpp_int.hpp>

#include <boost/multiprecision/cpp_dec_float.hpp>

#define NUMERIC 14000

using namespace boost::multiprecision;

typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<NUMERIC>> cpp_u_float;

cpp_int abs(cpp_int a) { return (a > 0) ? a : -a;}

cpp_int gcd(cpp_int a,cpp_int b) {

if (a < b) { cpp_int tmp = a; a = b; b = tmp; }

cpp_int r = a%b;

while (r != 0) { a = b; b = r; r = a % b; }

return b;

}

void Bernoulli(unsigned long long N) {

std::vector<std::vector<cpp_int>> List;

std::vector<cpp_int> TList(2);

TList[0] = TList[1] = 1;//0:分子　1:分母

List.push_back(TList);

for (unsigned long long i = 1; i <= N; i++) {

cpp_int numer = List[0][0] * Combination(i + 1, 0);

cpp_int denom = List[0][1];

for (unsigned long long j = 1; j < i; j++) {

numer = numer * List[j][1] + denom * List[j][0] * Combination(i + 1, j);

denom = denom * List[j][1];

if (numer != 0) {

cpp_int g = gcd(numer, denom);

numer /= g;

denom /= g;

}

numer *= -1;

denom *= i + 1;

if ( (numer < 0 && denom < 0) || (numer > 0 && denom < 0)) {

numer = -numer;

denom = -denom;

}

if (numer == 0)denom = 1;

else {

cpp_int g = gcd(abs(numer), abs(denom));

numer /= g;

denom /= g;

}

TList[0] = numer;

TList[1] = denom;

List.push_back(TList);

}

if (List[N][0] != 0) std::cout << "B_" << std::to_string(N) << " = " <<List[N][0] << " / " << List[N][1] << std::endl;

else std::cout << "B_" << std::to_string(N) << " = " << List[N][0] << std::endl;

}

以下200までの出力

B_0 = 1 / 1
B_1 = -1 / 2
B_2 = 1 / 6
B_3 = 0
B_4 = -1 / 30
B_5 = 0
B_6 = 1 / 42
B_7 = 0
B_8 = -1 / 30
B_9 = 0
B_10 = 5 / 66
B_11 = 0
B_12 = -691 / 2730
B_13 = 0
B_14 = 7 / 6
B_15 = 0
B_16 = -3617 / 510
B_17 = 0
B_18 = 43867 / 798
B_19 = 0
B_20 = -174611 / 330
B_21 = 0
B_22 = 854513 / 138
B_23 = 0
B_24 = -236364091 / 2730
B_25 = 0
B_26 = 8553103 / 6
B_27 = 0
B_28 = -23749461029 / 870
B_29 = 0
B_30 = 8615841276005 / 14322
B_31 = 0
B_32 = -7709321041217 / 510
B_33 = 0
B_34 = 2577687858367 / 6
B_35 = 0
B_36 = -26315271553053477373 / 1919190
B_37 = 0
B_38 = 2929993913841559 / 6
B_39 = 0
B_40 = -261082718496449122051 / 13530
B_41 = 0
B_42 = 1520097643918070802691 / 1806
B_43 = 0
B_44 = -27833269579301024235023 / 690
B_45 = 0
B_46 = 596451111593912163277961 / 282
B_47 = 0
B_48 = -5609403368997817686249127547 / 46410
B_49 = 0
B_50 = 495057205241079648212477525 / 66
B_51 = 0
B_52 = -801165718135489957347924991853 / 1590
B_53 = 0
B_54 = 29149963634884862421418123812691 / 798
B_55 = 0
B_56 = -2479392929313226753685415739663229 / 870
B_57 = 0
B_58 = 84483613348880041862046775994036021 / 354
B_59 = 0
B_60 = -1215233140483755572040304994079820246041491 / 56786730
B_61 = 0
B_62 = 12300585434086858541953039857403386151 / 6
B_63 = 0
B_64 = -106783830147866529886385444979142647942017 / 510
B_65 = 0
B_66 = 1472600022126335654051619428551932342241899101 / 64722
B_67 = 0
B_68 = -78773130858718728141909149208474606244347001 / 30
B_69 = 0
B_70 = 1505381347333367003803076567377857208511438160235 / 4686
B_71 = 0
B_72 = -5827954961669944110438277244641067365282488301844260429 / 140100870
B_73 = 0
B_74 = 34152417289221168014330073731472635186688307783087 / 6
B_75 = 0
B_76 = -24655088825935372707687196040585199904365267828865801 / 30
B_77 = 0
B_78 = 414846365575400828295179035549542073492199375372400483487 / 3318
B_79 = 0
B_80 = -4603784299479457646935574969019046849794257872751288919656867 / 230010
B_81 = 0
B_82 = 1677014149185145836823154509786269900207736027570253414881613 / 498
B_83 = 0
B_84 = -2024576195935290360231131160111731009989917391198090877281083932477 / 3404310
B_85 = 0
B_86 = 660714619417678653573847847426261496277830686653388931761996983 / 6
B_87 = 0
B_88 = -1311426488674017507995511424019311843345750275572028644296919890574047 / 61410
B_89 = 0
B_90 = 1179057279021082799884123351249215083775254949669647116231545215727922535 / 272118
B_91 = 0
B_92 = -1295585948207537527989427828538576749659341483719435143023316326829946247 / 1410
B_93 = 0
B_94 = 1220813806579744469607301679413201203958508415202696621436215105284649447 / 6
B_95 = 0
B_96 = -211600449597266513097597728109824233673043954389060234150638733420050668349987259 / 4501770
B_97 = 0
B_98 = 67908260672905495624051117546403605607342195728504487509073961249992947058239 / 6
B_99 = 0
B_100 = -94598037819122125295227433069493721872702841533066936133385696204311395415197247711 / 33330
B_101 = 0
B_102 = 3204019410860907078243020782116241775491817197152717450679002501086861530836678158791 / 4326
B_103 = 0
B_104 = -319533631363830011287103352796174274671189606078272738327103470162849568365549721224053 / 1590
B_105 = 0
B_106 = 36373903172617414408151820151593427169231298640581690038930816378281879873386202346572901 / 642
B_107 = 0
B_108 = -3469342247847828789552088659323852541399766785760491146870005891371501266319724897592306597338057 / 209191710
B_109 = 0
B_110 = 7645992940484742892248134246724347500528752413412307906683593870759797606269585779977930217515 / 1518
B_111 = 0
B_112 = -2650879602155099713352597214685162014443151499192509896451788427680966756514875515366781203552600109 / 1671270
B_113 = 0
B_114 = 21737832319369163333310761086652991475721156679090831360806110114933605484234593650904188618562649 / 42
B_115 = 0
B_116 = -309553916571842976912513458033841416869004128064329844245504045721008957524571968271388199595754752259 / 1770
B_117 = 0
B_118 = 366963119969713111534947151585585006684606361080699204301059440676414485045806461889371776354517095799 / 6
B_119 = 0
B_120 = -51507486535079109061843996857849983274095170353262675213092869167199297474922985358811329367077682677803282070131 / 2328255930
B_121 = 0
B_122 = 49633666079262581912532637475990757438722790311060139770309311793150683214100431329033113678098037968564431 / 6
B_123 = 0
B_124 = -95876775334247128750774903107542444620578830013297336819553512729358593354435944413631943610268472689094609001 / 30
B_125 = 0
B_126 = 5556330281949274850616324408918951380525567307126747246796782304333594286400508981287241419934529638692081513802696639 / 4357878
B_127 = 0
B_128 = -267754707742548082886954405585282394779291459592551740629978686063357792734863530145362663093519862048495908453718017 / 510
B_129 = 0
B_130 = 1928215175136130915645299522271596435307611010164728458783733020528548622403504078595174411693893882739334735142562418015 / 8646
B_131 = 0
B_132 = -410951945846993378209020486523571938123258077870477502433469747962650070754704863812646392801863686694106805747335370312946831 / 4206930
B_133 = 0
B_134 = 264590171870717725633635737248879015151254525593168688411918554840667765591690540727987316391252434348664694639349484190167 / 6
B_135 = 0
B_136 = -84290226343367405131287578060366193649336612397547435767189206912230442242628212786558235455817749737691517685781164837036649737 / 4110
B_137 = 0
B_138 = 2694866548990880936043851683724113040849078494664282483862150893060478501559546243423633375693325757795709438325907154973590288136429 / 274386
B_139 = 0
B_140 = -3289490986435898803930699548851884006880537476931130981307467085162504802973618096693859598125274741604181467826651144393874696601946049 / 679470
B_141 = 0
B_142 = 14731853280888589565870080442453214239804217023990642676194878997407546061581643106569966189211748270209483494554402556608073385149191 / 6
B_143 = 0
B_144 = -3050244698373607565035155836901726357405007104256566761884191852434851033744761276392695669329626855965183503295793517411526056244431024612640493 / 2381714790
B_145 = 0
B_146 = 4120570026280114871526113315907864026165545608808541153973817680034790262683524284855810008621905238290240143481403022987037271683989824863 / 6
B_147 = 0
B_148 = -1691737145614018979865561095112166189607682852147301400816480675916957871178648433284821493606361235973346584667336181793937950344828557898347149 / 4470
B_149 = 0
B_150 = 463365579389162741443284425811806264982233725425295799852299807325379315501572305760030594769688296308375193913787703707693010224101613904227979066275 / 2162622
B_151 = 0
B_152 = -3737018141155108502105892888491282165837489531488932951768507127182409731328472084456653639812530140212355374618917309552824925858430886313795805601 / 30
B_153 = 0
B_154 = 10259718682038021051027794238379184461025738652460569233992776489750881337506863808448685054322627708245455888249006715516690124228801409697850408284121 / 138
B_155 = 0
B_156 = -81718086083262628510756459753673452313595710396116467582152090596092548699138346942995509488284650803976836337164670494733866559829768848363506624334818961419869 / 1794590070
B_157 = 0
B_158 = 171672676901153210072183083506103395137513922274029564150500135265308148197358551999205867870374013289728260984269623579880772408522396975250682773558018919 / 6
B_159 = 0
B_160 = -4240860794203310376065563492361156949989398087086373214710625778458441940477839981850928830420029285687066701804645453159767402961229305942765784122421197736180867 / 230010
B_161 = 0
B_162 = 1584451495144416428390934243279426140836596476080786316960222380784239380974799880364363647978168634590418215854419793716549388865905348534375629928732008786233507729 / 130074
B_163 = 0
B_164 = -20538064609143216265571979586692646837805331023148645068133372383930344948316600591203926388540940814833173322793804325084945094828524860626092013547281335356200073083 / 2490
B_165 = 0
B_166 = 5734032969370860921631095311392645731505222358555208498573088911303001784652122964703205752709194193095246308611264121678834250704468082648313788124754168671815815821441 / 1002
B_167 = 0
B_168 = -13844828515176396081238346585063517228531109156984345249260453934317772754836791258987516540324983611569758649525983347408589045734176589270143058509026392246407576578281097477 / 3404310
B_169 = 0
B_170 = 195334207626637530414976779238462234481410337350988427215139995707346979124686918267688171536352650572535330369818176979951931477427594872783018749894699157917782460035894085 / 66
B_171 = 0
B_172 = -11443702211333328447187179942991846613008046506032421731755258148665287832264931024781365962633301701773088470841621804328201008020129996955549467573217659587609679405537739509973 / 5190
B_173 = 0
B_174 = 4166161554662042831884959593250717297395614318182561412048180684077407803317591270831194619293832107482426945655143357909807251852859279483176373435697607639883085093246499347128331 / 2478
B_175 = 0
B_176 = -1369347910486705707645621362512824332220360774476594348356938715366608044588614657557436131706543948464159947970464346070253278291989696390096800799614617317655510118710460076077638883999 / 1043970
B_177 = 0
B_178 = 1124251816617941290026484851206299982774720467712867275292043701618829826708395745459654170718363182143418314514085426692857018428614935412736063946853033094328968069656979232446257101741 / 1074
B_179 = 0
B_180 = -6173136454016248924640522272263470960199559328290655337530202055853397791747341312347030141906500993752700612233695954532816018207721731818225290076670213481102834647254685911917265818955932383093313 / 7225713885390
B_181 = 0
B_182 = 4277269279349192541137304400628629348327468135828402291661683018622451659989595510712915810436238721139546963558655260384328988773219688091443529626531335687951612545946030357929306651006711 / 6
B_183 = 0
B_184 = -857321333523056180131194437347933216431403305730705359015465649285681432317514010686029079324479659634642384809061711319481020030715989009140595170556956196762318625529645723516532076273012244047 / 1410
B_185 = 0
B_186 = 22258646098436968050639602221816385181596567918515338169946670500599612225742487595012775838387331550474751212260636163500086787417640903770807353228157478339547041472679880890292167353534100797481 / 42
B_187 = 0
B_188 = -14158277750623758793309386870401397333112823632717478051426522029712001260747920789473711562165031101665618225654329210473605281619696918061316240634857984019071572591940586875558943580878119388321001 / 30
B_189 = 0
B_190 = 5411555842544259796131885546196787277987837486638756184149141588783989774511509608733429067517383750706299486822702171672522203106730993581242777825864203487238429479957280273093904025319950569633979493395 / 12606
B_191 = 0
B_192 = -346465752997582699690191405750952366871923192340955593486485715370392154894102000406980162521728492501917598012711402163530166516991115122131398542029056286959857727373568402417020319761912636411646719477318166587 / 868841610
B_193 = 0
B_194 = 2269186825161532962833665086968359967389321429297588337232986752409765414223476696863199759981611817660735753831323900456495253961837175924312108872915089534970310604331636484174526399721365966337809334021247 / 6
B_195 = 0
B_196 = -62753135110461193672553106699893713603153054153311895305590639107017824640241378480484625554578576142115835788960865534532214560982925549798683762705231316611716668749347221458005671217067357943416524984438771831113 / 171390
B_197 = 0
B_198 = 88527914861348004968400581010530565220544526400339548429439843908721196349579494069282285662653465989920237253162555666526385826449862863083834096823053048072002986184254693991336699593468906111158296442729034119206322233 / 244713882
B_199 = 0
B_200 = -498384049428333414764928632140399662108495887457206674968055822617263669621523687568865802302210999132601412697613279391058654527145340515840099290478026350382802884371712359337984274122861159800280019110197888555893671151 / 1366530