2021-03-20

AHC001（マラソンマッチ）の参加者の使用言語の分布を調べてみた

競技プログラミングコンテストサイトである AtCoder のマラソン部門のコンテストに参加してみました： atcoder.jp 本記事では、コンテスト参加者の使用しているプログラミング言語はどういうものがあるのかを順位から調査してみました。

2021-02-27

ABC192 D - Base n を Python で解いてみる

ABC-D AtCoder python プログラミングコンテスト

問題概要

atcoder.jp

数字列 $X$ を $n$ 進法表記と見なしたとき, $M$ 以下になるような値が何種類あるか求めよ.

ただし, $n$ は $X$ に含まれる数字より大きいものとする.

考えたこと

入力例 1

22
10

22 を 3 進数と見なすと, $2\times 3 + 2 = 8 \leqq 10$
22 を 4 進数と見なすと, $2\times 4 + 2 = 10 \leqq 10$
22 を 5 進数と見なすと, $2\times 5 + 2 = 12 > 10$
22 を 6 進数と見なすと, $2\times 6 + 2 = 14 > 10$

により, 最初の2つだけが条件を満たします.

このことから, $n$ が小さいときには条件を満たし, $n$ を大きくするにつれて $M$ を超えてしまう閾値を見つければよさそうです.

以下, $X$ の長さが $K$ であるとし, $k$ 桁目を $a_k$ ( $a_k\leqq d$ )と表します. このとき, $X$ を $n(>d)$ 進法表記と見なすと,

[tex: \displaystyle X_n = \sum_{k=0}^{K-1} n^k a_k ]

と表すことができ, これは $X$ が二桁以上の場合には, $n$ について単調増大であり, したがっていつかは $M$ よりも大きくなります.

一方で, $X$ が一桁の場合には, 何進法表記であっても

$\displaystyle X_n = a_0$

です. いくら $n$ を大きくしても $M$ を超えないことがあり得ますが, 問題文で問われているのは $M$ 以下になるような値が何種類あるかであるため, $X$ そのものが $M$ 以下であれば $1$ が, そうでなければ $0$ が解答になります

解いていく

まずは, 入力を受け取ります. 数字列 $X$ は整数として入力されますが, $n$ 進数展開を行うことを見越してリストとして受け取りましょう.

X = [int(k) for k in list(input())]
d = max(X)
M = int(input())

また, 数字列 $X$ を $n$ 進数展開したときの値を p_adic(X, n) により計算します.

def p_adic(X, p):
  curr = 0
  for k in X:
    curr += k
    curr *= p
  curr //= p
  return curr

以下のようなコードだと, TLEしてしまいます

for ループが $O(M)$ の例 1

ans = 0
for j in range(d+1, M+1):
  if p_adic(X,j) <= M:
    ans += 1

print(ans)

for ループが $O(M)$ の例 2

ans = 0
for j in range(d+1, M+1):
  if p_adic(X,j) <= M:
    ans = j - d
  else:
    break

print(ans)

今回の制約では $M \leqq 10^{18}$ のため, より高速なアルゴリズムを構築する必要があります.

二分探索

今回の問題では $n$ が小さいときには条件を満たし, $n$ を大きくするにつれて, ある閾値を超えると条件を満たさなくなるという単調性があります.

したがって, その閾値を探す問題に帰着することができます. その際には,

は条件を満たすか？
- 満たさないならば閾値は $0$ から $M/2$ の間にある
- 満たすならば閾値は $M/2$ から $M$ の間にある

という性質を使って, 条件の範囲を1ステップごとに半分にしていくことができます. 効率よく探すことで, $O(\log M)$ の計算量で以上の閾値を求めていきましょう.

一般の二分探索については qiita.com

特に Python の二分探索をうまく書く方法については www.forcia.com

がそれぞれ読みやすく, 参考になります. ポイントとしては,

左端が True, 右端が False となるように初期値 ok, ngを設定する
中央値 (ok + ng)//2 での条件を判定する
- 左端が True, 右端が False となる状態を維持しながら ok, ngを変更する
ok と ng の差が1になったとき, ok の値が「条件を満たす $j$ の上限」である

ことになります.

解いていく

以上を踏まえて, - $X$ が一桁の場合がコーナーケースになっていることに注意する - 「 $X$ を $j$ 進数展開した結果が $M$ 以下である」をis_ok(j)と表記し, 二分探索で最大の $j$ を求める - 答えとして $j - d$ を出力する

を行うことにより, ACすることができます

X = [int(k) for k in list(input())]
d = max(X)
M = int(input())

def p_adic(X, p):
  curr = 0
  for k in X:
    curr += k
    curr *= p
  curr //= p
  return curr

