京都大学(一般入試)数学の問題を計算機で解く(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