ポセットのメビウス関数とオイラー標数

  • このPDF
  • かいつまむと
  • ポセットのチェイン(順序がたどれるパス)を単体とみなす(たとえば、x1,x2,x3がチェインなら、x1,x2,x3の軸上の点(1,0,0),(0,1,0),(0,0,1)を結んだ三角形がこのチェインの表す単体
  • ポセット全体は構成ノード数の次元の空間に置かれた単体的複体となる
  • この単体的複体は山あり谷あり
  • ポセットの各チェインに対して1を返す関数を、数論的関数(整数に対して複素数を返す関数)とすると、その数論的関数を係数とするディリクレ級数的母関数は、ゼータ関数に相当する(普通の自然数に対して、必ず1を返すような数論的関数に対応する、ディリクレ級数的母関数がリーマンのゼータ関数である)
  • ゼータ関数と、ディリクレ級数的畳み込みにおいて、逆関数になっているのがメビウス関数(メビウス関数も数論的関数。今はポセットを考えているから、各チェインに対して値を与える関数になっている)
  • メビウス関数を、ポセットのような、集合のinclusion-exclusion関係に対して考えるとき、+1,-1の交代が現れる
  • 今、ポセットの一番下と一番上を閉じてやると、メビウス関数\mu = \zeta^{-1}(\zeta-1)^kの項に分解した表現ができて、その各項に係数がつく
  • この係数がオイラー標数(点の数、辺の数、面の数、四面体の数…、という構成単体の規模別の数)の±1倍となっており、メビウス関数自体が、オイラー標数の式に一致することが示せるという

京都大学(一般入試)数学の問題を計算機で解く(2019)


問題と解答例はこちらなどで

第一問

  • 問1
    • \cos{\theta}無理数\cos{2\theta}\cos{3\theta}とが有理数となるような\theta
    • この手の計算の桁落ちなどの影響が出るタイプの問題は、整数割る整数などを、整数を用いたまま計算するタイプのパッケージなどを使わないとうまくいかないが、あえてやってみる
    • \cos{2\theta}\cos{3\theta}とが有理数となる場合のグリッド探索
    • 答えは\theta=\pi/6
# 有理数 q = m/n
# 適当にたくさんの0以上1以下(cos(2t),cos(3t)の値なので)の有理数を作る。1/3,2/6なども区別せずにさっさと作る
N <- 100
n <- 1:N
m <- 0:N
mn <- expand.grid(m,n)
mn <- mn[which(mn[,1] < mn[,2]),]
q <- mn[,1]/mn[,2]
# 対応する角度を計算する
phi <- acos(q)
# 対応する角度に2t,3tとなるペアがあるかどうかは、3倍した角度と2倍した角度が一致するかで見れば良い
phi.pair <- expand.grid(phi,phi)
theta.cand <- which(abs(phi.pair[,1] * 3 - phi.pair[,2] * 2) ==0)
theta.cand.pi <- phi.pair[theta.cand,]/pi
unique(theta.cand.pi)
# 2t = pi/3, 3t = pi/2 が候補であることがわかる
> unique(theta.cand.pi)
       Var1 Var2
3 0.3333333  0.5
  • 問2
    • 積分計算。計算機には計算機の流儀がある
      • (1) \int_0^{\frac{\pi}{4}} \frac{x}{\cos^2{x}} dx
        • 答えは\frac{\pi}{4}-\frac{1}{2}\log{2}
      • (2)\int_0^{\frac{\pi}{4}} \frac{dx}{\cos{x}}
        • 答えは\log{\sqrt{2}+1}
f1 <- function(x){
  x/cos(x)^2
}
integrate(f1,upper=pi/4,lower=0)
pi/4 - 1/2 * log(2)
> integrate(f1,upper=pi/4,lower=0)
0.4388246 with absolute error < 4.9e-15
> pi/4 - 1/2 * log(2)
[1] 0.4388246
f2 <- function(x){
  1/cos(x)
}
integrate(f2,upper=pi/4,lower=0)
log(sqrt(2)+1)
> integrate(f2,upper=pi/4,lower=0)
0.8813736 with absolute error < 9.8e-15
> log(sqrt(2)+1)
[1] 0.8813736

第二問

  • またもや、整数・素数問題…(京大の)数学科が計算機科学と違う方向を向いていることを感じさせます
  • f(n) = n^3 + 2 n^2 + 2の絶対値が連続して素数になるような整数nを求める問題
  • 答えは-3,-2,-1,0
  • 列挙でやることにします
# 素数を扱う数論のパッケージをインストールして読み込む
install.packages("numbers")
library(numbers)
isPrime(1:10) # 1から10が素数かどうかの判定
# たくさんの整数を並べる
n <- (-1000):1000
# 指定の関数の値を計算する
fn <- n^3 + 2 * n^2 + 2
absfn <- abs(fn)
# 素数であるものの番地を取り出す
oks <- which(isPrime(absfn))
# |f(n)|と|f(n+1)|とが共に素数なので、連続する番地がoksに含まれるはず
diff.oks <- diff(oks) # この値が1になって入れば、absfnの隣同士が素数だったことがわかる
ans <- which(diff.oks==1)
oks[ans]
# これが答え
n[oks[ans]]
absfn[oks[ans]]
> oks[ans]
[1]  998  999 1000 1001
> 
> n[oks[ans]] # 答え
[1] -3 -2 -1  0
> absfn[oks[ans]] # |f(n)|の値
[1] 7 2 3 2

第三問

  • 鋭角三角形の内部の曲線を1パラメタで指定し、その曲線が区切る領域の面積と三角形の面積との比に関する問題
  • 鋭角三角形の頂点Bを原点に、Cを(1,0)とし、頂点Aの座標を(r,r\tan{\theta}),r>0,0 < \theta < \pi/2としても一般性を失わない
  • まず描図

f:id:ryamada:20190226100805j:plain

  • 赤い曲線の下の面積が三角形の面積に占める割合を計算機的に求めるために、三角形内部に一様な乱点を発生させ、その点が曲線の上にあるか下にあるかで分類する
  • 曲線が正確にも止まっているわけではないので、曲線を構成する点のうち、乱点とx座標が最も近い点を見つけ、それのy座標と乱点のy座標を比較することで、曲線の上下の判定をすることは、点の数が多ければまあ、良いだろうと言う作戦にする
  • 割合が答え(1/3)に近いことを示す
# 辺BCの長さを1に固定して、鋭角三角形を一般的に作成し、プロットする
B <- c(0,0)
C <- c(1,0)
r <- runif(1)
theta <- runif(1) * pi/2
A <- c(r,r*tan(theta))

plot(rbind(B,C,A))
segments(c(B[1],C[1],A[1]),c(B[2],C[2],A[2]),c(C[1],A[1],B[1]),c(C[2],A[2],B[2]))
# パラメタtを用いてくだんの曲線を描く
# parameter t
t <- seq(from=0,to=1,length=10000)
Q <- cbind(t * (C-B)[1]+ (1-t) * (A-B)[1], t*(C-B)[2] + (1-t) * (A-B)[2])
P <- t * Q
points(P,col=2,type="l")

# 三角形内部の一様乱点はディリクレ乱数で作ることにする
library(MCMCpack)
N <- 10000
r <- rdirichlet(N,rep(1,3))
R <- cbind(A[1] * r[,1] + B[1] * r[,2] + C[1] * r[,3], A[2] * r[,1] + B[2] * r[,2] + C[2] * r[,3])
# 乱点を表示する
points(R,pch=20,cex=0.5,col=4)

# Rのx座標と一番近い、曲線のx座標を探す
# そしてy座標を比較し、Rのそれが曲線のそれより小さければ、対象領域にあると判断する
# 曲線との上下関係を逐一、判定する
inOut <- rep(0,N)
for(i in 1:N){
  rx <- R[i,1]
# x座標が最も近い曲線構成点を見つける
  tmp <- which(abs(rx-P[,1])==min(abs(rx-P[,1])))
  ry <- R[i,2]
  neighbory <- P[tmp,2]
# y座標を比較する
  if(ry < neighbory){
    inOut[i] <- 1
  }
}
sum(inOut)/N # この割合が求めたい値。確かに約1/3
plot(rbind(B,C,A))
segments(c(B[1],C[1],A[1]),c(B[2],C[2],A[2]),c(C[1],A[1],B[1]),c(C[2],A[2],B[2]))
points(P,col=2,type="l")
points(R,pch=20,cex=0.5,col=4+inOut)

f:id:ryamada:20190226103957j:plain

> sum(inOut)/N # この割合が求めたい値。確かに約1/3
[1] 0.3332

第四問

  • 確率事象は計算機には乱数実験なので楽しい
  • サイコロを何度も振りつつ、その途中にある回数目までは4以下しか出ず、その次に初めて5以上がでるようなことは、サイコロを振る回数増えれば、ほとんど起きなくなるのは予想がつく
  • \frac{((n-1)*2^n+1}{3^n}が答えだと言う
  • Rで計算機実験する
n <- 50 # サイコロを続けて振る回数
N.iter <- 1000 # 実験を繰り返す回数
# サイコロの目をランダムに発生
Xs <- matrix(sample(1:6,n*N.iter,replace=TRUE,prob=rep(1/6,6)),ncol=n)
# 0回目に振った時の値を0とする、と言う問題条件を追加
Xs0 <- cbind(rep(0,N.iter),Xs)
# k-1回目まで…、k回目に…、と言う条件に合うかの判定
meet.cond <- matrix(0,N.iter,n)

for(i in 1:N.iter){
  for(k in 1:n){
    Xk_1 <- Xs0[i,k]
    Xk <- Xs[i,k]
    if(Xk_1 <= 4 & Xk >=5){
      meet.cond[i,k] <- 1
    }
  }
}
# n回振って、その経過中に1つのkが条件を満足する、と言うことを計算するために
# 条件を満足した回数の累和をとる
cum.meet.cond <- t(apply(meet.cond,1,cumsum))
# 都合、1度きり、の条件集計
num.ok.events <- apply(cum.meet.cond,2,function(x){length(which(x==1))})
# 以下が、計算機実験での、確率(振り続けた回数ごとに算出)
prob.ans <- num.ok.events/N.iter 

# いわゆる数学の試験的な答えとしての式があるので、それで計算してやる
nn <- 1:n
exact.prob <- ((nn-1)*2^(nn)+1)/3^nn
# 正確確率と実験的確率の比較
plot(1:n,exact.prob,type="h") # 棒が正確確率
points(1:n,prob.ans,pch=20,col=2) # 赤丸が計算機実験値

f:id:ryamada:20190226112456j:plain

第五問

  • 球面上の5点が作る四角錐の体積。底面は正方形
  • 底面をなす4点のz座標を同じにして話をすることにする。その座標をzbとすれば、正方形の面積は\sqrt{1-|zb|^2} ^2 \times 2
  • 底面からの高さは高ければ高いほど体積は大きくなるので、z軸上にピラミッドの頂上の点をとることにする。その点の座標を(0,0,1)と固定し、zbを-1から1まで動かそう
zb <- seq(from=-1,to=1,length=1000)

a <- 2 * (1-abs(zb)^2)
h <- 1-zb

# a*hの最大値をもたらすzbがわかれば良い
plot(zb,a*h,type="l")
abline(v=-1/3)
abline(h=max(a*h))
  • 答えはz=-1/3に正方形がある場合で、確かにそうなっている。実際の体積計算には、錐の体積計算のための係数を掛けるなどする必要がある

f:id:ryamada:20190226114516j:plain

第六問

  • (1+i)^n + (1-i)^n > 10^{10}を満たす最小の正の整数…
  • 馬鹿でかい数だけれども、最近の計算機はこのくらいはものともしない
  • 左辺の虚部が消えることから、計算して、実部だけを取り出そう
  • 全列挙作戦
n <- 1:100
v <- Re((1+1i)^n + (1-1i)^n)
which(v > 10^10)
# 絶対値を底10の対数にする方がプロットしやすいので
logv <- sign(v) * log10(abs(v))
plot(logv,type="h")
abline(h=10,col=2)

f:id:ryamada:20190226120008j:plain

  • 答えは71
which(v > 10^10)
 [1] 71 72 73 79 80 81 87 88 89 95 96 97

余代数 Coalgebra

  • Coalgebra, 余代数、って何?
  • 代数は、集合があって、集合の要素x,yに演算を定めて、x _times y =zとなった時のaが集合に含まれる
  • 余代数はその逆で、集合の要素xがあった時に、演算が定めてあって、xによって、集合の要素のどれをどう使うかが決まっていて、その結果、現れるものが、やはり集合の要素である、と言うようなもの
  • ポセットのインターバルに余代数が入る、と言うのは、あるインターバルが決まると、それに帰属するインターバルたちが決まり、それらに対してどう言う演算をするかが決まっているのだが、その演算の結果が、やっぱりインターバルだ、と言うそう言うことらしい…
  • こちらを参考

接続ホップ代数 Incidence Hopf algebras

この文書(Incidence Hopf algebras)を読んでみる

ディリクレ級数

  • 1以上の自然数を考える
  • 自然数を引数とする複素関数f(n)を考える
  • F(s) = \sum_{n \ge 1} \frac{f(n)}{n^s}をディリクレ級数と言う。sは複素数だから、これは複素関数
    • 関数f(n)の性質を複素関数的に表現したもの、と言う意味にとっても良い
  • リーマンの\zeta関数はf(n)=1のことで、それのディリクレ級数
    • \zeta(s) = \sum_{n \ge 1} \frac{1}{n^s} = \prod_{p prime} \frac{1}{1-1/p^s}と言う関係にあることは有名
  • ゼータ関数と言うのは、(一般的に?)要素を1に対応づける関数のこと(?)
    • 少なくともポセットにおけるゼータ関数の議論は、ポセットが作る全てのインターバルに1を対応づける関数をゼータ関数と呼び、それに対応するメビウス関数のことなどを扱う

メビウスの反転公式

  • 自然数nを引数とする複素関数g(n)があるとする
  • 今、自然数nを引数とする複素関数f(n)をf(n) = \sum_{d|n} g(d)とする。ただし、(d|n)は「nがdで割り切れる」ことを意味する
  • このときメビウス関数\mu(n)を使ってg(n) = \sum_{d|n} f(d)\mu(n/d)となる
  • ここで注意しておきたいことは、f(n) =\sum_{d|n} g(d) = \sum_{d|n} g(d) \zeta(n|d); \zeta(m) = 1, m=1,2,...であること
  • 以下の2式は、Dirichlet畳み込みと呼ばれる
    • f(n) =(g * \zeta)(n)= \sum_{d|n} g(d) \zeta(n|d) = \sum_{i\times j = n} g(i) \zeta(j)
    • g(n) = (f * \mu)(n) = \sum_{d|n} f(d)\mu(n/d) = \sum_{i\times j = n} f(i) \mu(j)
  • なぜDirichlet畳み込みと呼ばれるかと言うと、(g*\zeta)(n)と言うnを引数とする関数のディリクレ級数が、f(n)のディリクレ級数\zeta(n)のディリクレ級数の積になり、また,(f*\mu)(n)と言うもう一つのnを引数とする関数のディリクレ級数f(n)のディリクレ級数\mu(n)のディリクレ級数との積になるため
  • ちなみに、ディリクレ級数の式の作りがn^sと関数本体[(tex:f(n)]など)の積になっているため、掛け合わせてnになるi,jの組み合わせに関する足し合わせ(\sum_{i\times j = n})が、ちょうどそのような対応関係を(nが1,2,...と無限に続く限り/全ての場合を尽くす限り)満足させるからである
  • ここにも書いていたが、メビウス関数の書き方が悪かった…
  • Rでやっておく


ディリクレ畳み込みとメビウス関係・ゼータ関数の関係

  • (f*g)(n) = \sum_{i \times j = n} f(i) g(j)と言う関数の構成方法をディリクレ畳み込みと言う
  • 先に見たf(n) = \sum_{d|n} g(n)g(n) = \sum_{d|n} f(d)\mu(n/d)の関係は、f = g * \zeta \Rightarrow g = f * \muと言う畳み込みの点で逆変換になっている

ポセット

  • 自然数にポセット構造を入れようとすると、単純に、自然数の大小関係により(N, \le)と言うものが思い浮かぶ
  • しかしながら、違う構造も入れられる。nがdで割り切れるとき、nはd以上である、と言うルールもポセットを作る。(N, |)。ただし[tex:|は割り切れる、と言う演算を表すものとする(この辺りから、素数、公約数の匂いがしてくる…)
  • ポセットでは、「インターバル」と言うものが対象になるが、グラフっぽくイメージしたい時には、自然数に自然に入る大小ポセット(N, \le)を考えよう。この時、1はインターバル、2もインターバル、[1,2]もインターバル、[4,5,6]もインターバル
  • 割り切れるかどうかで自然数にポセットを作れば、[1,2,4,8]はインターバル、[1,2,6]も[1,3,6]もインターバル
  • 連続なポセット構造で言うなら、分岐のある道があって、どの道も一方通行である時に、ある点から別の点に行き着けるとき、その歩き方のそれぞれがインターバル

余代数 coalgebraと接続代数

  • ポセットが持つ、インターバル全体の集合が作るベクトル空間はに余代数が入ると言う
  • 余代数は代数の逆さまみたいな作りのもののことだが、それについてはこちらに別途、メモをした
  • ポセットの余代数における、余積は\Delta([x,y]) := \sum_{z \in [x,y]} [x,z] \otimes [z,y]だと言う
  • そして、そこにおける余積演算をしても変えない元であるcounitは\delta([x,y]) :=1 \text{if } x=y; \text{otherwise } 0
  • 余代数は要素の余積(一つの要素をとって、それが指定する複数の要素の演算結果を返す)があるが、それに対応する代数とその代数の積(複数の要素の演算結果が一つの要素を返す)とに対応づくと言う特徴がある
  • この代数と余代数をつなぐのが、畳み込み演算
  • ちょっとこの辺りから怪しいのだが…
  • 余積の書式は畳み込み演算に見える
  • その際の畳み込みは、行列の積になる
  • 単なる自然数の大小関係ポセットの場合には、自然数列x自然数列の行列を作れば、値の入り方は三角行列的になる

Reduced incidence algebra

  • 接続代数に少し制約を入れる
  • [x,y][x',y']とがisomorphicな時に\phi([x,y]) = \phi([x',y'])であるような場合に限るという制約
  • こうすると[x,y][0,y-x]とに同じ値を与える関数を考えることができるようになって、全てのインターバルに対する値を、0起源の値とすることができて、自然数の数列として扱えるようになる

接続代数でのゼータ関数メビウス関数

  • 全てのインターバルに1を対応づけるのが接続代数でのゼータ関数
  • ディリクレ級数的にf = g * \zeta \Rightarrow g = f * \muを満足するのが、メビウス関数
  • f([x,y]) = \sum_{z \in [x,y]} g([x,z]) \zeta([z,y])=\sum_{z \in [x,y]} g([x,z]) ならば g([x,y]) = \sum_{z \in [x,y]} f([x,y])\mu([z,y])
  • 足し合わせ要素の関係が(d|n)と言うような割り切るか否か、になっていないので、対応するメビウス関数の方もそうはならず
  • \mu([x,y]) = 1; \text{if } x=y そうでなければ\mu([x,y]) = - \sum_{z _in [x,y]} \mu([x,z])と言うような\muに関する再帰的定義を持つ関数になる
  • これは集合のinclusion-exclusion的な足したり引いたり関係に相当する
  • さらに、reduced incidence algebraになって、全てのインターバルが、「インターバル距離」のみで決まるようになると
  • f(y) = \sum_{z \le y} g(z) \Rightarrow g(y) = \sum_{z \le y} f(z) \mu(k); k = y-zのようになる
  • このとき、\mu = 1 (k=0), \mu = -1 (k=1), \mu=0 (それ以外) のように簡単になる

リーマンのゼータ関数の時のメビウス関数と、接続代数のメビウス関数

  • 自然数を単に大小関係でポセット化するなら、reduced incidence algebraで得られたメビウス関数を用いる事になる
  • 自然数素因数分解してやると、畳み込みのおかげで、素数に関するメビウス関数値の積を取ることで、割り算型(d|n)に対応したメビウス関数が得られる。k=1のとき(1度だけ含まれる素数の場合)には-1を掛け、k=0(登場しない素数の場合)には、1を掛け、k>1(2度以上、登場する素数の場合)には、メビウス関数値を0にする、ことになり、これが両者の関係

行列の冪がばか大きくなる…

  • 行列をべき乗した時の要素の大小を相対的に検討したい
  • 単純にべき乗するとものすごく大きい数になって「計算できません」と言われてしまうことがあるので
  • 固有値分解して、固有値を対数計算することにする
  • 固有値の対数化をしてもまだ、ダメなので、固有値を一定の値で除しておいて、その除する値を別途つけて、相対的な値を持たせた冪業行列を返すことにする