ゆるふわブログ

プログラミング,大学の勉強,日常生活で感じたことをゆるふわに書いていけたらなと思います.技術的に拙いところがあっても温かい目で見守っていただければ幸いです.

包除原理が好き

こんにちは.本日は包除原理の入門を書いてみようと思います. 競技プログラミングをやられている方には自明な内容だと思います.

本日の目標

目標としては以下が当たり前だと感じられるようになることです.

全射の個数
素数 m であるような有限集合  A から要素数 n であるような有限集合  B への全射  f: A \rightarrow B の個数は,  \begin{align*}
\sum_{i=0}^n (-1)^{n-i} \binom{n}{i} i^m
\end{align*}

 

包除原理について

有限集合  A の要素数 |A| とあらわすことにしましょう.高等学校で以下のような公式を習ったことと思います.

 |A \cup B| = |A|+|B|-|A \cap B|  |A \cup B \cup C| = |A|+|B|+|C|-|A \cap B|-|B \cap C|-|C \cap A|+|A \cap B \cap C|

これは集合の合併の要素数を集合の 1 個以上の交叉の要素数で記述する公式で,まさにこれが包除原理の特殊ケースです.合併をとる集合の個数を一般化して,以下が成り立ちます.

包除原理
 T = \{1,2,\ldots,n\} とすると,  \begin{align*}
|\cup_{i=1}^n A_i| = \sum_{S \subseteq T} (-1)^{|S|-1} |\cap_{i \in S} A_i|
\end{align*}

証明は, n=3 の場合にベン図を書いて,disjoint な領域に分け,それぞれが数え上げられている回数を考えれば良いです.一般には,

 \begin{align*}
(A_{i_1} \cap \cdots \cap A_{i_k}) \setminus (\cup_{j \in T \setminus \{i_1, \ldots, i_k\} } A_j) \quad (i_1 < \cdots < i_k)
\end{align*}

が右辺において何回数え上げられているかを考えます.これを含むのは  {i_1, \ldots, i_k} の部分集合に対応する交叉だけですから,二項定理より,

 \begin{align*}
\sum_{m=1}^k (-1)^{m-1} \binom{k}{m} = 1-\sum_{m=0}^k (-1)^{m} \binom{k}{m} = 1-(1-1)^k = 1
\end{align*}

です.

参考として以下を上げておきます.

mathtrain.jp

余事象バージョン

さて,競技プログラミングではいくつかの集合 (条件) の交叉の数え上げを行うことがよくあります.ド・モルガンの法則により,以下の余事象バージョンが成り立つからです.

包除原理 (余事象バージョン)
 T = \{1,2,\ldots,n\} とすると,  \begin{align*}
|\cap_{i=1}^n A_i| = \sum_{S \subseteq T} (-1)^{|S|} |\cap_{i \in S} A_i^c|
\end{align*} ただし, A^c は集合  A の補集合を表すとする.(全体集合が暗に仮定されている.空の添え字集合に対する  \cap は全体集合を表すものとする.)

故に,「いくつかの条件が「かつ」でつながっている時の数え上げは,それぞれの条件の余事象の「かつ」が楽に計算できるとうれしい.」 となります.

また,添え字集合の対称性を利用する場合もあります.愚直な包除原理は部分集合に関する和を取っているため,計算時間は exponential になりがちです.しかしながら,対称性があるならば,「いくつの集合の交叉であるか」のみが本質であるため,二項係数で同じ数の交叉をまとめることができ,線形時間で計算が可能です.( \sum_{i=0}^n \binom{n}{i} = 2^n は二項定理でも二進数の popcount を考えても明らかですが,それと要素数  n の部分集合の個数が  2^n であることが整合します.)

全射の数

さて,全射とは, B の各要素の逆像が空でない,すなわち各  B の要素についてそれへ対応づけられている  A の要素が少なくとも 1 つ存在することを言います.少なくとも 1 つと見たら余事象です.これは  A_i を 「 i 番目の  B の要素に少なくとも 1 つ対応づけられているものが存在する」として,「かつ」の条件で書けます.故に, A_i^c すなわち「 i 番目の要素に対応づけられるものが 1 つもない」がいくつか「かつ」でつながったものを数えられれば良いです.これは簡単で, n-i 個に絶対入らないようにしたいなら,各  A の要素について  i 個の行き先の選択肢しかないものとすればよく, i^m で求められます.(ここで注意したいのは条件がかかっているのは「(固定された)  n-i 個に絶対入らない」だけであって,他には何も課されていないことです.故に 1 つも対応づけられていない要素が  n-i 個以上であっても良いのです.(ここに慣れるのが重要) ) 対称性を利用して二項係数を用いて余事象バージョンの包除原理を考えれば,冒頭の式は明らかです.

その他

その他,競技プログラミングでは,動的計画法 (dynamic programming; dp) と組み合わせて数え上げる問題,高速メビウス変換を用いる問題,( (-1)^n より偶奇のみが重要であるため) 偶奇のみに着目する場合等があるようです.完全順列 (撹乱順列) も包除原理の典型的な応用例として有名です.確か包除原理は集合の包含関係から定まる半順序において,拡張されたメビウス関数についてのメビウスの反転公式から導かれるという話があったと思います.それについては後日余裕があれば書きます.

先日の ABC (AtCoder Beginner Contest) で包除原理の問題が出ました.

atcoder.jp

また,僕の記憶が正しければ以下はドンピシャで全射の数の話です.

codeforces.com

今年の年末

