2016-08-16

Convex-Hull Trick

概要

直線集合 $L=\{f_i(x)=a_ix+b_i\}$ に対して、「 $L$ への直線 $f_k(x)=a_kx+b_k$ の追加」(以下、追加クエリ)と、「ある $x$ に対し、 $L$ の中で最小値(あるいは最大値)を取るような直線の値を求める」(以下、最小値クエリ)という2つのクエリを効率的に行うことが出来ます。

使える状況

DPの漸化式を整理したときなどにおいて、 $\displaystyle min_{1\leq i\leq n}(a_ix+b_i)$ といった式が出てきたときに、Convex-Hull Trickを用いることで効率的に値を求めることが出来ます。

説明

ここでは最小値を求めるときのみを説明します(最大値を求めるときは上⇔下、増加⇔減少など、文章を補って読んでください)。

クエリ

まず、それぞれのクエリを見ていきます。

直線集合 $L$ に直線 $f_k(x)=a_kx+b_k$ を追加する

f:id:satanic0258:20160815045239p:plain
上図は、始め直線集合 $L$ が $L=\{y=2x+2, y=x+3, y=-x+10\}$ であった状態(左)から、そこに新たな直線 $y=-2x+15$ を追加した状態(右)へと遷移している様子を表しています。
(なお、グラフの描画にはGrapesというソフトを用いています。*1 )

ある $x$ に対し、 $L$ の中で最小値(あるいは最大値)を取るような直線の値を求める

上図左において、 $x=0$ のとき、 $3$ 本の直線の中で最も小さい値を取るものは $y=2x+2$ であり、その時の値は $2$ となります。
同様に $x=2$ では $5$ 、 $x=6$ では $4$ ということも分かります。

その後、赤色の直線を追加すると、 $x=0$ や $x=2$ のところでは変わらず $2, 5$ となっていますが、 $x=6$ のところでは $3$ と、直線を追加する前とは値が異なっていることが分かります。

このクエリを愚直に実装しようとすると、直線の本数を $N$ 本、最小値クエリの回数を $Q$ 回としたとき、最小値クエリを求める度に各直線における値を調べる必要があるので、計算量は $O(NQ)$ となります。

しかし、より速くクエリをこなす方法があるため、それについて考えていきます。

性質

もう一度先ほどの図をよく見てみると、2つの性質が隠れていることが分かります。

1.最小値を取る直線は左から順に傾きが減少している

f:id:satanic0258:20160815052207p:plain

ちなみにこのことを簡単に証明すると、上図の状態で傾き $2$ の直線の区間と傾き $1$ の直線の区間の間に傾き $0$ の直線の区間があるようなときに反例となりますが、

直線 $y=3$ を追加したとき、傾きは $2\to 0\to -2$ と遷移する
直線 $y=5$ を追加したとき、傾きは $2\to 1\to 0\to -2$ と遷移する

となります。よって、最小値を取る直線は左から順に傾きが減少している事がわかりました。

2.不必要な直線が無いとき、ある $x$ での直線の取る値は、その値を直線の傾き順に並べると下に凸となる

文章にすると少しわかりにくいですが、図にすると次のような感じです。
f:id:satanic0258:20160815163318p:plain
今、 $x=3$ のときについて考えています。このとき、直線の傾きを横軸、 $x=3$ のときの直線の取る値を縦軸とすると、次のようになります。
f:id:satanic0258:20160815163328p:plain

この性質についても簡単に証明すると、次の図左のような、直線集合 $\{y=-2x+3, y=-x+1, y=x+2, y=2x\}$ を考えたとき、 $x=0$ において直線の取る値をそのまま図にすると、図右のようになります。
f:id:satanic0258:20160815171620p:plain
しかし、図左に緑色で示した直線 $y=x+2$ は、 $x$ をどのようにとっても直線集合の中で最小値を取ることがありません。したがって、この緑色の直線を直線集合から取り除くことによって、図右の青色の折れ線のように下に凸となります。故に、不必要な直線を取り除いてしまえば、直線の取る値は傾き順に並べると下に凸となります。

これらの性質から、直線集合を、傾きをキーとしてソートしておくことにより、最小値クエリに二分探索を用いることが出来るため、 $N$ 本の直線に対する $Q$ 回の最小値クエリの計算量は $O(Q\log N)$ となります。

