サイトトップ | Software TipsR tips

統計処理ソフトウェアRについてのTips/Rで数学を

最終更新:2016年 3月 17日 (木曜日)

このページはRの基本操作を説明する。

Topics




リーマン・ゼータ関数まわりのこと

■川端裕人『算数宇宙の冒険 アリスメトリック!』実業之日本社,ISBN 978-4-408-53563-0(Amazon | bk1)を読み,解説を書かれている柏野雄太さんのサポートページを見て触発されたので,Rでリーマン・ゼータ関数まわりのことがどれくらいできるのか試してみた。

■とりあえずゼータ関数を探してみたら,VGAMライブラリ[manual pdf]とgslライブラリ[日本語版マニュアル……このgslというライブラリはGNU Scientific LibraryをRに実装したもののようだ]に入っていた。どちらも関数名はzeta()であった。

■しかし,VGAMでもgslでもリーマン・ゼータ関数は整数あるいは実数に対してしか定義されていないので,解析接続(注:この概念が今一つ良く理解できていないのがいけないんだけれども)ができず,柏野さんのページで示されているような図は作れない。Rでも複素数(complex)の四則演算などはできるのだが(例えば,x <- 1+1iと打ってから,x*2と打てば[1] 2+2iと表示されるし,x^2と打てば[1] 0+2iと表示される),zeta()は複素数を引数にとれず,gslでは虚部が捨てられて計算されてしまうし,VGAMではNAになってしまうのだった。

■実数範囲でできることは,例えばζ(2)=π2/6[式1]を確かめてみるとか,ζ(2),ζ(3),...,ζ(n)とnを増やしていったときの挙動を図示してみるといったことである。これは,Rで普通に試せる。

■最初の例は,ζ(2)=1+1/(2^2)+1/(3^2)+...という無限級数だから,とりあえず一万項くらいやってみるには,

sum(1/(1:10000)^2)

とすればいい。出力結果は,[1] 1.644834となる。一方,円周率はpiとして組み込まれているので,[式1]の右辺は,

pi^2/6

でいい。出力結果は,[1] 1.644934となる。つまり,ちょうど一万分の一の誤差があるわけである。誤差が見えなくなるまで項数を1桁ずつ増やしていったら,千万項が必要なことがわかった。つまり,

sum(1/(1:10000000)^2)

の結果は,[1] 1.644934となった。なお,後で青木先生からご指摘を受けた通り,コンピュータでの数値計算の常識として,小さい方から足さないと計算誤差が出るので,足し算は小さい方からやるべきと考えたら,

sum(1/(10000000:1)^2)

と書く方がいい。けれども,実はRではどちらでも小さい方から計算されているようだ。つまり,明示的に繰り返しループを書くと,

a <- 0; for (i in 10000000:1) a <- a+1/i^2; a

とした方がいいということだ。この場合はoptions(digits=20)くらいにしないと違いがわからないのだけれども,確かに

a <- 0; for (i in 1:10000000) a <- a+1/i^2; a

とすると,結果に微妙な誤差が出る。

■(2009年12月11日追記)黒川信重編著『数学セミナー増刊:リーマン予想がわかる』日本評論社の巻頭,桜井進「高校生からわかる――超入門:リーマン予想」に書かれていた,オイラーが求めた収束が速い関数というのは,本当に収束が速かった。つまり,

sum(1/((1:50)^2*2^(0:49)))+log(2)^2

の出力結果は,options(digits=20)のとき[1] 1.644934066848226となって,完璧にpi^2/6と一致した。

approxzetatwo <- function(x) { sum(1/((1:x)^2*2^(1:x-1)))+log(2)^2 }

と関数定義してから,sapply(1:50,approxzetatwo)とすると,実は41項目で同じ結果が得られていることがわかる。

■計算する項数を増やすにつれて級数がpi^2/6という極限に近付いていく様子を図示するには,例えば,

plot(cumsum(1/(1:1000)^2),type="l",log="x")
lines(c(1,1000),c(pi^2/6,pi^2/6),col="red")