年末はずっと家に引きこもって寝るか卒論を執筆するかしている.12 月初旬に中間発表を終え,特に新規性も出ないままはてどうするかととまどっていたのであるが,同じ分野をやっていた研究室の先輩の卒論をじっくり読み返していたら奇跡的にその拡張を思いついた.個人的にはそこまで非自明な結果ではないが,あまりにもきれいな結果であるので,他分野のどこかで誰かが同等のことを思いついているであろう.私がやっている分野では知る限り誰も思いついていない.恋人がベットで寝ている夜明けに夢中で数式を書きなぐった.思いついたものをすぐに実装し,いくつかの観点で「自分の出した結果は確かか?」を検証した.するとどうやら間違いないようであり,理論的に最適なものが導けたようである.確信を得てから指導教員に報告に行った.

そこからはずっと夢中だった.うつ気味でやる気のしない時もあったが,シャーペンを持って紙に向かっている時,プログラムを実装している時,論文を書いている時,夢中になれた.

先日卒論の第 1 原稿を持って指導教員とミーティングをしたのであるが,先生は「数学の能力も高そうであるし,研究者になる道も考えてみてはどうか」と言ってくださった.うれしかった.無論,自分に特筆するような才能はないことはわかっているし,先生は新たな研究者を輩出することを目的とした教育者であるからそのようなことをおっしゃったのだろうが,それでもなんだか無性に喜ばしかった.

卒論を学術論文として外部に出してはどうかと勧められている.卒論も英語で執筆することにした.いろいろなことに挑戦できている気がする.来年からは同大学院の修士に進学する.どこまでやれるかわからないが,希望が少し見えてきた.

立体切断を多角的に捉える

こんにちは、Ysmr-Ry です。

先日このようなツイートを見かけました。

このような立体切断の問題、小学校の頃何もわからず解いた記憶がありますが、こういう発想をするには 空間認知能力がなければいけないのでしょうか?

そこで今回は空間認識能力皆無な私がこの立体の切断について理解を深める過程を見ていきたいと思います。

まず思うこと

まず思うのは有名な立体切断です。立方体の重心を通り、各辺の中点を通るように切断すると正六角形が現れます。
このツイートは、切り口を六角形の6方向に伸ばすべく、立方体を各面に貼り合わせて計 7 個にしたものです。
立体切断の際は、軸を伸ばして三角形を描くように言われるのですが、まさにここではその三角形が切り口としてあらわになるようにうまく立体を足した形になっています。

このことは、Visualizer を作ったのでイメージのわかない方はそちらで確かめてみてください。初期位置が六芒星の切断となっているので、矢印キーでカメラをぐりぐり動かしてみてください。X, C キーは押さないほうが良いです。(サイトに飛んだ方がいいです)

codepen.io
CodePen - HexagonCube

空間認知能力がないならば

ないならば方程式で処理してしまえばいいんですね〜 代数的にゴリゴリするだけなら誰でもできます。やってみましょう。

まず、3 次元座標を設定し、その中心が考える立体の "中心" の立方体の重心になるようにします。
この立体を不等式で表してみましょう。

ここで知っていると見通しの良いのは、正方形が一種の "円" であるということです。
通常の意味での閉円板、すなわち円の周及び内部は中心  a、半径  r として、 D(a,r) = \{x | d_2(a,x) \leq r\} で表されます。
ただし、 d_2( (x_1,y_1),(x_2,y_2) ) = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2} は Euclid 距離です。この距離の部分を別の距離に置き換えることで、円板を拡張できます。閉円板を正方形にするには Chebychev 距離  d_\infty( (x_1,y_1),(x_2,y_2) ) = \max\{|x_1-x_2|,|y_1-y_2|\} を使います。

考える立体を 3 つに分けましょう。すなわち、2x2x6 の直方体が x, y, z 軸にそうように重なっているとします。(軸の方向はそうとるものとし、1 つの立方体の寸法は 2x2x2 とする) xy 平面上で 2x2 の正方形であり、z 軸方向へは無限に伸びる四角柱は、先の議論により、 d_\infty(O,(x,y) )=max\{|x|,|y|\} \leq 1 と表されます。これに、z 軸方向が [-3,3] の範囲であることを組み込んで、 max\{|x|,|y|,|\frac{z}{3}|\} \leq 1 となります。他の軸については、x, y, z を cyclic に置換して行けば良いです。すなわち、

 max\{|x|,|y|,|\frac{z}{3}|\} \leq 1 \cup max\{|y|,|z|,|\frac{x}{3}|\} \leq 1 \cup max\{|z|,|x|,|\frac{y}{3}|\} \leq 1

が求める立体です。

さて、2 次元のグラフ描画ツールなら desmos という便利なのを知っているのですが、私はなにぶん Geogebra とかできないので 3 次元の描画はめんどい (というかわかりづらい) のです。そこで、3 次元なものを扱うときの重要な観点なのですが、平面で切断した切り口など 2 次元の情報に落とし込んで考察する ことを考えます。いやまあ、立体切断なんだから当たり前なんですが。先に正六角形の切断がベースになっていると述べました。この切断面は定式化すると、原点を通り、法線ベクトル (平面と直交している向き) が (1,1,1) であるような平面  x+y+z=0 です。これを平面の向きはそのままに動かしていくので、パラメータ  k を使って  x+y+z=k としておきましょう。

切断平面上に 2 次元座標 (s,t) を設定し、s, t の方程式として切り口を立式できれば、グラフツールで st 平面上にグラフを書けば切り口を Visualize できそうです。その際、切断平面上の正規直交基底 (長さ 1 で互いに直交する "ものさし") を (1,1,1) に直交する 2 ベクトルから選び出す必要があります。ここでは、 x+y+z=0 の任意の点  P に対し、

 \overrightarrow{OP} = s\begin{pmatrix} \frac{1}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}} \\ 0 \end{pmatrix} + t\begin{pmatrix} \frac{1}{\sqrt{6}} \\ \frac{1}{\sqrt{6}} \\ -\frac{\sqrt{2}}{\sqrt{3}} \end{pmatrix}