また、直線の追加クエリは、ソート済みの直線集合に追加する、ということから、先ほどと同様に直線の挿入位置を調べるのに二分探索を用いることが出来るため、全体の直線の数を $N$ 本とすると、全ての追加クエリの計算量は $O(N\log N)$ となります。

不必要な直線を取り除く

他に、先ほどの性質を満たすためには「不必要な直線を取り除く」という動作も適宜行う必要があります。

それでは、どのようなとき、どうなっていれば「不必要」と判断することが出来るのか、それについて考えていきます。

まず、前提として直線集合がすでに不必要な直線を全て取り除いてあると仮定します。
このように仮定してしまえば、次に直線が不必要になるときは、直線を追加する際であると分かります。
よって、「直線の追加クエリ」を行うときに、「追加する事によって不必要になる直線を調べる」ようにすれば良いことが分かりました。

次に、どのような直線が不必要かについてですが、ここで次の図を考えます。
f:id:satanic0258:20160815190352p:plain
これは不必要な直線が無い状況です。なお、直線 $l_i:y = a_ix+b_i$ とすると傾き $a_i$ は $a_1\geq a_2\geq a_3$ となっています。

このとき、直線 $l_1$ と直線 $l_2$ の交点を赤の点、直線 $l_2$ と直線 $l_3$ の交点を青の点、直線 $l_3$ と直線 $l_1$ の交点を緑の点とすると、これらの交点の $x$ 座標は、
赤の点 $\leq$ 緑の点 $\leq$ 青の点
となっていることが分かります。

一方、次の図を考えます。
f:id:satanic0258:20160815191214p:plain
これは不必要な直線がある状況です。直線の設定は先ほどと変わりません。

このとき、交点の $x$ 座標は、
青の点 $\leq$ 緑の点 $\leq$ 赤の点
となっていることが分かります。

このように、不必要な直線があるときとないときで、相異なる $2$ 直線の交点の位置関係が変わっています。

よって、 $3$ 直線のうち相異なる $2$ 点を選んで、それらの $x$ 座標を比較することにより不必要かどうかを判断できることが分かりました。

ここでは、赤の点と青の点を選んで考えます(どの $2$ 点を選んでも判断はできますが、蟻本*2と同じ選び方としておきます)。

まず、赤の点の $x$ 座標を求めます。
赤の点は直線 $y=a_1x+b_1$ と直線 $y=a_2x+b_2$ の交点なので、次の連立方程式を解けば求めることが出来ます。

$\displaystyle \begin{cases}y=a_1x+b_1\\ y=a_2x+b_2 \end{cases}$

これを解くと(簡単のために $x$ についてのみ)、

$\displaystyle \begin{eqnarray} a_1x+b_1 &=& a_2x+b_2 \\ (a_1-a_2)x &=& b_2-b_1 \\ x &=& \cfrac{b_2-b_1}{a_1-a_2} \end{eqnarray}$

よって、赤の点の $x$ 座標は $x = \cfrac{b_2-b_1}{a_1-a_2}$ であることが分かりました。

同様に、青の点の $x$ 座標を求めます。
青の点は直線 $y=a_2x+b_2$ と直線 $y=a_3x+b_3$ の交点なので、次の連立方程式を解けば求めることが出来ます。

$\displaystyle \begin{cases}y=a_2x+b_2\\ y=a_3x+b_3 \end{cases}$

これを解くと(簡単のために $x$ についてのみ)、

$\displaystyle \begin{eqnarray} a_2x+b_2 &=& a_3x+b_3 \\ (a_2-a_3)x &=& b_3-b_2 \\ x &=& \cfrac{b_3-b_2}{a_2-a_3} \end{eqnarray}$

よって、青の点の $x$ 座標は $x = \cfrac{b_3-b_2}{a_2-a_3}$ であることが分かりました。

あとはこれらの大小関係について調べればよいことになります。ここでは、「不必要であるかどうか」という不等式を作ります。
そのためには、先ほどの議論により、
青の点の $x$ 座標 $\leq$ 赤の点の $x$ 座標
とすれば良いです。

今、 $a_1\geq a_2\geq a_3$ であることから、 $a_1-a_2\geq 0$ 及び $a_2-a_3\geq 0$ であるので、

