ベルヌイ数(1)

または、ベルヌーイ数。

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

で定義される(定義は2通りあるが、ここではこの定義を用いる)。

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

だから、

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

を展開して、xn+1の係数を比較すると、

 \sum_{k=0}^{n}{\frac{1}{(n-k+1)!}\frac{B_k}{k!}} = \frac{1}{n!}
 B_n = 1 - \frac{1}{n+1}\sum_{k=0}^{n-1}{_{n+1}C_{k}B_k}

Perlで適当に書いて、ベルヌイ数を50項まで求めた。
(本当は200項まで求めたが割愛)
3以上の奇数の項は省略。

 \frac{xe^x}{e^x-1} - \frac{1}{2}x = \frac{1}{2}x\frac{e^x+1}{e^x-1}

が偶関数で0になるから。


B0 = 1
B1 = 1/2
B2 = 1/6
B4 = -1/30
B6 = 1/42
B8 = -1/30
B10 = 5/66
B12 = -691/2730
B14 = 7/6
B16 = -3617/510
B18 = 43867/798
B20 = -174611/330
B22 = 854513/138
B24 = -236364091/2730
B26 = 8553103/6
B28 = -23749461029/870
B30 = 8615841276005/14322
B32 = -7709321041217/510
B34 = 2577687858367/6
B36 = -26315271553053477373/1919190
B38 = 2929993913841559/6
B40 = -261082718496449122051/13530
B42 = 1520097643918070802691/1806
B44 = -27833269579301024235023/690
B46 = 596451111593912163277961/282
B48 = -5609403368997817686249127547/46410
B50 = 495057205241079648212477525/66

よく見たら、最初以外は(奇数)/(偶数)になっている。


以下は、ソース。



use strict;
use bigint;

my $n = 200;

my @B;
$B[0] = [ 1, 1 ];
for my $i(1..$n) {
$B[$i] = [ 1, 1 ];
my $C = [ 1, $i + 1 ];
for(my $k = 0; $k < $i; $k++) {
$B[$i] = rat_subtract($B[$i], rat_multiply($B[$k], $C));
my $tmp = [ $i - $k + 1, $k + 1 ];
rat_normalize($tmp);
$C = rat_multiply($C, $tmp);
}
printf "B%d = %s\n", $i, rat_toString($B[$i]);
}

sub factorial {
if($_[0] == 0) {
return 1;
}
else {
return factorial($_[0] - 1) * $_[0];
}
}

sub rat_add {
my @r;
my ($a, $d) = @_;
if(ref $a eq "ARRAY") {
if(ref $d eq "ARRAY") {
$r[0] = $a->[0] * $d->[1] + $a->[1] * $d->[0];
$r[1] = $a->[1] * $d->[1];
rat_normalize(\@r);
}
else {
$r[0] = $a->[0] + $d * $a->[1];
$r[1] = $a->[1];
}
}
else {
if(ref $d eq "ARRAY") {
$r[0] = $a * $d->[1] + $d->[0];
$r[1] = $d->[1];
}
else {
$r[0] = $a + $d;
$r[1] = 1;
}
}
return \@r;
}

sub rat_subtract {
my @r;
my ($a, $d) = @_;
if(!$a) {
$a->[0] = 0;
$a->[1] = 1;
}
if(!$d) {
$d->[0] = 0;
$d->[1] = 1;
}
if(ref $a eq "ARRAY") {
if(ref $d eq "ARRAY") {
$r[0] = $a->[0] * $d->[1] - $a->[1] * $d->[0];
$r[1] = $a->[1] * $d->[1];
rat_normalize(\@r);
}
else {
$r[0] = $a->[0] - $d * $a->[1];
$r[1] = $a->[1];
}
}
else {
if(ref $d eq "ARRAY") {
$r[0] = $a * $d->[1] - $d->[0];
$r[1] = $d->[1];
}
else {
$r[0] = $a - $d;
$r[1] = 1;
}
}
return \@r;
}

sub rat_multiply {
my @r;
my ($a, $d) = @_;
if(!$a) {
$a->[0] = 0;
$a->[1] = 1;
}
if(!$d) {
$d->[0] = 0;
$d->[1] = 1;
}

if(ref $a eq "ARRAY") {
if(ref $d eq "ARRAY") {
$r[0] = $a->[0] * $d->[0];
$r[1] = $a->[1] * $d->[1];
}
else {
$r[0] = $a->[0] * $d;
$r[1] = $a->[1];
}
}
else {
if(ref $d eq "ARRAY") {
$r[0] = $a * $d->[0];
$r[1] = $d->[1];
}
else {
$r[0] = $a * $d;
$r[1] = 1;
}
}
rat_normalize(\@r);
return \@r;
}

sub rat_divide {
my @r;
my ($a, $d) = @_;
if(ref $a eq "ARRAY") {
if(ref $d eq "ARRAY") {
$r[0] = $a->[0] * $d->[1];
$r[1] = $a->[1] * $d->[0];
}
else {
$r[0] = $a->[0];
$r[1] = $a->[1] * $d;
}
}
else {
if(ref $d eq "ARRAY") {
$r[0] = $a * $d->[1];
$r[1] = $d->[0];
}
else {
$r[0] = $a;
$r[1] = $d;
}
}
rat_normalize(\@r);
return \@r;
}

sub rat_normalize {
my ($a) = @_;
my $a1 = abs($a->[0]);
my $b1 = abs($a->[1]);
my ($a2, $f);
while(1) {
$a2 = $a1 % $b1;
if($a2 == 0) {
$f = $b1;
last;
}
$b1 %= $a1;
if($b1 == 0) {
$f = $a1;
last;
}
$a1 = $a2;
}
$a->[0] /= $f;
$a->[1] /= $f;
if($a->[1] < 0) {
$a->[0] = -$a->[0];
$a->[1] = -$a->[1];
}
}

sub rat_toString {
my ($a) = @_;
if($a->[1] == 1) {
return $a->[0];
}
else {
return $a->[0] . "/" . $a->[1];
}
}

sub max {
return $_[0] > $_[1] ? $_[0] : $_[1];
}

sub min {
return $_[0] > $_[1] ? $_[1] : $_[0];
}