としましょう。ここで、一般には  x+y+z=k ですから、座標原点の変位ベクトルを (k/3,k/3,k/3) として足しあわせましょう。すると、

 \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} \frac{s}{\sqrt{2}}+\frac{t}{\sqrt{6}}+\frac{k}{3} \\ -\frac{s}{\sqrt{2}}+\frac{t}{\sqrt{6}}+\frac{k}{3} \\ -\sqrt{\frac{2}{3}}t+\frac{k}{3} \end{pmatrix}

となります。これを求める立体の不等式に代入すれば完成です。cyclic な 3 つの直方体のうち、一つを示せば十分でしょう。

 \max\{|x|,|y|,|\frac{z}{3}|\} \leq 1 \iff \max\left\{\left|\frac{s}{\sqrt{2}}+\frac{t}{\sqrt{6}}+\frac{k}{3}\right|,\left|-\frac{s}{\sqrt{2}}+\frac{t}{\sqrt{6}}+\frac{k}{3}\right|,\left|\frac{-\sqrt{\frac{2}{3}}t+\frac{k}{3}}{3}\right|\right\}\le1

 \max\{a,b,c\} \leq 1 \iff a \leq 1 \land b \leq 1 \land c \leq 1 であることに注意すると、これは幾つかの線形な不等式群のなす領域であることがわかります。(すなわち、幾つかの直線に囲まれた領域。線形計画法の実行可能領域)