$\displaystyle \begin{eqnarray} \cfrac{b_3-b_2}{a_2-a_3} &\leq & \cfrac{b_2-b_1}{a_1-a_2} \\ (b_3-b_2)(a_1-a_2) &\leq & (b_2-b_1)(a_2-a_3) \\ (b_3-b_2)(a_2-a_1) &\geq & (b_2-b_1)(a_3-a_2) \end{eqnarray}$

よって、この式を満たしているとき、直線 $l_2:y=a_2x+b_2$ が不必要であることが分かりました。

以上をまとめると $2$ つのクエリは、

直線の追加するクエリ $\cdots$ 直線を追加する位置を二分探索で求める。その際、追加によって不必要となる直線は取り除く。
最小値を求めるクエリ $\cdots$ 最小値を取る直線を二分探索で求める。

とすることで、直線 $N$ 本、最小値クエリ $Q$ 回とすると、全体の計算量は $O((N+Q)\log N)$ となることが分かりました。

特殊な状況における更なる高速化

上記では、どのようなクエリ群であっても最小値を求められるものを説明しました。

ここでは、そのクエリ群がある条件を満たしているようなときに、より高速化が可能であるため、それについて解説していきます。

追加クエリにおける直線の傾きが単調であるとき

「追加クエリにおける直線の傾きが単調」とは、始め直線集合が空である状態から、直線 $l_i:y=a_ix+b_i$ を、 $i=1,2,3,\cdots$ と追加していくときに、その傾き $a_i$ が $a_1\geq a_2\geq a_3\geq \cdots$ (あるいは $a_1\leq a_2\leq a_3\leq \cdots$ )であることを指します。

このとき、先ほど述べたように直線の挿入位置を二分探索で調べる必要はなく、常に先頭、あるいは末尾に追加していくことになるため、直線の本数 $N$ 本とすると、追加クエリの計算量は $O(N)$ となります。

最小値クエリにおける $x$ が単調であるとき

「最小値クエリにおける $x$ が単調」とは、最小値クエリにおいて与えられる $x_i$ が $x_1\geq x_2\geq x_3\geq \cdots$ (あるいは $x_1\leq x_2\leq x_3\leq \cdots$ )であることを指します。

このとき、最小値クエリにおける $x$ を進めると、直線集合の先頭の直線が最小値を取らなくなるときがあります。すると、 $x$ は後戻りしないため、もうこの直線は使わないことが分かります。よって、直線集合の先頭からも直線を取り除くようにすることで、最小値クエリを行うときは常に直線集合の先頭の直線が最小値を取ることになるため、最小値クエリ $1$ 回にかかる計算量は $O(1)$ となります。
f:id:satanic0258:20160815211119p:plain
(図における色付きの直線は最小値を取らなくなったもの)

上記 $2$ パターンを両方満たしているとき

つまり、「追加クエリにおける直線の傾きも単調」であり「最小値クエリにおける $x$ も単調」であるようなときを指します。

このときは、直線の本数 $N$ 本、最小値クエリの回数を $Q$ 回とすると、その計算量は $O(N+Q)$ となるため、条件が無い時に比べて $\log N$ 分速くクエリをこなすことが出来ます。

コード(C++)

※注：このコードは追加クエリが単調であると仮定したものです。
~~(力不足で一般的な状況のものが出来ていません…。書きかけのものを一応ここに残しておきます)~~
[2018-02-27 追記] 一般的なクエリに対応したものをココに書きましたが，若干遅そうなのでまだα版です

#include <algorithm>
#include <functional>
#include <utility>
#include <vector>

template<typename T>
class ConvecHullTrick {
private:
	// 直線群(配列)
	std::vector<std::pair<T, T>> lines;
	// 最小値(最大値)を求めるxが単調であるか
	bool isMonotonicX;
	// 最小/最大を判断する関数
	std::function<bool(T l, T r)> comp;

public:
	// コンストラクタ ( クエリが単調であった場合はflag = trueとする )
	ConvecHullTrick(bool flagX = false, std::function<bool(T l, T r)> compFunc = [](T l, T r) {return l >= r; })
		:isMonotonicX(flagX), comp(compFunc)  {
		lines.emplace_back(0, 0);
	};

