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

競プロ用の整数周りのライブラリ(言語は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;
}