Cork

banbooooのブログ

統計、機械学習、数学に関する記事を投稿します。

ARC094/ABC093 D:Worst Case

競技プログラミング

700点問題のD問題。むずかった。
D - Worst Case
高橋君のとった得点を A,B ( A \leq B) とする。
高橋君よりスコアの小さい最大の人数を求めたいので、1つ目のコンテストが高いスコアならば、2回目は低いスコア(逆もしかり)というように掛け算のペアを作ってあげる。(考察1)
(わからなければ解説放送へ)
www.youtube.com
高橋君よりスコアが小さい人のスコア a,b ( a \leq b) は必ず a < A または b < B を満たすので
まず、 a < A の場合を考えてみる。
(考察1)より、以下のようにペアを考えることができる。
f:id:banboooo:20180408143002j:plain
すると、この場合はすべてのペアの積がABより小さくなることがわかる。
次に b < B の場合を考えてみる。先ほどのように 1 \sim (B-1) を元に、対応する A より大きい数を考えても、必ずABより大きいことが保証される数列が見つからない。
ここで発想を転換して、(A+1) \sim を元に、対応する数を考える。すると、積の最大値がABを超えない範囲で数列の長さを最大にする問題と捉えることができる。
ここまでを整理すると以下の図になる。
f:id:banboooo:20180408152026j:plain
ここからは2つほど解法があると思ったのだが、2つ目の方は実装で挫折してしまった。1つ目は、editorialの方法。2つ目は2分探索と3分探索を組み合わせる方法。
考え方などで間違っていたら教えていただきたいです。

1.editorialの方法
https://beta.atcoder.jp/contests/arc094/submissions/2325009
C^2 < AB を満たす C を取ると、 C(C+1) または C^2 が積の最大値となることを利用する方法。
これは掛け算のペアの和が常に一定ならば
ペアの数が同じくらいの大きさの時に、積の値が最大となるからである。
さらに例外的に、 A = B の時と A+1 = B の時は、 C = A となってしまうので、別で処理しないといけない。
これを踏まえると以下のように考えることができる。


2.三分探索 & 二分探索の方法 (仮)
https://beta.atcoder.jp/contests/arc094/submissions/2325717 <- 不正解になってしまった提出。
C の値をとることに気づかない場合は、こっち。
積の最大値を求めたいため、探索をするが逐次探索では当然間に合わない。
そこで、積の値の関数の形に注目してみると上凸の関数になっているはずである。

こういう時は、三分探索を使えば、極値、今回で言えば積の最大値を非常に効率よく求められる。
極値を Xの関数と見れば、この関数は単調非減少の関数なので、こちらは二分探索でXを定めてあげれば答えがもとまる。

自作ライブラリ(整数編)

競技プログラミング

競プロ用の整数周りのライブラリ(言語はc++14用)
使用可能なメソッド
1.素数関連

  • エラトステネスのふるい
  • 素数かどうか判定
  • 素数列挙
  • n以下の素数の個数

2.約数倍数関連

3.階乗累乗関連
※数が大きいので10^9 + 7を法としてmodをとっている。

  • n階乗
  • a^x
  • mod割り算
  • nCr,nPr
template <typename T>
class Integer {
private:
	const long mod = 1000000007;
	bool primeflg;
	bool factorialflg;
	vector<T> prime_list;
	vector<T> factorial;
	map<T,long> factoring;
	vector<long> prime_num;
	vector<bool> prime_table;

	void factorialModInit(T maxval) {
		factorialflg = true;
		factorial = vector<T>(maxval+1); 
		factorial[0] = factorial[1] = 1;
		for (T i = 2;i < (maxval+1);i++) {
			factorial[i] = (factorial[i-1]*i)%mod;
		}
	}

public:
	Integer()  : primeflg(false),factorialflg(false) {}