	// 直線l1, l2, l3のうちl2が不必要であるかどうか
	bool check(std::pair<T, T> l1, std::pair<T, T> l2, std::pair<T, T> l3) {
		if (l1 < l3) std::swap(l1, l3);
		return (l3.second - l2.second) * (l2.first - l1.first) >= (l2.second - l1.second) * (l3.first - l2.first);
	}

	// 直線y=ax+bを追加する
	void add(T a, T b) {
		std::pair<T, T> line(a, b);
		while (lines.size() >= 2 && check(*(lines.end() - 2), lines.back(), line))
			lines.pop_back();
		lines.emplace_back(line);
	}

	// i番目の直線f_i(x)に対するxの時の値を返す
	T f(int i, T x) {
		return lines[i].first * x + lines[i].second;
	}

	// i番目の直線f_i(x)に対するxの時の値を返す
	T f(std::pair<T, T> line, T x) {
		return line.first * x + line.second;
	}

	// 直線群の中でxの時に最小(最大)となる値を返す
	T get(T x) {
		// 最小値(最大値)クエリにおけるxが単調
		if (isMonotonicX) {
			static int head = 0;
			while (lines.size() - head >= 2 && comp(f(head, x), f(head + 1, x)))
				++head;
			return f(head, x);
		}
		else {
			int low = -1, high = lines.size() - 1;
			while (high - low > 1) {
				int mid = (high + low) / 2;
				(comp(f(mid, x), f(mid + 1, x)) ? low : high) = mid;
			}
			return f(high, x);
		}
	}
};

コードの説明

ほぼ先ほどの議論をそのままコードに書いています。
なお、関数はそれぞれ、

${\tt check}$ 関数 $\cdots$ 直線の不必要性を判断する
${\tt add}$ 関数 $\cdots$ 直線集合に直線を追加する
${\tt f}$ 関数 $\cdots$ 直線 $l_i$ に対する $f_i(x)$ を返す
${\tt get}$ 関数 $\cdots$ 直線集合内で $x$ のときの最小値を返す

となっています。

また、使うときには

ConvexHullTrick<int> cht(true);

といった感じで宣言します。
このとき、第一引数を ${\tt true}$ か ${\tt false}$ にすることでそれぞれ最小値(最大値)クエリが単調であるか否かを指定することが出来ます(省略可)。
また、第二引数に、

ConvexHullTrick<int> cht(true, std::less<int>());

としたり、

ConvexHullTrick<int> cht(true, [](T l, T r) {return l <= r; });

とすることで最大値クエリをこなすことが出来ます。

具体例

実際に直線を追加していって、その後の最小値クエリの違いを見ていきます。

ソースコード

※ヘッダ、 ${\tt ConvexHullTrick}$ クラスは省略

int main() {
	ConvecHullTrick<int> cht(false);
	cht.add(2, 0);
	std::cout << "Add y=2x" << std::endl;
	std::cout << "min(f(-2)) = " << cht.get(-2) << std::endl;
	std::cout << "min(f(2)) = " << cht.get(2) << std::endl << std::endl;

	cht.add(0, -1);
	std::cout << "Add y=-1" << std::endl;
	std::cout << "min(f(-2)) = " << cht.get(-2) << std::endl;
	std::cout << "min(f(2)) = " << cht.get(2) << std::endl << std::endl;

	cht.add(1, 1);
	std::cout << "Add y=x+1" << std::endl;
	std::cout << "min(f(-2)) = " << cht.get(-2) << std::endl;
	std::cout << "min(f(2)) = " << cht.get(2) << std::endl << std::endl;

	cht.add(-1, 0);
	std::cout << "Add y=-x" << std::endl;
	std::cout << "min(f(-2)) = " << cht.get(-2) << std::endl;
	std::cout << "min(f(2)) = " << cht.get(2) << std::endl << std::endl;
	return 0;
}

出力

Add y=2x
min(f(-2)) = -4
min(f(2)) = 4

Add y=-1
min(f(-2)) = -4
min(f(2)) = -1

Add y=x+1
min(f(-2)) = -4
min(f(2)) = -1

Add y=-x
min(f(-2)) = -4
min(f(2)) = -2

