ベルヌイ数(3)

ベルヌイ数を使うと、tanのテーラー展開を表せる。

 \tan{x} = \frac{\sin{x}}{\cos{x}} = -i\frac{e^{2ix}-1}{e^{2ix}+1}
 = \frac{1}{x}(\frac{2ixe^{2ix}}{e^{2ix}-1}-\frac{4ixe^{4ix}}{e^{4ix}-1}+1)
 = \frac{1}{x}(\sum_{n=0}^{\infty}{\frac{B_n}{n!}((2ix)^n-(4ix)^n)}+1)

計算すると、n=0と1は消えることが分かるので、
残るのは2以上の偶数項である。

 \tan{x} = \sum_{n=1}^{\infty}{(-1)^n2^{2n}(1-2^{2n})\frac{B_2n}{(2n)!}x^{2n-1}}

これを計算すると、


x+1/3x^3+2/15x^5+17/315x^7+62/2835x^9+1382/155925x^11+21844/6081075x^13
+929569/638512875x^15+6404582/10854718875x^17+443861162/1856156927625x^19
+18888466084/194896477400625x^21+113927491862/2900518163668125x^23
+58870668456604/3698160658676859375x^25
+8374643517010684/1298054391195577640625x^27
+689005380505609448/263505041412702261046875x^29
+129848163681107301953/122529844256906551386796875x^31
+1736640792209901647222/4043484860477916195764296875x^33
+418781231495293038913922/2405873491984360136479756640625x^35
+56518638202982204522669764/801155872830791925447758961328125x^37
+32207686319158956594455462/1126482925555250126673224649609375x^39
+1410211493828985228276049834684/121699582862361447435141825020548828125x^41
+516098083439704913515348955653804/109894723324712387033933067993555591796875x^43
+103537504005512749467288942622106408/54397888045732631581796868656810017939453125x^45
+45361105584983995647044252937847808918/58804116977436974739922415018011629392548828125x^47
+87176517890549500795745183943750553204/278845328893007589895761129278958371635634765625x^49



use strict;
use bigint;

my $n = 50;

my @B;
$B[0] = [ 1, 1 ];
for my $i(1..$n) {
next if $i % 2 == 1 && $i != 1;
$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);
}

if($i % 2 == 1) {
next;
}
my $c = rat_copy($B[$i]);
if($i % 4 == 2) {
$c->[0] = -$c->[0];
}
my $a = 2**$i;
$c = rat_multiply($c, $a * (1 - $a));
my $fac = 1;
for my $k(2..$i) {
$fac *= $k;
}
$c = rat_divide($c, $fac);
if($c->[0] > 0) {
if($i != 2) {
print "+";
}
}
if($c->[0] == 1 && $c->[1] == 1) {
;
}
elsif($c->[0] == -1 && $c->[1] == 1) {
print "-";
}
else {
print rat_toString($c);
}
print "x";
if($i != 2) {
print "^" . ($i - 1);
}
}

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_copy {
return [ @{$_[0]} ];
}

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];
}