	//素数関連
	void Eratosthenes(T n) { // n以下の数でエラトステネスのふるいを作る。 
		primeflg = true;
		prime_table = vector<bool>(n+1,true);
		prime_num = vector<long>(n+1,0);
		T maxiter = sqrt(n);
		for(T i = 2;i < maxiter+1;i++){
			if (prime_table[i]) {
				for (T j = (i + i);j <= n;j += i){
					prime_table[j] = false;
				}
			}
		}
		T prnum = 0;
		for(T i = 2;i < n+1;i++){
			if (prime_table[i]) {
				prnum++;
				prime_list.push_back(i);
			}
			prime_num[i] = prnum;
		}
	}

	bool isPrimeNum(T n) { //素数かどうか判定
		if (!primeflg) Eratosthenes(max(n,(T)100)); 
		return prime_table[n];
	}

	vector<T>* getPrimeList(T n = 100){ //n以下の素数リスト
		if (!primeflg) Eratosthenes(n); 
		return (&prime_list);
	}

	long getPrimeNumUnder_N(T n) { //n以下の素数の個数。
		if (!primeflg) Eratosthenes(max(n,(T)100)); 
		return prime_num[n];
	}

	long getPrimeNum_NM(T k,T m){  //k以上m以下の素数の個数。(kはN以下ならok)
		if (!primeflg) Eratosthenes(max(m,(T)100)); 
		return prime_num[m] - prime_num[k];
	}

	//約数倍数関連
	map<T,long>* getFactoring(T n) { // 素因数分解 key->因数 val->因数の数
		T copy_n = n;
		T maxiter = sqrt(n);
		for (T i = 2;i < maxiter+1;i ++) {
			if (copy_n == 1) break;
			while (!(copy_n % i)) {
				if (IN(i,factoring)) factoring[i]++;
				else factoring[i] = 1;
				copy_n /= i;
			}
		}
		return (&factoring);
	}

	T gcd(T a,T b) {
 		 return b != 0 ? gcd(b, a % b) : a;
	}

	T lcm(T a,T b) {
  		return (a / gcd(a, b))*b;
	}

	long getFactorial(T n) { // nの階乗の値を求める。
		if (!factorialflg) factorialModInit(n);
		return factorial[n];
	}

	long long square_pow(long long a,long x){ // (a^x) % mod (繰り返し二乗)
		T p = a;
		T q = 1LL;
		while (x != 0){
			if (x & 1) q = (q * p) % mod;
			p = (p * p) % mod;
			x >>= 1;
		}	
		return q;
	}

	long long mod_inv(long long a) { // aの逆元 % mod
		return square_pow(a,mod-2);
	}

       long long nCr(T n,T r) { //普通のnCr,modを取らない場合はこっち
		if (r > n / 2)  r = n - r;
		long long retval = 1;
		for (long i = 1;i < r+1;i ++){
			retval *= n - r + i;
			retval /= i;
		}
		return retval;
	}

	long long nCrMod(T n,T r) { // nCr 使用前に初期化必要
		if (!factorialflg) factorialModInit(n);
		return (((factorial[n]%mod * mod_inv(factorial[r]))%mod) * mod_inv(factorial[n-r]))%mod;
	}

	long long nPrMod(T n,T r) {
		if (!factorialflg) factorialModInit(n);
		return ((factorial[n]%mod) * mod_inv(factorial[n-r]))%mod;
	}
};

使用例

int main(void){
	Integer<long> it; //インスタンス生成

	it.Eratosthenes(100L);
	cout << it.isPrimeNum(53) << endl;                  // -> 1 (true)
	cout << it.getPrimeNumUnder_N(50) << endl; // -> 15

	map<long,long>* m = it.getFactoring(12);
	for(auto itr = m->begin(); itr != m->end(); ++itr) { 
		cout << itr->first << " ^ "  << itr->second << endl; 
	} // -> 2 ^ 2 , 3 ^ 1

	cout << it.gcd(18,12) << endl; // -> 6
	cout << it.lcm(18,12) << endl; // -> 36

	cout << it.getFactorial(8) << endl; // -> 40320 (= 8!)
	cout << it.square_pow(2LL,10) << endl; // -> 1024
	cout << it.nCrMod(8,3) << endl; // -> 56
	return 0;
}