def is_ok(j):
  curr = p_adic(X, j)
  if curr <= M:
    return True
  else:
    return False

if len(X) == 1:
  if X[0] <= M:
    print(1)
  else:
    print(0)
  import sys
  sys.exit()

ok = d
ng = M + 1
while ng - ok > 1:
  mid = (ok + ng) // 2 # 平均(小数切り捨て)
  if is_ok(mid):
    ok = mid
  else:
    ng = mid

print(ok-d)

atcoder.jp

2021-02-23

SIRモデル方程式を数値的に解いて、gifで表示させてみる

数値計算数学 python

伝染病の拡大過程を記述するSIRモデル
$\begin{equation} \left\{ \begin{split} S′(t)&=−βS(t)I(t),\\ I′(t)&=βS(t)I(t)−γI(t),\\ R′(t)&=γI(t) \end{split} \right. \end{equation}$
の初期値問題を Python で解き, 得られた解をgifにします.
f:id:wakabame:20210223221730g:plain

PCとインターネットがあればどなたでも, 以下のノートブックを逐次実行することでgifファイルを生成することができます.
colab.research.google.com

2020-10-05

水素原子の基底状態解は原点での値を用いずに特徴付けられる

偏微分方程式数学

3次元空間における水素原子ハミルトニアンに対する固有値問題 *1

\begin{equation}
-\dfrac{1}{2}\Delta \psi - \dfrac{1}{|x|} \psi = E \psi,\quad x \in \mathbb{R}^3
\end{equation}

について, 1s軌道に対応する基底状態解を考察します.

よく知られているように, 基底状態解は球対称であるため,
$u(|x|) = \psi(x)$ と変換することにより, 未知関数 $u: [0,+\infty)\rightarrow \mathbb{C}$ についての方程式

\begin{equation}
\left(-\dfrac{1}{2}\dfrac{d^2}{dr^2} -\dfrac{1}{r}\dfrac{d}{dr} -\dfrac{1}{r}\right)u = E u,\quad r \in (0,\infty) \tag{1}
\end{equation}

と変形できます.
さらに, (1) に対して, $r^2 \bar{\phi}$ をかけて部分積分することにより弱形式

\begin{equation}
\left\langle \dfrac{r}{2}\dfrac{du}{dr}, \dfrac{r}{2}\dfrac{d\phi}{dr}\right\rangle
- \left\langle \sqrt{r}u, \sqrt{r}\phi \right\rangle
=E \left\langle ru, r\phi\right\rangle \tag{2}
\end{equation}

を得ます*2.

本記事では, 1次元における固有値問題 (1) に焦点を当てて, 適当な関数空間での変分問題と見なすことで, 特に陽に境界条件を与えることなく基底状態解が得られることを解説していきます.

また, 本記事は有限要素法による数値計算について紹介している, dc1394さんの記事
qiita.com
についての考察を兼ねています.
モデリングや数値計算の箇所を読んでいて, 境界条件を設定していないように見えるのに, 正しく基底状態解が得られていることが気になり調べてみました.
なお, 弱形式の導出などについても紹介記事にて丁寧に行われているため, 気になる方はそちらを参考にしてみてください.

*1:単に, シュレディンガー方程式とも呼びます

*2:部分積分から境界項が出ますが, その項が 0 になるよう $\bar{\phi}$ が良い性質を持っていることを仮定しますが, ここでは細かいことは気にしない.

2020-09-22

変分法の直接法で、ラプラシアンの第一固有関数と第一固有値を手に入れる

偏微分方程式数学

この記事では, 滑らかな有界領域 $\Omega \subset \mathbb{R}^n$ に対して Dirichlet Laplacian の第一固有関数 $u_1$ と第一固有値 $\lambda_1$ を変分的な手法で得ます.
すなわち, 次の Dirichlet 条件付き楕円型偏微分方程式を満たす固有値 $\lambda\in \mathbb{R}$ と固有関数 $u=u(x)$ を求める問題を考えましょう.

\begin{equation}
\begin{cases}
- \Delta u = \lambda u, &\text{in } \Omega,\\
u = 0, &\text{on } \partial \Omega.
\end{cases}\tag{P}
\end{equation}

簡単のため, 固有値の中でも最小の第一固有値と, 第一固有関数に着目して解いていきます.

レイリー・リッツの特徴づけ

第一固有値 $\lambda_1$ は次のように特徴づけられる.
\begin{equation}
\lambda_1= \inf\left\{ \frac{\|\nabla v\|^2_2}{\| v\|^2_2} ;v\in H_0^1(\Omega)\setminus\{0\}\right\} \tag{M}
\end{equation}

ただし, $H_0^1(\Omega)$ はソボレフ空間, $\|\cdot\|_2$ は $L^2$ ノルムを表します.
以下では, まずはこの最小化問題が解を持つことを解説していき, 次いでその解がもとの固有値問題の解になっていることを述べます.

2020-09-02

最適輸送「モンジュの問題」と「カントロビッチの問題」

最適輸送理論

最適輸送問題とは？

ある砂山の形 $\mu$ を別の砂山の形 $\nu$ に移すことを考える.
砂山の移動に掛かる労力が「質量」と「始点と終点の位置」で決まっているとき, トータルの労力が最小になるような移し方はどんなものだろうか？

数学的に定式化すると

砂がありえる場所の全体を領域 $\Omega \subset \mathbb{R}^d$ として表し, 砂山の形をそれぞれ確率測度 $\mu, \nu \in P(\Omega)$ で表し, 関数 $c:\Omega\times\Omega \rightarrow [0,+\infty$ ] として単位体積を始点から終点へ輸送したときの労力を表す.

モンジュの問題

次の最小化問題を満たす $T$ を見つけ出せ;
\begin{equation}
\inf\left\{\int c(x,T(x))\,d\mu(x)\,|\, T:\Omega\rightarrow\Omega :\text{可測}, T_{\#}\mu =\nu\right\}.
\end{equation}
ここで, $T_{\#}\mu$ は
\begin{equation}
T_{\#}\mu(A) = \mu(T^{-1}(A)),\quad A\in \mathcal{B}(\Omega)
\end{equation}
で与えられる測度であり, $\mu$ の $T$ による押し出し (push forward) という.

また, この最小化問題の解 $T$ は ( $\mu$ から $\nu$ への) 最適輸送写像と呼ばれる.

カントロビッチの問題

次の最小化問題を満たす $\gamma$ を見つけ出せ;
\begin{equation}
\inf\left\{\int\int_{\Omega\times\Omega} c\,d\gamma\,|\, \gamma\in\Pi(\mu,\nu)\right\}.
\end{equation}
ここで, $\Pi(\mu,\nu)$ は
\begin{equation}
\Pi(\mu,\nu) = \left\{ \gamma\in P(\Omega\times\Omega)\,|\, (\pi_1)_{\#}\gamma = \mu, (\pi_2)_{\#}\gamma = \nu\right\}
\end{equation}
により与えられる測度の族であり, この要素を ( $\mu$ から $\nu$ への) 輸送計画と呼ぶ.
ただしここで, $\pi_1, \pi_2$ はそれぞれ $\pi_1(A\times B) = A, \pi_2(A\times B) = B$ なる射影である.

また, この最小化問題の解 $\gamma$ は ( $\mu$ から $\nu$ への) 最適輸送計画と呼ばれる.

カントロビッチの問題は割と解ける

関数 $c :X\times Y\rightarrow \mathbb{R}\cup \{+\infty\}$ が下に有界かつ下半連続であるとする.
このとき, カントロビッチの問題は解を持つ.

モンジュの問題が解けるならば, カントロビッチの問題の最適輸送計画はそれ以下のコストを持つ

輸送写像 $T$ がモンジュの問題の解とすると,
\begin{equation}
\gamma_T = (\text{id}, T)_{\#}\mu
\end{equation}
は $\gamma_T \in \Pi(\mu,\nu)$ により,
\begin{align}
\int c(x,T(x))\,d\mu(x)
&= \int\int_{\Omega\times\Omega} c(x,y)\,d\gamma_T(x,y)\\
&\ge \inf\left\{\int\int_{\Omega\times\Omega} c\,d\gamma\,|\, \gamma\in\Pi(\mu,\nu)\right\}.
\end{align}

モンジュの問題は解けないが, カントロビッチの問題は解ける場合がある

これを見るためには,
\begin{align}
\mu &= \delta_0,\\
\nu &= \dfrac{\delta_0 + \delta_1}{2},\\
c(x,y)&= |x-y|
\end{align}
の条件で, それぞれの問題を考察すればよい.
つまり, $0$ にあった砂のうちの半分を $1$ に輸送することを考える.
このとき, モンジュの問題は $T(0)$ の値がいずれであっても輸送写像が作れず, 一方でカントロビッチの問題は
\begin{equation}
\gamma(\{0\}\times \{0\}) = \gamma(\{0\}\times\{1\}) = 1/2
\end{equation}
となるようなものを取れば, これが最適輸送計画となる.

2020-08-10

無限次元ヒルベルト空間の有界点列は収束部分列を持たないが、弱収束部分列なら持つ

数学関数解析

この記事では, 無限次元ヒルベルト空間 $H$ について,

有界な点列であっても収束するような部分列がない場合があること
弱収束を考えれば収束するような部分列があること

を述べていきます.