読者です 読者をやめる 読者になる 読者になる

競技プログラミングをするんだよ

ICPC国内予選突破を目標に一日一問題以上解いていきます。

UVa 106 - Fermat vs. Pythagoras

UVa

C++でinlineというテクニックを覚えました。
再帰をしないが何度も使用する関数に用いるとメモリの許す限り高速化が図れるということなので意識していきたいと思います。
繰り返しのマクロ文が何が適切なのか定まらなくてもやもやします。
今まではrep(i,a,b) でa以上b未満としてきたのですが今回は問題の性質上a以上b以下もほしかったので、E=equal or less than,T=less thanとしました。

- 問題概要

フェルマーの最終定理というものが存在し
自然数a,b,cについてa^N + b^N = c^N(N > 2)を満たすものは存在しないというものです。
これはしっかりと証明もされております。
ここまでは問題背景のお話。

ある整数Nが入力で与えられます
0 < a < b < c \le N,a^2 + b^2 = c^2を満たす自然数の組み合わせ(a,b,c)について、

  1. a,b,cの値が互いに素である組み合わせの総数
  2. N以下で全ての組み合わせに含まれない自然数の個数

をそれぞれ出力せよ、という問題です。

- 解法

0 < a < b < Nを満たす全ての(a,b)について探索をしようとすると計算量はO(N^2)となります。
この解法はNの最大が10^6なので間に合わずTLEとなります。
ここで、x < yを満たす自然数x,yを用いて、a,b,cを以下のように表します。
a=y^2-x^2\\
b=2xy\\
c=y^2+x^2
すると、a^2 + b^2 = c^2が成り立ちます。詳しい解説は別のサイトを参照してください。
これを利用して問題を解きます。
整数Mをm^2 < Nを満たす最大のmと置くと、0 < x < y \le Mを満たす整数(x,y)について上記の式を使って三平方を求めます。
ここでcの値がNを超えるものは捨てる必要があることに注意が必要です。
また、必ずしもa < bが成り立つわけではないので場合分けしてスワップを行う必要があります。
この方法で列挙した(a,b,c)の組み合わせの集合Uは、N以下の互いに素な(a,b,c)の組み合わせをすべて含みます。
しかしUは互いに素でない(a,b,c)の組み合わせもすべて含むというわけではないので注意が必要です。
M = \sqrt{N}が成り立つので全てのx,yの列挙に必要な計算量はO(\sqrt{N}^2)=O(N)です。
見つかった(a,b,c)の組み合わせに対して、c \le kc \le Nを満たす(ka,kb,kc)の組み合わせも三平方を満たすので繰り返し文を使って列挙しましょう。
ka,kb,kcについて配列を用いてその値が出たかどうかを更新します。
ある整数が使われたかどうかの管理を二分探索木で行うと計算量が大きくなるようです。
この繰り返しはO(N)であるように思えますが全体の繰り返し回数の総計は10^7以下なので間に合います。
これに関しては計算ができなかったので実測値からの結果論です。whileの中にカウンタを入れてカウントしました。
また、(a,b,c)の組み合わせが互いに素であるかどうかは(a,b)(b,c)(c,a)の全ての組み合わせに対してユークリッドの互助法を用いて最大公約数が1であるかどうか調べればよいです。
ユークリッドの互除法O(\log N)で最大公約数を計算できます。
また、互いに素な組み合わせの管理は二分探索木を使って保持しましょう。今回はsetを利用しました。
setを使うことにより要素の追加もO(\log N)です。
したがって、計算量はO(N\log N)となります。
N以下で三平方に使われていない総数の計算は全ての計算が終わった後に配列を全て参照して使われていない数字の個数をカウントしましょう。

- ソースコード

#include<iostream>
#include<algorithm>
#include<set>
#include<utility>
#define repE(i,a,b) for(auto (i)=(a);(i)<=(b);(i)++)
#define repT(i,a,b) for(auto (i)=(a);(i)<(b);(i)++)
#define itr(i,b) for(auto (i)=(b).begin();(i)!=(b).end();++(i))
#define MAX_N 1000000
#define MAX_LOGN 1000
using namespace std;

int N, M, a, b, c;
int used[MAX_N + 1];
int db[MAX_LOGN + 1];
set<pair<int, int> > s;
int gcd(int x, int y){
	return(y == 0 ? x : gcd(y,x%y));
}

inline void init(){
	repE(i, 1, MAX_LOGN)
		db[i] = i*i;
}

inline void clear(){
	s.clear();
	repE(i, 1, N)
		used[i] = 1;
	M = lower_bound(db, db + MAX_LOGN, N)-db-1;
}


inline int check(){
	return (gcd(a, b) == 1 && gcd(b, c) == 1 && gcd(c, a) == 1);
}
int main(void){
	init();
	while (cin >> N){
		clear();
//		cout << M << endl;
		repE(i,2,M)
			repT(j, 1, i){
			a = db[i] - db[j];
			b = 2 * i * j;
			c = db[i] + db[j];
			if (c <= N&&check()){
				int k = 1;
				if (a > b)swap(a, b);
				auto tmp = make_pair(a, b);
				s.insert(tmp);
				while (c*k <= N){
					used[a*k] = used[b*k] = used[c*k] = 0;
					k++;
				}
			}
		}
		M = 0;

		repE(i, 1, N){
			M += used[i];
		}
		cout << s.size() <<" "<<M<< endl;
	}
}