Educational Codeforces Round 38 D問題

競技プログラミング

競プロ界ではとても有名(?)なダイクストラ法。
典型問題からちょっとひねった問題で個人的には良問だと思った。
codeforces.com

問題文の要約 (訳に脚色入ってます。)

n個の町とm本の道があって道には通行料が存在。コンサートが全ての町で行われるが、チケット料も町ごとに異なる。町iにいる人が、コンサートを一番安く見る時、その料金(他の町へ行くなら通行料往復分込み)を表示せよ。

アルゴリズム

まずグラフ構造に直す。n頂点m辺、辺の重みは往復分の通行料w_{ij} \times 2とできる。
さて、チケット料の情報をどう扱うかが今回最大の工夫ポイント。
まず、一つの頂点 v_s について考える。通常のダイクストラでは最初優先度付きキューに
(0,始点)をpushするが、今回は(チケット料 a , v_s)をpushする。
ここで下図のような例を考えて、頂点v_s からの最短距離の意味を考える。

このように、頂点v_sから頂点v_tまでの最短距離が、v_tからv_sへコンサートを見に行った時の最小の料金となっている。後は同様に全ての頂点に対して行えば良いのだが、最小の料金 = 最短距離 に対応しているので最初に(チケット料,v)を全ての頂点分pushしてダイクストラ法をやれば distance配列 (通常は始点からの距離を格納する配列)に本問題の答えが格納されている。
計算量は、最初にpushする分だけ増えて O ( V + |E| \log |V|) \approx 10^6
(間違ってたらごめんなさい)

ソースコード

#include <bits/stdc++.h>
#define REP(i,n) for (long i=0;i<(n);i++)
#define FOR(i,a,b) for (long i=(a);i<(b);i++)
#define RREP(i,n) for(long i=n;i>=0;i--)
#define RFOR(i,a,b) for(long i=(a);i>(b);i--)
#define dump1d_arr(array) REP(i,array.size()) cerr << #array << "[" << (i) << "] ==> " << (array[i]) << endl;
#define dump2d_arr(array) REP(i,array.size()) REP(j,array[i].size()) cerr << #array << "[" << (i) << "]" << "[" << (j) << "] ==> " << (array[i][j]) << endl;
#define dump(x)  cerr << #x << " => " << (x) << endl;
#define CLR(vec) { REP(i,vec.size()) vec[i] = 0; } 
#define llINF (long long) 9223372036854775807
#define loINF (long) 2147483647
#define shINF (short) 32767
#define SORT(c) sort((c).begin(),(c).end())
#define MIN(vec) *min_element(vec.begin(), vec.end());
#define MAX(vec) *max_element(vec.begin(), vec.end());
using namespace std;
typedef long long LL;
typedef vector<LL> VL;
typedef vector<VI> VVL;
typedef pair<LL,LL> pr;
struct ver {LL di,node,prev;};
struct Order {
	bool operator() (ver const& a,ver const& b) const {
		return a.di > b.di || ((a.di == b.di) && (a.node > b.node));
	}
};
typedef priority_queue<ver,vector<ver>,Order> pq;

// ダイクストラ法
class Dijekstra {
private:
	vector<vector<pr>> E;
	vector<LL> dist;
	size_t V;
public:
	Dijekstra() : V(0) { }
	Dijekstra(size_t v) :
		V(v),E(v,vector<pr>(0)),dist(v,llINF) {}

	size_t size() { return V;}

	void add_edge_directed(LL from,LL to,LL dis){
		E[from].push_back(pr(to,dis));
	}

	void add_edge_undirected(LL from,LL to,LL dis){
		E[from].push_back(pr(to,dis));
		E[to].push_back(pr(from,dis));
	}	