とすれば下図が描ける。ここでは,見やすくするため項数を対数軸表示にした。なお,pi^2/6は,library(gsl)の後でなら,zeta(2)と置き換え可能である。参考までに書いておくと,ここまでの内容は,いつの間にかMaximaと融合し,Gnuplotも内蔵した高機能なものになったEulerというソフトでも,同じような感じで実行できる。

ζ(2)への漸近

■続いて,ゼータ関数の引数を2以上の整数で徐々に大きくしていくとどうなるかを試してみる。もちろん,賢いgslライブラリを使えばplot(2:1000,zeta(2:1000),type="l",log="x")などということも簡単にできるが,定義式通りにx乗分の1を計算させてzeta()並みの精度を出すために一千万項まで計算させてしまうと,ζ(100)まででもハングアップしたかと思うくらいの途方もない時間がかかるのでやめておく方が無難である。

■そこで,計算させる項数を1万項くらいにとどめてζ(2)からζ(100)までの近似をプロットさせるには,べき乗関数powを定義してから,sapplyの結果の行方向の和をとればいいので,Rのコードは次の通りである。

pow<-function(x,y) x^y
plot(2:100,rowSums(sapply(1:10000,pow,-(2:100))),type="l",log="x",xlab="ζ(n)の引数n",ylab="ζ(n)")
require(gsl)
lines(2:100,zeta(2:100),col="red",lty=2)
legend("topright",lty=c(1,2),col=c("black","red"),c("近似","ζ"))
ζ(2)~ζ(100)

■上図のように横軸に2から100まで,縦軸にそれに対応するゼータ関数の近似値をとったグラフが描けるわけである。この状態でlines(2:100,zeta(2:100),col="red",lty=2)とすると,gslライブラリのゼータ関数で計算した値が重ね描きされる。パソコン画面くらいの解像度では,2つのグラフのズレはほとんどわからないくらいである。

■さていよいよ,複素数を引数にとれるゼータ関数だが,webを検索していたら,Javaで書かれたものが見つかった。コンパイル済みのもの(*.jar)しか公開されていないようだが,このソフトのソースを何とか入手してC++に移植すればRから呼べるかもしれない。しかし公開されていないわけだから,たぶんそれは難しいだろう。それで,このページに書かれていた式をそのままRに置き換えてZeta(s)を関数定義したら(但しもちろん無限和はできないので柏野さんに倣って500項の和にして),複素数を受け付けるゼータ関数ができた。もっとも,Zeta(-1)が-0.0833313となってしまうし(本当は-1/12だから-0.08333333となるはずだし,不思議なことにgslライブラリのzeta(-1)はそうなる。なぜ不思議かといえば,ζ(-1)=-1/12は解析接続をしたことで得られる解だからである),Zeta(2)は1.644930074848となってしまうので近似精度は低い。しかも,Rのガンマ関数gamma()が複素数を引数にとれないので実部が負の場合はまだ計算できない。

それでも,実部が正の場合,例えばZeta(1/2+1i)は計算できて, 0.144659-0.6974934iとなる。もし,これが正しければ,柏野さんのサポートサイトに載っていたグラフが作れるはずである。

複素ガンマ関数は,わざわざ作らなくても,gslライブラリにlngamma_complex()という対数をとった形で入っていることがわかった。exp()は複素数を通すので,これで複素数を通すζ関数の近似ができるはずである。また,大浦さんという方が開発されたライブラリにはCとFortranで複素ガンマ関数のコードが提供されているので,これをRに移植すればgslライブラリを使う必要がなくなるように思う。ぱっと見た感じではわかりやすいソースコードなので,そう苦労しなくても移植できそうな気がした。

■大浦さんのコードをRに移植し,複素数を受け付けるゼータ関数改良してみた。500項の和だったところを1万項の和にしたら,Zeta(2)は1.644934となってくれたのだが,Zeta(-1)は実部は-0.0833333となってくれるものの,なぜか虚部がついてしまう。もう少し検討が必要か。