上記例で行ったクエリを図にしてみると、次のようになっています。
f:id:satanic0258:20160816170731p:plain
これで、確かに正しく出力出来ていることが分かりました。

備考

前述のように、上記コードは全ての状況に対応出来るものではないので注意してください。

あとがき

使える機会は式を見て割と分かりやすいですが、実装は慣れていないと難しそうです。

そのため、問題を色々調べてみました(まだ全て解いてないです)。
3709 -- K-Anonymous Sequence
No.409 ダイエット - yukicoder
I: Live Programming - Japan Alumni Group Summer Camp 2015 Day 4 | AtCoder
Contest Page | CodeChef
Contest Page | CodeChef
Problem - E - Codeforces

このConvex-Hull Trickを用いることにより、計算量を $O(NQ)$ から $O((N+Q)\log N)$ へと高速化出来るため、覚えておくとどこかで役に立つと思います。

*1:http://www.osaka-kyoiku.ac.jp/~tomodak/grapes/

*2:

プログラミングコンテストチャレンジブック [第2版]　?問題解決のアルゴリズム活用力とコーディングテクニックを鍛える?

作者: 秋葉拓哉,岩田陽一,北川宜稔
出版社/メーカー: マイナビ
発売日: 2012/01/28
メディア: 単行本（ソフトカバー）
購入: 25人クリック: 473回
この商品を含むブログ (35件) を見る

2016-05-31

SRM691(Div1)に参加しました！が…

競プロ参加競プロ参加-SRM