	void dijekstra(VI weight){
		pq que;
		REP(i,weight.size()) que.push(ver{weight[i],i,i});

		while(! que.empty()){
			ver next = que.top();
			que.pop();
			LL next_v = next.node;
			if (dist[next_v] < next.di) continue; //決定済み
			dist[next_v] = next.di;
			
			REP(i,E[next_v].size()){
				pr e = E[next_v][i];
				if (dist[e.first] > next.di+e.second) {
                                       que.push(ver{next.di+e.second,e.first,next_v});
                                 }
			}
		}
	}

	vector<LL> result(void){
		return dist;
	}

	long result_query(long goal){
		return dist[goal];
	}
};


int main(void){
	long n,m,v,u;
	LL w;
	cin >> n >> m;

	Dijekstra dj(n);

	REP(i,m){
		cin >> u >> v >> w;
		dj.add_edge_undirected(u-1,v-1,w*2); //添え字に注意
	}

	VI weight(n);
	REP(i,n) cin >> weight[i];

	dj.dijekstra(weight);
	VI ans = dj.result();
	REP(i,n) cout << ans[i] << " ";
	cout << "\n";

	return 0;
}

注意

数は結構大きくなるのでlong long型じゃないと通らなかった。

AGC021 ~B問題 Holes ~

競技プログラミング

初投稿です。個人的にはアルゴリズムも実装も大変だったので書き残しておきます。
agc021.contest.atcoder.jp

乱数で決まる(x,y)は半径Rの円の内部。穴は1辺 2 \times 10^6 の正方形の内部の領域に存在する。
Rがめちゃめちゃ大きいので、穴の近くには滅多に落ちない。

全体図

穴を頂点とする凸多角形で、かつ全ての穴を内部含むような図形(凸包というらしい)を作る。乱数で決まる点は主に、この凸多角形の外側の領域に落ちるので、凸多角形の頂点になってる穴のどれかに落ちる。

内部に落ちた時のみ、頂点に使われてない穴にも落ちる可能性があるが、確率は0とみなせる。(精度的に)つまり、凸多角形の頂点に対してのみ確率を考えればよい !!

隣同士の2つの穴について考えると、2つの穴を結んだ直線の垂直二等分線を基準に、どちらの穴に落ちるかがわかるので、これを凸多角形の全ての辺に対してやってみる。

すると、全部の垂直二等分線が中心に集まり、1点で交わるわけではないので小さい領域ができてしまうが、これらを無視すれば、下の図のように領域分割すれば良さそう。

領域分割

中心角 = (垂直二等分線の交わる角度) = \pi- (辺の交わる角度)
例: 赤い領域の中心角 = \pi - \angle BAE

アルゴリズムはここまでだが、実装でも苦労した点が2点ほど。

  • そもそも凸包はどうやって求めるのか。

凸包を求めるアルゴリズムで「グラハムスキャン」というのがあるらしい (蟻本p224)
2辺のベクトル(\vec{v1},\vec{v2})の符号付き角度が、det(v1,v2)の正負に一致するので、これが常に正になるように点を反時計回りに選ぶことがポイント。

  • 長さの値をどうやって角度の値に直すか。

最初は\arctanを使って角度に変換することを考えたが、場合分けが必要でめんどくさい。
そこでベクトルをおいて、\theta = \arccos{\dfrac{\vec{a} \vec{b}}{|a||b|}}で求めると楽だった。

ソースコード
まず二次元ベクトルの内積やノルムなどを簡単に求められるようにクラスを構築しておく。
なお、今回の問題では点ごとに確率を出力しないといけないので、indexの値を保持できるようにしておき、辺のベクトルを表すときは、indexの値は(-1)にしておいた。