試しにgslライブラリを使って,Zeta()の関数定義中の複素変数のガンマ関数をexp(lngamma_complex(t))とした関数Zeta2()を定義したら,Zeta2(-1)は[1] -0.08333333+0iとなったので,移植したコードで虚部がついてしまったのは,計算精度の問題だったようだ。

■柏野さんのサイトで言及していただいたのにできていないのは悔しいので,とりあえずgslのlngamma_complex(t)を使ったZeta2()で図示してみる。

source("http://minato.sip21c.org/swtips/Zeta.R")
AbsZeta <- function(x) Mod(Zeta2(1/2+x*1i))
curve(AbsZeta(x),0,50)
lines(c(0,50),c(0,0))

とすると,条件判定のところで警告メッセージが出る部分があるが,下図ができる。残念ながら計算精度が低いのか,非自明なゼロ点がきちんとゼロになってくれないのだが,だいたい近い形のグラフにはなる。

Abs(zeta(1/2+t*1i))のカーブ

■(2010年5月19日追記)大きな勘違いをしていたのだが,上記欠点は,計算精度が低いというよりも,curve()だとステップが粗いため,非自明なゼロ点を与えるxが計算されていないということが原因だった。また,条件判定の警告は,Zeta2()関数がベクトル化されていないためであった。一方,Zeta2()を定義する部分は,100000項でなくて10000項でもZeta2(-1)が[1] -0.08333333+0iとなってくれたし,1000項でも[1] -0.08333328+0iとなった。以上を踏まえ,以下のように書き換えたコード(Zeta2.R)を試してみた。

require(gsl)
Zeta2 <- function(s) {
 if (s==1) { ret <- Inf } else if (Re(s)>0) {
  a <- 1/1^s
  for (i in 1000:2) { a <- a + ((-1)^(i-1))/(i^s) }
  ret <- 1/(1-2^(1-s))*a
 }  else {
  t <- 1-s
  ret <- 2*(2*pi)^(-t)*cos(t*pi/2)*exp(lngamma_complex(t))*Zeta2(t)
 }
ret
}

AbsZeta <- function(x) Mod(Zeta2(1/2+x*1i))
x <- 0:50000/1000
y <- sapply(x,AbsZeta)
plot(x,y,type="l")
lines(c(0,50),c(0,0),col="red")

かなりの計算時間がかかったが,下図の通り,ある程度満足のいくグラフになったように思う。速いコンピュータがあれば,Zeta2()の計算の項数を増やしたりAbsZeta()を描画するステップを細かくすれば,もっときれいなグラフが描けるだろう。

Abs(zeta(1/2+t*1i))のカーブ・改良版

多倍長演算(2014年6月1日)

■滅多に使わないため,自分でもよく忘れるので,このページのアーカイヴからここにコピペしておく。

(多倍長演算)Rで多倍長演算をサポートするパッケージは,gmpが便利である。gmpパッケージをロードすると,as.bigz()で多倍長整数,as.bigq()で多倍長実数を扱える(過去のメモ参照)。しかし,高精度計算と表示には,Rmpfr(R Multiple Precision Floating-Point Reliable)というパッケージの方が向いている。πなどの定数も高精度で含まれているからだ(pi,gamma,catalan,log2については,Const("pi",100)などとすれば高精度で得られる)。mpfr()という関数で高精度オブジェクトができ,そのまま計算ができる。なお,Const("pi",100)の結果はasin(mpfr(1,100))*2と一致するが,mpfr(pi,100)とは一致しないことには注意が必要である。mpfrになった値に関数を適用すると自動的に高精度計算になるが,逆では単なる型変換になってしまうので下位の桁は適当に補われるだけになる。100桁のexp(1)が欲しい時は,mpfr(exp(1),100)ではなくて,exp(mpfr(1,100))でなくてはいけない。


リンクと引用について

Correspondence to: minato-nakazawa[at]umin.net.