初めてのDiv1だ｡ﾟ(ﾟ´ω`ﾟ)ﾟ。
— satanic@C++ (@satanic0258) 2016年5月30日

今回は色々意外なことがあったので日記として書いておきます。
(この記事に解説は含まれないので予めご了承下さい)

私が前回参加したのはSRM688だった、ということで気づいたら2回分参加忘れていたようです…

そういえばSRMのシステムあんまり覚えてない(白目)
— satanic@C++ (@satanic0258) 2016年5月30日

そのSRM688では、割といい成績を残せ(てしまっ)たので、SRM2回目にしてDiv1という魔境(？)に放り込まれることに。
ヒヤヒヤしていたらすぐに時間になってしまいました。

Coding Phase

Easy(300)

問題を解読するのにやや手こずりながらも色々考察をしてみる…
f:id:satanic0258:20160531224646j:plain
(テストケースを書いてみたり、DPで出来…ないなぁと色々やっています)

…その考察の結果、

テストケース通らず悲しみ
— satanic@C++ (@satanic0258) 2016年5月30日

結局、時間内に通すことは出来ませんでした…

個人的にグラフ理論関係の問題がまだ不慣れな部分もあってか、あまりどういうことしたら良いかを見出すことが出来ませんでしたね…(今後の目標ポイント)

Medium(500)

Hard(900)

Easyに下手に時間を掛けてしまい、気づいたら時間が終わってしまったので開けてすらいません。
点数が減ることを気にする以前に、解けなかったら意味がないのである程度考えてもわからなさそうだったら、さっさと他の開けちゃった方がいいですね…(今後の目標ポイント２)

(ちなみに、問題などを初めから順に埋めていこうという癖があるようで、大学受験時のセンター試験なども解く順番を替えていたらもう少し点が取れたのかもしれなかったですね…)

Intermission

結局何もできず終わりの鈴がなってしまいました…。

room内でeasy3人しか通してないのやばすぎでしょ
— satanic@C++ (@satanic0258) 2016年5月30日

しかし、roomの状況を見ると以外と他の人もあまりいい調子ではなさそう。

これで次のChallenge成功したら大逆転のチャンスが…？

…ということで、なんとか爪痕を残そうと改めて意気込みました(半ば諦めてはいましたが…↓)。

challengeで点稼ご(無理)
— satanic@C++ (@satanic0258) 2016年5月30日

Challenge Phase

Challenge Phaseが始まって数秒、突然Hardを通していた4人のコードが次々と撃ち落されていく光景が…

「なんだここは…恐ろしいroomだ…」

と思いながら、他の人のeasyはまだ通ったままになっていたので、粗探しをしまくりました(悪い顔)。
具体的には、Challenge Phaseが始まってからすぐにその人たちのコードを急いで写経して、テストケースをいくつか投げてみたりしていました。

すると、手元で書いていたケースが通らないコードを発見！
前回の反省を活かして、ちゃんとChallengeは出来るようにしておいたので、早速Challenge！

challenge初めてにして成功した！
— satanic@C++ (@satanic0258) 2016年5月30日

見事Challengeに初成功しました！

そのおかげで順位が一気に上がり、ガッツポーズをしながらChallenge Phaseを終えました…

最終結果

スコア

問題	満点スコア	獲得スコア	かかった時間
Easy	300	0.00	--:--
Medium	500	0.00	--:--
Hard	900	0.00	--:--

Challenge数	Challenge成功数
1	1

計:50.00

レート

ここでレートなんですが…

1231(青)→1524(黄)

ﾌｧｯ!? pic.twitter.com/WoG8AwIUsA
— satanic@C++ (@satanic0258) 2016年5月30日

…もうこのとき「えっ！？嘘でしょ！？」とか素で叫んでしまいました…。
まさか無→青→黄と2回でトントン遷移してしまうとは…。

恐らく、0完で終わってそのまま0点だった人が多かったんじゃないですかね…？
そのため順位上では黄色コーダーの人を何人か抜き、こういう結果になったんじゃないかと予想しています…。

あとがき

たとえレートが上がろうが解法全然わかんないんですけど…
— satanic@C++ (@satanic0258) 2016年5月30日

しかし、レートは確かに上がりましたが、これは自分が圧倒的成長した訳でもなく、運が良かっただけなのでなんだか複雑な心境に…
解けなかったことには変わりないので、まだまだ頑張る必要があるようですね…！
数字の結果に惑わされず、これからも精進していきたいと思います👍

2016-04-29

繰り返し自乗法

アルゴリズムアルゴリズム-か

概要

$N^P \pmod{M}$ の値を効率的に解きます。

使える状況

「 $N^P$ を $M$ で割ったあまり」が欲しいときに、 $P$ があまりにも大きいとその計算だけで時間がかかってしまうことがあります。そのような時には、この繰り返し自乗法を用いることで素早く求めることが出来ます。

説明

$N^P \pmod{M}$ の値を求めたいとき、愚直な方法として次のような手法を考えることが出来ます。

using ll = long long int;
ll PowMod(ll N, ll P, ll M){
    ll ans = 1;
    for(ll i=0; i<P; ++i){
        ans *= N;
        ans %= M;
    }
    return ans;
}

しかし、上記のような方法だと $P$ 回ループしてしまうため、 $P=10^9$ などの場合では非常に時間がかかってしまいます。

そこで使うのが、繰り返し自乗法です。

繰り返し自乗法は、

「 $N$ を一回一回掛けてmodを取るのでは時間がかかるので、まとめて掛けてしまおう！」

という考えをもとにしたアルゴリズムです。

この「まとめて掛ける」というのがどういうことか、具体的に $P=20$ のときについて考えながら説明します。

以下、求める値を $X$ 、 $N^{(i)}:=N^i \pmod{M}$ とし、また、" $\times$ "の計算をした時には $M$ で割った余りを求める操作も一緒に行っているものとします。

ボトムアップ(小さい値から大きい値へと計算していく)な考え方

先ほどの愚直な方法では、次のような計算をしていることになります。
$\displaystyle \begin{eqnarray} X &=& N^{(1)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\\ &=& N^{(2)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\\ &=& N^{(3)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\\ &&\vdots\\ &=& N^{(19)}\times N^{(1)}\\ &=& N^{(20)} \end{eqnarray}$
確かに、 $P$ 回、つまり $20$ 回分の計算をしていることがわかりますね。

ここで、上で計算しているように、
$N^{(2)} = N^{(1)}\times N^{(1)}$
となっていることから、
$\displaystyle \begin{eqnarray} X &=& N^{(1)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\times N^{(1)}\\ &=& (N^{(1)}\times N^{(1)})\times (N^{(1)}\times N^{(1)})\times \cdots \times (N^{(1)}\times N^{(1)})\\ &=& N^{(2)}\times N^{(2)}\times \cdots \times N^{(2)}\\ \end{eqnarray}$

とまとめてしまい、あとは、
$\displaystyle \begin{eqnarray} X &=& N^{(2)}\times N^{(2)}\times N^{(2)}\times N^{(2)}\times \cdots \times N^{(2)}\\ &=& N^{(4)}\times N^{(2)}\times N^{(2)}\times \cdots \times N^{(2)}\\ &=& N^{(6)}\times N^{(2)}\times \cdots \times N^{(2)}\\ &&\vdots\\ &=& N^{(18)}\times N^{(2)}\\ &=& N^{(20)} \end{eqnarray}$

と出来るので、計算回数は最初の $N^{(2)}$ を計算する $1$ 回と、まとめた後の $10$ 回となっているので、見事に半分にすることが出来ていますね。

これをさらに、
$N^{(4)}=N^{(2)}\times N^{(2)}$
とすることにより $N^{(4)}$ を求め、同様に $N^{(8)}$ を求め、…と繰り返すことで、最終的には
$\displaystyle \begin{eqnarray} X &=& N^{(16)} \times N^{(4)}\\ &=& N^{(20)} \end{eqnarray}$

とまで式を短くすることが出来ます。
このとき、次の $5$ 回しか計算をしていないことが分かりますね。

$N^{(2)}$ を求める( $N^{(1)}\times N^{(1)}$ )
$N^{(4)}$ を求める( $N^{(2)}\times N^{(2)}$ )
$N^{(8)}$ を求める( $N^{(4)}\times N^{(4)}$ )
$N^{(16)}$ を求める( $N^{(8)}\times N^{(8)}$ )
$N^{(20)}$ を求める( $N^{(16)}\times N^{(4)}$ )

このように、 $N$ の自乗を繰り返し計算した結果を用いて計算を行っているので、「繰り返し自乗法」と言うんですね。

トップダウン(大きい値から小さい値へと計算していく)な考え方

先ほどの方法とほとんど変わりませんが、実装が(少しだけ)簡単な方法もあります。
おそらく、一般的に「繰り返し自乗法」と言ったときにはこちらの方を指すのではないでしょうか。

この方法は、次のように再帰的な式で書き表すことが出来ます。
$\displaystyle N^P = \left\{ \begin{array}{ll} 1 & (P=0)\\ \left(N^{P/2}\right)^2 & (P:偶数)\\ N \times N^{P-1} & (P:奇数) \end{array} \right.$

この式を用いて、先ほどの $P=20$ の時についてを計算してみましょう。
$\displaystyle \begin{eqnarray} N^{(20)} &=& \left(N^{(10)}\right)^2\\ &=& \left(\left(N^{(5)}\right)^2\right)^2\\ &=& \left(\left(N^{(1)}\times N^{(4)}\right)^2\right)^2\\ &=& \left(\left(N^{(1)}\times \left(N^{(2)}\right)^2\right)^2\right)^2\\ &=& \left(\left(N^{(1)}\times \left(\left(N^{(1)}\right)^2\right)^2\right)^2\right)^2\\ \end{eqnarray}$

よって、先ほどと同様に $5$ 回の計算で $N^{(20)}$ を求めることが出来ています。

コード(C++)

using ll = long long int;
ll RepeatSquaring(ll N, ll P, ll M){
    if(P==0) return 1;
    if(P%2==0){
        ll t = RepeatSquaring(N, P/2, M);
        return t*t % M;
    }
    return N * RepeatSquaring(N, P-1, M);
}

コードの説明

トップダウンな考え方の節で示した式を、そのまま関数に落とし込んでいるだけです。

具体例

$2^{1000000000} \pmod{1000000007}$ の値を求める。

今回は、愚直な方法(関数PowMod_Simple)と繰り返し自乗法(PowMod_RepeatSquaring)でどのくらい性能の差があるかを、

再帰回数(再帰orループするたびにカウント)
処理時間(std::chronoを用いて計測)

の2点について見てみます。

#include <iostream>
#include <algorithm>
#include <chrono>

using ll = long long int;

ll count_ps = 0;
ll PowMod_Simple(ll N, ll P, ll M){
    ll ans = 1;
    for(ll i=0; i<P; ++i){
        ++count_ps;
        ans *= N;
        ans %= M;
    }
    return ans;
}

ll count_prs = 0;
ll PowMod_RepeatSquaring(ll N, ll P, ll M){
    ++count_prs;
    if(P==0) return 1;
    if(P%2==0){
        ll t = PowMod_RepeatSquaring(N, P/2, M);
        return t*t % M;
    }
    return N * PowMod_RepeatSquaring(N, P-1, M);
}

int main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(0);
    
    //愚直な方法-----------------------------
    std::cout << "愚直な方法\n";
    auto startTime = std::chrono::system_clock::now();
    //処理部------
    std::cout << "2^1000000000 mod 1000000007 = " << PowMod_Simple(2, 1000000000, 1000000007) << "\n";
    std::cout << "再帰回数\t:" << count_ps << "\n";
    //------------
    auto endTime = std::chrono::system_clock::now();
    std::cout << "処理時間:" << std::chrono::duration_cast<std::chrono::nanoseconds>(endTime - startTime).count() << "ns\n";
    //---------------------------------------
    
    std::cout << "------\n";
    
    //愚直な方法-----------------------------
    std::cout << "繰り返し自乗法\n";
    startTime = std::chrono::system_clock::now();
    //処理部------
    std::cout << "2^1000000000 mod 1000000007 = " << PowMod_RepeatSquaring(2, 1000000000, 1000000007) << "\n";
    std::cout << "再帰回数\t:" << count_prs << "\n";
    //------------
    endTime = std::chrono::system_clock::now();
    std::cout << "処理時間:" << std::chrono::duration_cast<std::chrono::nanoseconds>(endTime - startTime).count() << "ns\n";
    //---------------------------------------
    
    return 0;
}

実行結果(一例)

愚直な方法
2^1000000000 mod 1000000007 = 140625001
再帰回数	:1000000000
処理時間	:22322817871ns
------
繰り返し自乗法
2^1000000000 mod 1000000007 = 140625001
再帰回数	:43
処理時間	:2709ns

この結果からも分かるように、再帰回数、処理時間ともに、文字通り桁違いで小さく計算することが出来ていますね。
ちなみに、処理時間の「ns」は「ナノ秒(ナノ $=10^{-9}$ )」を表しているため、

愚直な方法 :22.322817871秒
繰り返し自乗法 : 0.000002709秒

の時間がかかっていることを表しています。

備考

愚直な方法は計算量が $O(P)$ であるのに対し、繰り返し自乗法では $O(\log P)$ にまで抑えられています。

あとがき

繰り返し自乗法は使いどころ(累乗の割った余りを求めたいとき)が分かりやすく、また実装も軽めなので、習得しやすいアルゴリズムだと思います。
一度この関連の問題を解いたら、また同じ場面に出会ったときに「あの再帰のやつだ…！」と思い出すことが出来ると思うので、ATC002Bなどで経験しておくと良いでしょう。

概要

使える状況

説明

クエリ

直線集合に直線を追加する

あるに対し、の中で最小値(あるいは最大値)を取るような直線の値を求める

性質

1.最小値を取る直線は左から順に傾きが減少している

2.不必要な直線が無いとき、あるでの直線の取る値は、その値を直線の傾き順に並べると下に凸となる

不必要な直線を取り除く

特殊な状況における更なる高速化

追加クエリにおける直線の傾きが単調であるとき

最小値クエリにおけるが単調であるとき

上記パターンを両方満たしているとき

コード(C++)

コードの説明

具体例

ソースコード

出力

備考

あとがき

Coding Phase

Easy(300)

Medium(500)

Hard(900)

Intermission

Challenge Phase

最終結果

スコア

レート

あとがき

概要

使える状況

説明

ボトムアップ(小さい値から大きい値へと計算していく)な考え方

トップダウン(大きい値から小さい値へと計算していく)な考え方

コード(C++)

コードの説明

具体例

の値を求める。

実行結果(一例)

備考

あとがき

直線集合 $L$ に直線 $f_k(x)=a_kx+b_k$ を追加する

ある $x$ に対し、 $L$ の中で最小値(あるいは最大値)を取るような直線の値を求める

2.不必要な直線が無いとき、ある $x$ での直線の取る値は、その値を直線の傾き順に並べると下に凸となる

最小値クエリにおける $x$ が単調であるとき

上記 $2$ パターンを両方満たしているとき

$2^{1000000000} \pmod{1000000007}$ の値を求める。