このような領域を図示したものが以下になります。(ページに飛んでもらった方がいじりやすいです。

www.desmos.com

色分けとしては z, x, y 軸にそう直方体がそれぞれ紫、青、緑であり (それぞれの領域を  S, T, U とする)、その共通部分が橙です ( V とする)。橙の領域は立体を "中心" の立方体に限った時の切り口です。(ゆえに  k=0 の時は正六角形になっている)

左側にある  k の値をスライドさせていじると領域が動きます。 k=0 の時が六芒星 k=1 の時は、"中心" の立方体において 3 面の対角線のなす正三角形が切り口となっている状況、 k=3 の時は "中心" の立方体の (1,1,1) 方向の頂点をギリギリかすめている位置にあり、ここから先は橙の領域がなくなります。(ここら辺は冒頭の 3d Visualizer と合わせて確認してみてください)

さて、 k をぐりぐり動かしているといろいろと見えてきます。 S, T, U は基本的に平行四辺形 (正確にはひし形) であり、ある t 軸に垂直な水平線が上下にあって (図中の一番上 / 下の赤い線です)、そこを超えてしまうと削られていきます。また、 S, T, U は原点を中心として 120° ずつ回転させた関係になっています。

求積

さて、数学 III あたりを履修したことのある人は、立体の切り口というと積分を思い浮かべると思います。この立体は 2x2x2 の立方体が 7 個なので体積は簡単に求まるのですが、切り口の面積を解析的に計算し積分することで等しくなることを確かめましょう。

面積を求める上で少し違う視点に立つとやりやすいです。すなわち、切り口を  V\sqcup(S-V)\sqcup(T-V)\sqcup(U-V) と捉えてみましょう。すると、先ほどの水平線の間にさらに 2 本の水平線があり (図中の中央 2 本の赤い線です)、橙の部分はそれに切り取られていることがわかります。 S-V, T-V, U-V の面積は等しいので、 V, S-V の面積が求まれば  k における切り口の面積  S(k) が求まります。平行四辺形が水平線で切り取られている部分は相似を使うと容易に求めることができ(相似比の 2 乗ってやつです)、まとめると、

 S(k) = \begin{cases}
6\sqrt{3}\left(\left(\frac{k}{2}+\frac{1}{2}\right)^2+\left(-\frac{k}{2}+\frac{1}{2}\right)^2\right)+2\sqrt{3}\left(2-\left(\frac{1}{2}+\frac{k^2}{2}\right)\right)\ \quad (0\le k\le1) \\
6\sqrt{3}\left(2-\left(\frac{k}{2}-\frac{3}{2}\right)^2-\left(\frac{k}{2}-\frac{1}{2}\right)^2\right)+2\sqrt{3}\left(\frac{k}{2}-\frac{3}{2}\right)^2\ \quad (1\le k\le3) \\
6\sqrt{3}\left(\frac{5}{2}-\frac{k}{2}\right)^2\ \quad (3\le k\le5)
\end{cases}

と、滑らかに接続された放物線群となります。( k\geq 0 の時) これを積分します。
この際注意しなければならないのは、 k の 1 目盛が実際の長さ 1 ではないということです。変位ベクトルの長さは  |(\frac{1}{3},\frac{1}{3},\frac{1}{3})|=\frac{1}{\sqrt{3}} であるので、これをかけねばなりません。さらに、 k が正の時と負の時では同じ挙動をしめす ( S(k) は偶関数であり、元ツイートの右側の立体は全立体のちょうど半分) なので、体積は

 2\int_0^5 S(k) \frac{dk}{\sqrt{3}} = 56 = 2^3\cdot 7

となり、一致しました。

お疲れ様でした。

終わりに

冒頭の 3d Visualizer で X/C キーで k を増減できるので試してみてください!

他にも切断平面の法線ベクトルを自由に変えたらどうなるかとか考えているので続編を乞うご期待!

AtCoder Beginner Contest 106 敗戦記 (Haskell/C++)

はい…

A Garden

やるだけ

import Control.Applicative ((<$>))

main :: IO ()
main = do
  [a,b] <- map read . words <$> getLine :: IO [Int]
  print $ a*b-a-b+1
B 105

やるだけ。約数列挙 O(N) でも間に合いますね。

main :: IO ()
main = do
  n <- readLn :: IO Int
  print $ length $ filter odd8 [1..n]

factors :: Int -> [Int]
factors n = filter ((==0).(n`mod`)) [1..n]

odd8 :: Int -> Bool
odd8 n
  | even n    = False
  | otherwise = length (factors n) == 8
C To Infinity

1 でないなら 5000 兆乗するので指数関数的に爆発していきます。
よって、先頭から見ていって最初に現れる 1 でないものが基本的には答えです。
k 文字 1 が続くなら 1 です。(ここら辺のコーナーケースで何回か失敗した記憶) (水色失格)

import Data.List
import Data.Maybe

main :: IO ()
main = do
  s <- getLine
  k <- readLn :: IO Integer
  let none = filter (/='1') s
      noneId = findIndex (/='1') s

  putChar $ if null none
    then '1'
    else if fromJust (noneId) < fromIntegral k 
      then head none
      else '1'

  putStrLn ""
D AtCoder Express 2

できませんでした()
まあ結論をいうと区間を 2 次元座標と捉えて 2 次元累積和をします (この発想がなかった)
[p, q] のクエリがきた時は座標が (p,p), (p,p+1), ..., (p, q), (p+1, p), (p+1, p+1), ... (q,q) のところ ([(p,p),(q,q)] の矩形) の和がわかれば良いです。

#include <cstdio>
#include <algorithm>
#include <vector>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

int N, M, Q;
int I[510][510], S[510][510];

int main()
{
  scanf( "%d%d%d", &N, &M, &Q );
  rep( i, M )
  {
    int L, R;
    scanf( "%d%d", &L, &R );

    ++S[L][R];
  }

  rep( i, N+1 ) rep( j, N+1 )
    S[i][j+1] += S[i][j];

  rep( j, N+1 ) rep( i, N+1 )
    S[i+1][j] += S[i][j];

  rep( t, Q )
  {
    int p, q;
    scanf( "%d%d", &p, &q );

    printf( "%d\n", S[q][q]-S[q][p-1]-S[p-1][q]+S[p-1][p-1] );
  }
  
  return 0;
}
おわりに

言い訳を述べますが、競技プログラミングにおいてこういう瞬間はあるんだと思います。
自分の中にそのアイデアというか概念がなければ、それが即わからない・解けないにつながってしまうのです。
〜色といっても(知らんですが)〜点問題を〜%の正答率で解けるみたいな感じだと思うので、それを 100% に近づけるのは点数低めであっても難しいものなのだろうと。
"自分がいかにわかっていないか" を思い知る瞬間というのは何を勉強していてもあると思うし、それが非常に重要であると思うので、こういう瞬間に巡り会う機会を増やす (コンテストになるべく出る)、巡り合った時に確実に潰す、この 2 つを大事にしていきたいと思いました。
水色といってもまだまだ ABC を確実に全完できるレベルには達していないことが今回わかったし、初心者を早く脱っさなければならないとお尻に火がつきました。
これからも ABC といっても油断せず参加したいと思います。

AtCoder Beginner Contest 105 (by Haskell)

はい。

A AtCoder Crackers

問題よく読んでなかったですが 0 か 1 しかない感。

import Control.Applicative ((<$>))

fi = fromIntegral

main :: IO ()
main = do
  [n,k] <- map read . words <$> getLine :: IO [Int]
  print $ if min (ceiling (fi n/fi k)*k-n) (n-(n`div`k)*k) > 0
    then 1
    else 0

B Cakes and Donuts

ループ回して探索すれば良いと思います。

main :: IO ()
main = do
  n <- readLn :: IO Int
  putStrLn $ if null $ filter (\i -> i*4 <= n && (n-i*4)`mod`7 == 0) [0..n]
    then "No"
    else "Yes"

C Base -2 Number

なんかよく考えたら 2 進数もどきなのでどうとでもなった気がするが、頭が死んでいたので愚直にやった。divMod と quotRem 的な罠にはまった() anamorphism です。

import Data.List
import Data.Char

div' :: Int -> Int -> Int
div' a b
  | a`mod`b == 0 = a`div`b
  | otherwise    = a`div`b+1
  
mod' :: Int -> Int -> Int
mod' a b
  | a`mod`b==0 = 0
  | otherwise  = a`mod`b-b

f :: Int -> Maybe (Char, Int)
f i
  | i == 0    = Nothing
  | i == 1    = Just ('1',0)
  | otherwise = Just (chr (ord '0'+(i`mod'`(-2))), i`div'`(-2))

main :: IO ()
main = do
  n <- readLn :: IO Int
  let cs = unfoldr f n

  putStrLn $ if null cs
    then "0"
    else foldl' (\ans c -> c : ans)  "" cs

D Candy Distribution

一瞬アレっとなりました。C++ で AC してから Haskell で書き直した…
まあ累積和を取るところ、mod M で考えるところまではパッと出てくると思います。
制約をみると O(N^2) で探索させる気はないので、線形オーダーでどうにかできればうれしい。
区間 [i,j) の和が S[j]-S[i] で表せることを思えば、S[i] を mod M の値で分類し、頻度を数えれば良いことがわかります。(S[i] と S[j] の mod M が同じ時その差は M の倍数。mod M が同じものから異なる 2 つを選ぶ場合の数)
mod M が 0 の場合は自身だけを選ぶ場合も含まれることに注意。

import Control.Applicative ((<$>))
import Data.List

main :: IO ()
main = do
  [n,m] <- map read . words <$> getLine :: IO [Integer]
  as <- map read . words <$> getLine :: IO [Integer]
  let
    ss = map (`mod` m) $ scanl1 (+) as
    cs = map (\i -> (head i,genericLength i)) $ group $ sort ss

  print $ sum $ map (\(i,c) -> c*(c-1)`div`2 + (if i == 0 then c else 0)) cs

なんとなく行列を思い浮かべて、対角成分と上三角成分のようなイメージが湧いた()
直交行列の群 O(n) が Cayley 変換により n(n-1)/2 次元の Lie 群とみなせることを思い出して終了しました((((

Mujin Programming Contest

直した

A コンテスト名

やるだけ

#include <iostream>
#include <string>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

std::string S;

int main()
{
  std::cin >> S;
  
  std::cout << (S.substr(0,5) == "MUJIN" ? "Yes" : "No") << std::endl;
  
  return 0;
}
B セキュリティ

やるだけ

#include <iostream>
#include <string>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

int A;
std::string S;

int main()
{
  std::cin >> A >> S;

  if( A == 0 )
  {
    puts("Yes");

    return 0;
  }

  rep( i, S.size() )
  {
    if( S[i] == '+' )
      ++A;
    else
      --A;

    if( A == 0 )
    {
      puts("Yes");

      return 0;
    }
  }

  puts("No");
  
  return 0;
}
C 右折

やるだけ (これに結構時間かかった) (色々と考えて通せたのはうれしかった) (でも 400 点)

#include <cstdio>
#include <vector>
#include <map>
#include <utility>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

using P = std::pair<int, int>;
using ll = long long;

int N, M;
char fld[2010][2010];
int horB[2010][2010], horE[2010][2010];
bool fl[2010][2010];
ll ans;

bool safe( int i, int j )
{ return i >= 0 && i < N && j >= 0 && j < M; }

int main()
{
  scanf( "%d%d", &N, &M );
  rep( i, N )
    scanf( "%s", fld[i] );

  rep( i, N )
  {
    int b = -1;

    rep( j, M )
    {
      if( fld[i][j] == '#' )
        b = -1;
      if( b == -1 && fld[i][j] != '#' )
        b = j;

      horB[i][j] = b;
    }

    b = -1;

    for( int j = M-1; j >= 0; --j )
    {
      if( fld[i][j] == '#' )
        b = -1;

      if( b == -1 && fld[i][j] != '#' )
        b = j;

      horE[i][j] = b;
    }
  }  

  /*rep( i, M ) rep( j, N )
    printf( "%d%c", horB[i][j], j==N-1?'\n':' ' );
  rep( i, M ) rep( j, N )
    printf( "%d%c", horE[i][j], j==N-1?'\n':' ' );
*/
  rep( j, M )
  {
    int b = -1;
    ll sum = 0;

    rep( i, N ) 
    {
      //printf( "(%d,%d)\n", i, j );

      if( b == -1 && fld[i][j] != '#' )
        b = i;

      if( fld[i][j] != '#' && safe( i, j+1 ) && fld[i][j+1] != '#' )
        sum += horE[i][j]-j;//, printf( "+sum: %d\n", horE[i][j]-j );

      if( fld[i][j] != '#' && safe( i, j-1 ) && fld[i][j-1] != '#' )
        sum += j-horB[i][j];//, printf( "+sum: %d\n", j-horB[i][j] );

      if( b != -1 && ( i == N-1 || fld[i+1][j] == '#' ) )
      {
        //printf( "i-b: %d, i: %d, b: %d\n", i-b, i, b );
        //printf( "sum: %lld\n", sum );
        ans += (i-b)*sum;
        sum = 0;
        b = -1;
      }
    }
  }

  printf( "%lld\n", ans );
  
  return 0;
}
D うほょじご

互除法の亜型。規則はない。制約が 1 <= M,N < 1000 なので全て試せる。手順を逆順に探索。(4 つの候補が出る) 手順の途中でも < 1000 は保たれることに注意。rev の逆写像に注意。(‘0’ を前から足せるだけ複数解のある多価関数。逆に rev^{-1}(90) など後ろに 0 があると解はない。(後ろに 0 があるということは rev する前に前に 0 があるということだが、10 進法ではそのようには見ない。)

#include <algorithm>
#include <iostream>
#include <string>
#include <sstream>
#include <queue>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

using P = std::pair<int, int>;

std::vector<int> rev( int x )
{
  std::vector<int> ret;
  std::ostringstream oss;
  oss << x;
  std::string t = oss.str();
    
  //printf( "rev(%d)\n", x );

  if( x % 10 == 0 )
    return ret;

  do
  {
    int y = 0;

    for( int i = t.size()-1; i >= 0; --i )
      y *= 10, y += t[i]-'0';

    ret.push_back( y );

    //printf( "y: %d\n", y );

    if( !x )
      break;

    t = "0"+t;
  } while( t.size() <= 3 );

  return ret;
}

int N, M;
std::queue<P> que;
bool used[1010][1010];
int ans;

int main()
{
  scanf( "%d%d", &N, &M );

  rep( i, 999 )
    que.push( P( i+1, 0 ) );
  rep( j, 999 )
    que.push( P( 0, j+1 ) );

  while( !que.empty() )
  {
    P p = que.front(); que.pop();

    int x = p.first, y = p.second;
    std::vector<int> rys = rev(y), rxs = rev(x);

    if( used[x][y] ) 
      continue;

    //printf( "(%d,%d) %d\n", x, y, que.size() );

    used[x][y] = true;

    if( x+y < 1000 )
    {
      auto rxys = rev(x+y);
      
      for( auto rxy : rxys ) if( rxy < 1000 && y < 1000 && rxy <= y && !used[rxy][y] )
        que.push( P( rxy, y ) );//, printf( "push (%d,%d)\n", rxy, y );
      
      for( auto rxy : rxys ) if( x < 1000 && rxy < 1000 && rxy <= x && !used[x][rxy] )
        que.push( P( x, rxy ) );//, printf( "push (%d,%d)\n", x, rxy );
    }

    for( auto ry : rys ) if( ry < 1000 && x+y < 1000 && ry <= x+y && !used[x+y][ry] )
      que.push( P( x+y, ry ) );//, printf( "push (%d,%d)\n", x+y, ry );
    
    for( auto rx : rxs ) if( rx < 1000 && x+y < 1000 && rx <= x+y && !used[rx][x+y] )
      que.push( P( rx, x+y ) );//, printf( "push (%d,%d)\n", rx, x+y );
  }

  repi( i, 1, N+1 ) repi( j, 1, M+1 ) if( !used[i][j] )
    ++ans;//, printf( "(%d,%d)\n", i, j );

  printf( "%d\n", ans );

  return 0;
}
E 迷路

迷路の問題。bfs でなく dijkstra でやる方が高速。”待ち” の選択がポイント。そこまでの最短路を dist として、dist%K から方向 d (=‘U’ とか) で移動したいとして、どれだけ待てば良いかを nxt[d][dist%K] でもっておけば良い。(nxt[s[i]][i]=0 として、遡って 1 を足す。末尾から最初の 0 にかけては最初まで遡った後 cyclic に通る) INF は 1<<30 だと 1 つテストケース通らない。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <string>
#include <queue>
#include <limits>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()
#define clr(a,v) memset((a),(v),sizeof(a))

using ll = long long;
using P = std::pair<int, int>;
using State = std::pair<ll, P>;

const std::string URDL = "URDL";
constexpr ll dx[4] = { 0, 1, 0, -1 }, dy[4] = { -1, 0, 1, 0 };
constexpr ll INF = std::numeric_limits<ll>::max()>>2;

ll N, M, K;
char d[100010], s[1010][1010];
ll si, sj, gi, gj;
bool used[1010][1010];
std::priority_queue<State, std::vector<State>, std::greater<State> > pque;
ll nxt[4][100010];
ll di[1010][1010];

bool safe( int i, int j )
{ return i >= 0 && i < N && j >= 0 && j < M; }

int main()
{
  scanf( "%lld%lld%lld", &N, &M, &K );
  scanf( "%s", d );
  rep( i, N )
  {
    scanf( "%s", s[i] );

    rep( j, M )
    {
      if( s[i][j] == 'S' )
        si = i, sj = j;
      if( s[i][j] == 'G' )
        gi = i, gj = j;
    }
  }

  std::fill( (ll*)di, (ll*)(di+1010), INF );
  
  clr( nxt, -1 );

  rep( i, K )
    nxt[URDL.find(d[i])][i] = 0;
  
  rep( i, 4 ) 
  {
    bool fl = false;
    ll p = 1;

    for( int j = K-1; j >= 0; --j )
    {
      if( !nxt[i][j] )
        fl = true, p = 0;

      if( fl )
        nxt[i][j] = p++;
    }

    for( int j = K-1; j >= 0; --j )
    {
      if( !nxt[i][j] )
        break;

      if( fl )
        nxt[i][j] = p++;
    }
  }

  pque.push( State( 0, P( si, sj ) ) );
  di[si][sj] = 0;

  while( !pque.empty() )
  {
    State st = pque.top(); pque.pop();

    ll dist = st.first, ci = st.second.first, cj = st.second.second;

    if( di[ci][cj] < dist )
      continue;

    if( ci == gi && cj == gj )
    {
      printf( "%lld\n", dist );

      return 0;
    }

    rep( d, 4 ) if( nxt[d][dist%K] != -1 )
    {
      ll cost = 1+nxt[d][dist%K];
      ll ni = ci+dy[d], nj = cj+dx[d];

      if( safe( ni, nj ) && s[ni][nj] != '#' && di[ni][nj] > dist+cost )
      {
        di[ni][nj] = dist+cost;

        pque.push( State( di[ni][nj], P( ni, nj ) ) );
      }
    }
  }

  puts("-1");

  return 0;
}
F チーム分け

こういう dp しゅき (初見で解けると言ってない)
dp[i][j] : i 人以上のグループを作り、j 人余っている。a[i]==k なる個数 c[k] を求めておくと、dp[i][c[i]+j-k*i] += comb(c[i]+j,k*i)*(k*i)!/(i!)^k/k!*dp[i+1][j] (c[i]+j から k*i を選び、k*i を i グループ k 個に分ける。k*i を並び替え、前から i 個ずつグループを作り、グループ内外の並び替えを考える)

#include <cstdio>
#define repi(i,a,b) for(ll i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

using ll = long long;

constexpr int MAX_N = 1010;
constexpr ll mod = 998244353;

ll N;
ll a[MAX_N], c[MAX_N+1];
ll dp[MAX_N+1][MAX_N+1];
ll fact[MAX_N+1], inv[MAX_N+1];

ll powMod( ll x, ll y )
{
  ll ret = 1;

  while( y > 0 )
  {
    if( y & 1 )
      ret = (ret*x) % mod;

    y >>= 1;
    x = (x*x) % mod;
  }

  return ret;
}

int main()
{
  scanf( "%lld", &N );
  rep( i, N )
    scanf( "%lld", a+i ), ++c[a[i]];
  
  fact[0] = 1;
  rep( i, N )
    fact[i+1] = ((i+1)*fact[i]) % mod;

  rep( i, N+1 )
    inv[i] = powMod( fact[i], mod-2 );

  dp[N+1][0] = 1;

  for( ll i = N; i >= 1; --i ) rep( j, N+1 )
  {
    ll f = 1;

    // dp[1][0] = k=3i=1 dp[2][3]
    for( ll k = 0; k*i <= j+c[i]; ++k )
    {
      //printf( "dp[%d][%d](%lld) += %lld*dp[%d][%d](%lld)\n", i, j+c[i]-k*i, dp[i][j+c[i]-k*i], fact[j+c[i]]*f%mod*inv[j+c[i]-k*i]%mod*inv[k]%mod, i+1, j, dp[i+1][j] );
      dp[i][j+c[i]-k*i] = (dp[i][j+c[i]-k*i]+fact[j+c[i]]*f%mod*inv[j+c[i]-k*i]%mod*inv[k]%mod*dp[i+1][j]%mod) % mod;
      f = (f*inv[i])%mod;
    }
  }

  printf( "%lld\n", dp[1][0] );

  return 0;
}
G 移動

荷重重心。av_1+bv_2+cv_3 (0 <= a,b,c, a+b+c <= K) の作る点は Σ_{k=0}^K 3_H_k = (K+1)(K+2)(K+3)/6 である。それから重複を引く。重複は、xv_1+yv_2+zv_3 = 0 なる (x,y,z) に対して、0 <= a+x, b+y, c+z, (a+x)+(b+y)+(c+z) <= K であるような (a,b,c) の個数。max(0,-x) <= a, max(0,-y) <= b, max(0,-z) <= c であり、max(0,-x)+max(0,-y)+max(0,-z) <= a+b+c <= min(k,k-x-y-z) なので、K = min(k,k-x-y-z)-max(0,-x)-max(0,-y)-max(0,-z) として計算すれば良い。K < 0 の時は、当然和は 0。gcd 取る時は abs 取る。

#include <cstdio>
#include <cstdlib>
#include <algorithm>
#define repi(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,a) repi(i,0,a)
#define all(a) (a).begin(), (a).end()

using ll = long long;

constexpr ll mod = 998244353;

ll powMod( ll x, ll y )
{
  ll ret = 1;

  while( y > 0 )
  {
    if( y & 1 )
      ret = (ret*x) % mod;

    y >>= 1;
    x = (x*x)%mod;
  }

  return ret;
}

ll func( ll K )
{ return K < 0 ? 0 : (K+1)*(K+2)%mod*(K+3)%mod*powMod(6,mod-2)%mod; }

ll gcd( ll a, ll b )
{ return b ? gcd(b,a%b) : a; }

int Q;

int main()
{
  scanf( "%d", &Q );
  rep( i, Q )
  {
    ll a, b, c, d, e, f, k;

    scanf( "%lld%lld%lld%lld%lld%lld%lld", &a, &b, &c, &d, &e, &f, &k );

    ll x = d*e-c*f, y = a*f-b*e, z = -a*d+b*c;
    ll g = gcd(llabs(x),gcd(llabs(y),llabs(z)));

    x /= g, y /= g, z /= g;

    printf( "%lld\n", (func(k)-func(std::min(k,k-x-y-z)-std::max(0ll,-x)-std::max(0ll,-y)-std::max(0ll,-z))+mod)%mod );
  }
  
  return 0;
}
H タイル張り

わからん。とりあえずわかりやすそうなこの辺りのコードを眺めていたがよくわからん (情報量が不足)

mujin-pc-2018.contest.atcoder.jp

mujin-pc-2018.contest.atcoder.jp

最初のやつをメインに見てたけど、集合の冪集合が最大 91 通りしかないらしく、それの隣接行列 (?) (列を 1 列経る毎の状態遷移) を W 乗している??? (繰り返し 2 乗法で高速化)
なんで集合の集合を考えなきゃいけないのかしっくりこない。dfs では そのまま or 2 個連続でフラグが立っているところを潰している っぽいが、そもそもこのフラグが editorial に沿っているのかもわからない (つじつまが合わない?)

教えて強い人()

pV 図上の吸発熱の判定

こんにちは,Ysmr-Ry です.
今回は僕が高校のとき曖昧だったことをまとめてみたいと思います.
意外とコレがなかなか教えてくれなかったりする.

話題は高校の熱力学です.(高校の熱力学の基礎的な知識を仮定します)
pV 図というのを書くじゃないですか.囲む面積が仕事になるので…
大学に行くとエントロピーという状態量を習って,TS 図とかいうのもかけたりするんですね.
それは囲む面積が熱を表すんです.pV と TS どちらも次元がエネルギーになってるからなんですけど.

つまり高校ではエネルギーを供給する熱と仕事のうち仕事にフォーカスを当てて,そちらをわかりやすいようにしてるのです.物事というのは往々にしてトレードオフであり,そうすると熱がわかりにくくなってしまうんですね.

今回はそのわかりにくい熱のお話です.

熱力学とは

熱力学というのは系の熱平衡状態を仮定して,状態量というマクロな少数パラメータで熱力学的状態を記述し,その上での法則を考えようというものです.モル数を固定すると,状態量としては,圧力  p,体積  V絶対温度  T,内部エネルギー  U などがあります.(高校で出てくるのはたぶんこれだけ)
大学に行くともっとたくさん出てくるのですが… そのうちで最も重要なのがエントロピー  S です.まぁざっくり言うと熱を司り,エントロピーが変化しないのなら熱がなく断熱過程です.(ここではその程度の認識で十分)
で,状態量の自由度は 2 なのです.状態量としてはいっぱいあるのですが,そのうちの 2 つの量がわかれば,他の状態量は既知の状態量との関係式をもとに原理的に決定できるのです.(僕はこういう理解なのですが間違ってたら教えてください)
冒頭で pV 図や TS 図が出てきましたが,これらは 2 次元平面です.平面上の 1 点を与えれば,実数 2 つ分の情報量を与えることができ,それが熱力学的な状態と 1 対 1 に対応しているというわけです.
で,高校での主役は pV 図です.その上での熱の話をしたいのです.熱を司るのはエントロピーなのです.
あらわに目にする状態量はこの場合  p, V なわけですが,その裏で  S が働いているという構図です.
 p, V S の関係を熱力学的な法則から考えたいのです.
まずは,"あらわでない状態量が裏で働いている" 様を体感してもらいましょう.

等温過程

あらわでない状態量として,絶対温度  T を考えます.
例えば, T=298K の状態を集めてくると pV 図上ではどうなるのでしょうか.熱力学的状態は pV 図上の 1 点に対応しているので,pV 図上の領域になることは容易にわかります.
それに加えて,自由度が 2 であることも考えると,そのうちの 1 つである  T を固定しているため,実質的な自由度は 1 になります.この場合自由度を次元と置き換えることができ,1 次元,すなわち曲線になることが予想されます*1
高校物理の熱力学で考える作業物質は理想気体がほとんどですから,具体的な曲線は Boyle の法則が与えます.
Boyle の法則により,等温であるとき, pV が一定となります.
 xy 座標平面で言えば  xy=const. であり,これは反比例のグラフです.(直角双曲線) (曲線にちゃんとなった)
また,理想気体の状態方程式から, pV の値は  T に比例します.
よって, T をとびとびに増やしてやると,このように平面を覆っているのがわかります.(一つ一つが  T が一定な状態を集めた双曲線)

右上に行くほど一定である  T の値が大きくなります.
これが "あらわでない状態量が裏で働いている" 例です.

断熱過程

では,裏で働いている状態量がエントロピー  S である場合を考えましょう.
熱の議論をするためには基本的に熱力学第一法則から出発しなければなりません.
熱力学第一法則は内部エネルギーの変化を  dU,気体がする仕事を  d'W,熱を  d'Q として, d'Q=dU+d'W となります.高校では  \Delta を使って書いていましたが,厳密には微小量です*2.ここで  d'W = pdV です.作業物質は理想気体なので, dU = nc_VdT です.( c_V は定積モル比熱)
すると,
\begin{align*}
\frac{d'Q}{T} &= nc_V\frac{dT}{T}+\frac{p}{T}dV \\
&= nc_V\frac{dT}{T}+nR\frac{dV}{V}
\end{align*}
1 行目から 2 行目へは状態方程式  pV=nRT を使いました.
両辺を積分すると,
 \begin{align*}
\int \frac{d'Q}{T} &= nc_V\log T+nR\log V \\
&= nc_V\log TV^\frac{R}{c_V} \\
&= nc_V\log \frac{pV^\frac{c_V+R}{c_V}}{nR}
\end{align*}
これが断熱で 0 ならば,Mayer の関係  c_P = c_V+R と比熱比  \gamma=\frac{c_P}{c_V} を使って, pV^\gamma = const. です.
Poisson の式として教科書にぽつりと載っていたりします.この式の表す pV 図上の曲線を断熱曲線といいますが,そのある接点での接線の傾きは, \frac{dp}{dV} = -\gamma\frac{p}{V} になります.接点の表す状態から接線に沿って進むことは断熱過程を表します.その点からの過程においてこのラインが吸発熱の境界線になります.
ちなみにエントロピーの定義が  \int \frac{d'Q}{T} だったりするので,いわば "裏でエントロピーが働いている" 曲線なわけです.

変形の仕方を少し変えると,
\begin{align*}
d'Q &= nc_VdT+pdV \\
&= nc_V\frac{Vdp+pdV}{nR}+pdV \\
&= \frac{c_V}{R}Vdp+\left(\frac{c_V}{R}+1\right)pdV
\end{align*}
となります.1 行目から 2 行目へは,状態方程式から, dT = \frac{d(pV)}{nR} = \frac{Vdp+pdV}{nR} を使いました.
 d'Q > 0 で吸熱のときは,
 \begin{align*}
\frac{c_V}{R}Vdp+\left(\frac{c_V}{R}+1\right)pdV &> 0 \\
\frac{dp}{dV} &> -\gamma\frac{p}{V}
\end{align*}
となります.
 d'Q < 0 で発熱のときは  \frac{dp}{dV} < -\gamma\frac{p}{V} です.
つまり,ある 1 点の状態からの過程は,そこを接点とする断熱曲線の接線よりも右上ならば吸熱,左下ならば発熱なのです.

図の紫の曲線が断熱曲線です.オレンジの点の状態から,紫の接線に沿って進めば断熱過程,青い線に沿って進めば吸熱,赤い線に沿って進めば発熱です.

つまり,

このように,青い吸熱の領域と赤い発熱の領域に分けられているイメージです.(この場合はオレンジの 1 点に関しての話です.)

おまけ

以下のような熱サイクルを考えてみましょう.

f:id:Ysmr_Ry:20171231030407p:plain

先ほどの話から,左,上に進む過程はそれぞれ発熱,吸熱と決まりますが,斜めの過程はそうとは限りません.
特にこの場合,途中で切り替わるようにできています*3.切り替わるのはこの斜めの線と断熱曲線が接している点です.

右下の方のオレンジの点が切り替わる点です.試しにそれよりも左上のもう 1 つのオレンジの点を見てみましょう.
この点を通る断熱曲線の接線がオレンジの線です.この線よりも熱サイクルの青の斜めの線の方が右上にあるので,この点では吸熱です.接点をすぎた先の点では発熱に切り替わります.
これは, \frac{d^2p}{dV^2} = \gamma\frac{p}{V^2} > 0 より断熱曲線が下に凸なことによります.

理想気体の場合  \gamma=\frac{5}{3} なので*4,この斜め線の方程式と Poisson の式を連立すれば切り替わる接点の座標がわかります.計算すると  (\frac{15}{8}V_0,\frac{9}{8}p_0) であり,この熱サイクルの熱効率は  \frac{48}{97} となります.(間違ってたら教えてください)

*1:ここで直線といいたくなるかもしれませんが,途中折れ曲がっていてもよく,折れ線を限りなく細かくしていくと曲線になるので一般には曲線です.

*2:僕の先生は熱力学第一法則を  \Delta で書くやつを信用するなとまで言っていました.

*3:過去に東工大がこのサイクルを出題したそうです.途中で切り替わるのは出題ミスではと僕の先生が言ってました.

*4:ここらへんは統計力学の結論です.