１章 駆け足で読む『トポロジカル・インデックス』その２

トポロジカル・インデックス: フィボナッチ数からピタゴラスの三角形までをつなぐ新しい数学

• フィボナッチ数列
  • 

# 最初の２項 f0, f1 を指定
# 漸化式の係数a,bを指定
fn <- function(n, f0 = 1, f1 = 1,a = 1, b = 1){
 ret <- rep(0,n+1)
 ret[1] <- f0; ret[2] <- f1;
 if(n > 1){
   for(i in 2:n){
    ret[i+1] <- a * ret[i-1] + b * ret[i]
  }
 }
 ret
}
fn(20)

> fn(20)
 [1]     1     1     2     3     5     8    13    21    34
[10]    55    89   144   233   377   610   987  1597  2584
[19]  4181  6765 10946

• • 負の方向にも定義できる

• • • フィボナッチ数列の値が0のところを中心に正負で絶対値的に対称。ただし正方向はすべて正、負方向は符号の交代がおきる

inv.fn <- function(n, f0 = 1, f1 = 1,a = 1, b = 1){
	ret <- seq(0,n+2)
	ret[1] <- f1; ret[2] <- f0;
	if(n < 0){
		for(i in 1:(-n)){
			ret[i+2] <- (ret[i] - a * ret[i+1])/b
		}
	}
	ret
}
inv.fn(-2)

both.fn <- function(n, f0 = 1, f1 = 1, a =1, b = 1){
	# 正負の対称をとるために、fn(n-2),inv.fn(-n)としている
	# 値が0であるところが中心で絶対値的に左右対称である
	# その意味では i = 0 のときに0とするのが本当はよい
	posi <- fn(n-2,f0=f0,f1=f1,a=a,b=b)
	nega <- inv.fn(-n,f0=f0,f1=f1,a=a,b=b)
	
	c(nega[(n+2):3],posi)
}

both.fn(20)

plot(both.fn(30),type="b")

• • フィボナッチ数、ルカ数、ペル数、ペル-ルカ数は、漸化式の係数と初期値のバリエーション

# フィボナッチ
n <- 10
both.fn(n,1,1,1,1) -> ff
# ルカ
both.fn(n,2,1,1,1) -> ll
# ペル
both.fn(n,1,2,1,2) -> pp
# ペル-ルカ
both.fn(n,2,2,1,2) -> pl

matplot(cbind(ff,ll,pp,pl),type="l")

• 漸化式の行列表現
  • のとき
  • 
  • 行列の２行目は(1,0)。これは、に対応する

library(expm)
fn.2 <- function(ns, f0 = 1, f1 = 1, a = 1, b =1){
 m <- matrix(c(a,b,1,0),byrow=TRUE,2,2)
 ret <- matrix(0,length(ns),2)
 for(i in 1:length(ns)){
  ret[i,] <- m%^%ns[i] %*% c(1,0)
}
 ret
}
fn.2(0:30)

• • 行列のべき乗を固有値に意識してRで書くと
• 行列のべき乗と指数行列

# べき乗
pow.m <- function(A,n){
	# 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(eigen.out[[1]]^n)
	X <- V%*%B%*%U
	return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
}

# 指数行列
exp.m <- function(A,n){
	# 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(exp(eigen.out[[1]]*n))
	X <- V%*%B%*%U
	return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
}
pow.m(matrix(c(1,1,1,0),byrow=TRUE,2,2),2)

> pow.m(matrix(c(1,1,1,0),byrow=TRUE,2,2),2)
$matrix
     [,1] [,2]
[1,]    2    1
[2,]    1    1

[[2]]
[1]  1.618034 -0.618034

• 連分数
  • 

renbunsuu <- function(as){
	ret <- as[1]
	tmp <- 1/as[length(as)]
	if(length(as) > 2){
		for(i in (length(as)-1):2){
			tmp <- 1/(as[i] + tmp)
		}

	}
	ret + tmp
}

renbunsuu(c(1,1,1))

# フィボナッチ数の比は1の連分数
ff <- fn(10)

for(i in 1:(length(ff)-1)){
	print(ff[i+1]/ff[i])
	print(renbunsuu(rep(1,i)))
}
# ペル数の比は2の連分数
pp <- fn(10,1,2,1,2)

for(i in 1:(length(pp)-1)){
	print(pp[i+1]/pp[i])
	print(renbunsuu(rep(2,i)))
}

• • 連分多項式を使うと、連分多項式を行列に対応付けて、その行列式として計算できる

mat.renbun.takou <- function(xs){
	if(length(xs)==1){
		m <- matrix(xs,1,1)
	}else{
		m <- diag(xs)
		for(i in 2:length(xs)){
			m[i-1,i] <- 1
			m[i,i-1] <- -1
		}
	}
	m
}

renbun.takou <- function(xs){
	m <- mat.renbun.takou(xs)
	return(list(m = m, v = det(m)))
}


n <-10
# フィボナッチ
for(i in 1:n){
	print(renbun.takou(rep(1,i))$v)
}
# ルカ
for(i in 1:n){
	if(i == 1){
		print(renbun.takou(c(1))$v)
	}else if(i == 2){
		print(renbun.takou(c(1,2))$v)
	}else{
		print(renbun.takou(c(1,2,rep(1,i-2)))$v)
	}
	
}

# ペレ
for(i in 1:n){
	print(renbun.takou(rep(2,i))$v)	
}
# ペレ-ルカ
for(i in 1:n){
	if(i == 1){
		print(2*renbun.takou(c(1))$v)
	}else{
		print(2*renbun.takou(c(1,rep(2,i-1)))$v)
	}
	
}

• 数列を係数とする多項式、チェビシェフ多項式