今週のアルゴリズム：最大公約数と最小公倍数

問題：https://codeiq.jp/ace/thisweek_masuipeo/q665

- 全探索
- 全然時間が間に合いませんでした。

```cpp:codeiq665.cpp
#include <vector>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
using namespace std;
int N;
int gcd(int x,int y){return y?gcd(y,x%y):x;}
bool ok(int x,int y){return min(x,y)==gcd(x,y);}
bool frozen(vector<int> &v){ //全てのi,jについて、min(v[i],v[j])はmax(v[i],v[j])の約数
	for(int i=0;i<N-1;i++){
		for(int j=i+1;j<N;j++)if(!ok(v[i],v[j]))return false;
	}
	return true;
}
void output(vector<int> &v){
	for(int i=0;i<N;i++)printf("%d ",v[i]);
	puts("");
}
int m;
void dfs(vector<int> &v,int d){
	if(d>=m)return;
	if(frozen(v)){
		fprintf(stderr,"Found %d\n",d); //26
		fflush(stderr);
		m=d;
	}
	for(int i=2;i<N-1;i++){ //2は20の約数であるから無視して問題ない
		if(!d&&i>2)return;
		for(int j=i+1;j<N;j++){
			int g=gcd(v[i],v[j]);
			if(min(v[i],v[j])!=g){
				int x=v[i],y=v[j];
				v[i]=g,v[j]=x/g*y;
				dfs(v,d+1);
				v[i]=x,v[j]=y;
			}
		}
	}
}
int main(int argc, char **argv){
	if(argc<2){
		printf("give me array size\n");
		return 1;
	}
	N=strtol(argv[1],NULL,10);
	m=N; ///
	int i,j,n;
	vector<int>v(N);
	for(i=0;i<N;i++){
		v[i]=i+1;
	}
	dfs(v,0);
	printf("%d\n",m);
}
```

- モンテカルロ法
- 上記のdfsで愚直に計算していた部分について、変更する座標をmt19937_64でランダムに選ぶようにした。
- 確かめは20億回シミュレーションして行ったが、実際には1億回程度で12を出せるらしい。

```cpp:codeiq665_rand.cpp
#include <vector>
#include <random>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <unistd.h>
using namespace std;
int N;
int gcd(int x,int y){return y?gcd(y,x%y):x;}
bool ok(int x,int y){return min(x,y)==gcd(x,y);}
bool frozen(vector<int> &v){ //全てのi,jについて、min(v[i],v[j])はmax(v[i],v[j])の約数
	for(int i=0;i<N-1;i++){
		for(int j=i+1;j<N;j++)if(!ok(v[i],v[j]))return false;
	}
	return true;
}
void output(vector<int> &v){
	for(int i=0;i<N;i++)printf("%d ",v[i]);
	puts("");
}
int main(int argc, char **argv){
	if(argc<2){
		printf("give me array size\n");
		return 1;
	}
	N=strtol(argv[1],NULL,10);
	uniform_int_distribution<int> distribution(2,N-1);
	mt19937_64 engine((unsigned int)time(NULL)^(getpid()<<16));
	int m=N*N;
	int t,i,n;
	for(t=0;t<100000000;t++){ //monte carlo simulation (100 million times)
		vector<int>v(N);
		for(i=0;i<N;i++){
			v[i]=i+1;
		}
		for(n=0;n<m&&!frozen(v);n++){ //nは高々m
			for(;;){
				int a=distribution(engine);
				int b=distribution(engine);
				if(!ok(v[a],v[b])){
					int g=gcd(v[a],v[b]);
					v[b]*=v[a]/g;
					v[a]=g;
					break;
				}
			}
		}
		if(m>n)printf("Found %d in try %d\n",n,t);
		m=n;
	}
	printf("%d\n",m); //12回
}
/*
possibly:
0 0 1 1 2 2  3  4  4  5
6 6 7 8 9 9 10 11 12 12

Found 27 in try 0
Found 24 in try 1
Found 21 in try 2
Found 18 in try 4
Found 17 in try 142
Found 16 in try 217
Found 15 in try 4181
Found 14 in try 30536
Found 13 in try 644573
Found 12 in try 760597
*/
```

- 高速版
- 全ての座標について処理した後、N以下の数値を全て取り出して重複数を数えれば良いそうです。

```cpp:665_fast.cpp
#include <vector>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
using namespace std;
int N;
int gcd(int x,int y){return y?gcd(y,x%y):x;}
int main(int argc, char **argv){
	if(argc<2){
		printf("give me array size\n");
		return 1;
	}
	N=strtol(argv[1],NULL,10);
	int i,j;
	vector<int>v(N);
	for(i=0;i<N;i++){
		v[i]=i+1;
	}
	for(i=2;i<N-1;i++){ //N=3の場合はi=1にする必要がある
		for(j=i+1;j<N;j++){
			int g=gcd(v[i],v[j]);
			v[j]*=v[i]/g;
			v[i]=g;
		}
	}
	vector<int>ret;
	for(i=0;i<N;i++)if(v[i]<=N)ret.push_back(v[i]);
	sort(ret.begin(),ret.end());
	vector<int>::iterator it=unique(ret.begin(),ret.end());
	for(j=0;it!=ret.end();++it)j++;
	printf("%d\n",j); //12回
}
```