typedef struct vec2D {
  double x, y;
  long ind;
  vec2D() {}
  vec2D(double x, double y,long index) : x(x), y(y) ,ind(index) {} 
  vec2D operator + (vec2D p){
	   return vec2D(x + p.x, y + p.y,-1);
  }
  vec2D operator - (vec2D p){
	   return vec2D(x - p.x, y - p.y,-1);
  }
  vec2D operator * (double d){
	   return vec2D(d*x, d*y,-1);
  }

  double norm(void){
  	return sqrt(x*x + y*y);
  }

  double twiceNorm(void){
  	return x*x + y*y;
  }

  double dot(vec2D p){ //pとの内積
	   return x * p.x + y * p.y;
  }
  double det(vec2D p){ // (x,p)の行列式
	   return x * p.y - p.x * y;
  }
  double dist(vec2D p){ //pとの距離の2乗
	   return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y);
  }
  bool operator <(const struct vec2D &e) const{
     return x == e.x? (y < e.y) : x < e.x;
  }
  bool operator >(const struct vec2D &e) const{
     return x == e.x? (y > e.y) : x > e.x;
  }
} pt;

次に、多角形クラスを定義しておく。これでまた凸包を求める問題が出ても怖くない !
なお、グラハムスキャンのところはほぼ蟻本写してます。

bool cmp_x(const pt &p1,const pt &p2){ //ソート用の比較関数
	if (p1.x != p2.x) return p1.x < p2.x;
    return p1.y < p2.y;
}
class Polygon{
  private:
    vector<pt> pts;
    size_t n;
  public:
    Polygon() : n(0) {}
    Polygon(vector<pt> v) : n(v.size()) {
      pts = v;
    }

    void addPoint(pt P){
      pts.push_back(P);
    }

    vector<pt> convexHull(void){ // 凸包を構成する点の配列を返す。
      vector<pt> qs(n * 2);
      sort(pts.begin(),pts.end(),cmp_x);
      long k = 0;

      REP(i,n){ //下側凸包の作成
        while(k > 1 && (qs[k-1] - qs[k-2]).det(pts[i] - qs[k-1]) <= 0) k--;
        qs[k++] = pts[i];
      }

      for(long i = n - 2, t = k; i >= 0; i--){ //上側凸包の作成
        while(k > t && (qs[k-1] - qs[k-2]).det(pts[i] - qs[k-1]) <= 0) k--;
        qs[k++] = pts[i];
      }
      qs.resize(k-1);
      n = k-1;
      pts = qs;
      return qs;
    }
};

最後に、main関数を書く

int main(void){
	short N;
	double x,y;
	cin >> N;
	vector<pt> point,convex;

	REP(i,N) {
		cin >> x >> y;
		point.push_back(pt(x,y,i));
	}

	if (N == 2) { // コーナーケース
		cout << 0.5 << "\n" << 0.5 << endl;
		exit(0);
	}

	Polygon P(point);
	convex = P.convexHull(); //凸包を求める
	vector<double> ans(N,0.0); //答えの確率を格納
	long m = convex.size();

	double alpha;
        // 角度を計算。
	REP(i,m-2){//凸包の点1 ~ (m-2)までの確率を計算
		pt v1 = convex[i]-convex[i+1],v2 = convex[i+2]-convex[i+1];
		alpha = acos((v1).dot(v2) / sqrt(v1.twiceNorm() * v2.twiceNorm()));
		ans[convex[i+1].ind] = (M_PI - alpha)/ (2 * M_PI);
	}
        // 凸包の点m-1の確率
	pt v1 = convex[m-2] - convex[m-1],v2 = convex[0] - convex[m-1];
	alpha = acos((v1).dot(v2) / sqrt(v1.twiceNorm() * v2.twiceNorm()));
	ans[convex[m-1].ind] = (M_PI - alpha)/ (2 * M_PI); 
        // 凸包の点0の確率
	v1 = convex[m-1] - convex[0],v2 = convex[1] - convex[0];
	alpha = acos((v1).dot(v2) / sqrt(v1.twiceNorm() * v2.twiceNorm()));
	ans[convex[0].ind] = (M_PI - alpha)/ (2 * M_PI);

	REP(i,N){
		cout << ans[i] << endl;
	}

	return 0;
}