SPFを利用した素因数分解

spf_fact

ABC179のC問題で約数の個数を出していました。この時に”複数クエリの素因数分解に対して”高速な素因数分解アルゴリズムが必要になったので、メモとしてまとめておきます。

1.ABC179C問題

$A*B+C=N$ となる $A,B,C$ の組の通りを出力せよという問題です。解説では最終的に $A$ を全探索で解いています。しかし、私は $A*B=N-C$ と式変形し、 $N'=N-C$ とすれば、 $A*B=N'$ となるので、 $N'$ の約数の個数の総和を求めればいいと考えました。

2.約数列挙

AtCoder 版！マスター・オブ・整数 (素因数分解編)を参考にしました。

	#include <bits/stdc++.h>
	#define rep(i, n) for (int i = 0; i < (n); i++)
	#define rep2(i, s, n) for (int i = (s); i < (n); i++)
	using namespace std;
	using ll = long long;
	using P = pair<int, int>;

	template <typename T>
	vector<T> divisor(T n) {
	vector<T> res;
	for (T i = 1; i * i <= n; i++) {
	if (n % i == 0){
	res.push_back(i);
	if (n / i != i) res.push_back(n / i);
	}
	}
	sort(res.begin(), res.end());
	return res;
	}

	int main() {
	int n;
	cin >> n;
	auto result = divisor(n);
	for (auto res : result) cout << res << " ";
	cout << endl;
	return 0;
	}

view raw basic_divisor.cpp hosted with

by GitHub

約数を探すときは $\sqrt{n}$ まで探せば良くて、それと対になる数字も一緒に探しています。

そのため、計算量は $O(\sqrt{N})$ になります。これでも十分高速です。

3.SPFを利用したアルゴリズム

上の方法で書いた場合、全体の計算量は $O(N\sqrt{N})$ になります。間に合いそうですが、実装のやり方によってはTLEします。私はTLEになりました。if文がボトルネックになっているらしいです。

さて、それ以外で高速に約数を求める必要があります。そういうわけで、SPFを利用した素因数分解を用います。

3-1.SPF

SPFとはSmallest Prime Factorの略で、1以外の最小の素因数のことです。例えば10なら、SPFは2になります。

3-2.SPFの求め方

エラトステネスの篩と同様に求めていきます。

template <typename T>
vector<T> smallest_prime_factor(T n) {
  vector<T> spf(n + 1);
  for (T i = 0; i <= n; i++) spf[i] = i;
  for (T i = 2; i * i <= n; i++) {
    if (spf[i] == i) {
      for (T j = i * i; j <= n; j += i) {
        if (spf[j] == j) spf[j] = i;
      }
    }
  }
  return spf;
}
Copy

n=10ぐらいでやってみると、ちゃんと篩されていることがよくわかります。

3-3.素因数分解

template <typename T>
vector<T> factolization(T x, vector<T> &spf) {
  vector<T> v;
  while (x != 1) {
    v.push_back(spf[x]);
    x /= spf[x];
  }
  sort(v.begin(), v.end());
  return v;
}
Copy

SPFを先に計算されているとすると、素因数分解が簡単にできます。例えばx=48であった場合、

vにspf[48]=2をpush_back
x=48/2=24
vにspf[24]=2をpush_back
x=24/2=12
vにspf[12]=2をpush_back
x=12/2=6
vにspf[6]=2をpush_back
x=6/2=3
vにspf[3]=3をpush_back
x=3/3=1
x==1なので終了

という処理をします。

これで素因数分解をすると、 $2,2,2,2,3$ と出てきますが、約数の個数を求めるときは、各素因数の個数が必要になります。しかし、このvector内をもう一回見ていくのはめんどくさいです。そのため、mapで管理するとうまくいきます。

template <typename T>
map<T, T> prime_factolization(T x, vector<T> &spf) {
  map<T, T> mp;
  while (x != 1) {
    mp[spf[x]]++;
    x /= spf[x];
  }
  return mp;
}
Copy

上のコードと殆ど一緒です。ただ、mapで管理するようにしただけです。

これによって、n=48としたとき、2 4,3 1と出力されます。前が素因数で後ろがその素因数の個数です。

3-4.コード全体

	#include <bits/stdc++.h>
	#define rep(i, n) for (int i = 0; i < (n); i++)
	#define rep2(i, s, n) for (int i = (s); i < (n); i++)
	using namespace std;
	using ll = long long;
	using P = pair<int, int>;

	template <typename T>
	vector<T> smallest_prime_factor(T n) {
	vector<T> spf(n + 1);
	for (T i = 0; i <= n; i++) spf[i] = i;
	for (T i = 2; i * i <= n; i++) {
	if (spf[i] == i) {
	for (T j = i * i; j <= n; j += i) {
	if (spf[j] == j) spf[j] = i;
	}
	}
	}
	return spf;
	}

	// x = a * b * c * ...のように表示する時
	template <typename T>
	vector<T> factolization(T x, vector<T> &spf) {
	vector<T> v;
	while (x != 1) {
	v.push_back(spf[x]);
	x /= spf[x];
	}
	sort(v.begin(), v.end());
	return v;
	}

	// x = a^i * b^j * ...のように表示する時
	template <typename T>
	map<T, T> prime_factolization(T x, vector<T> &spf) {
	map<T, T> mp;
	while (x != 1) {
	mp[spf[x]]++;
	x /= spf[x];
	}
	return mp;
	}

	int main() {
	constexpr int max = 1000000;
	auto spf = smallest_prime_factor(max);

	int n;
	cin >> n;
	auto fact = factolization(n, spf);
	auto p_fact = prime_factolization(n, spf);

	for (auto num : fact) cout << num << " ";
	cout << endl;

	for (auto num : p_fact) cout << num.first << " " << num.second << endl;

	return 0;
	}

view raw spf_factolization.cpp hosted with

by GitHub

main関数の最初の2行で使うであろう範囲のspfを求めています。ここをどうにかしてmain関数外に出したいなあと思ったのですが、やり方が分からず。

4.結果

ABC179のC問題で通したコードは次のとおりでした。

	#include <bits/stdc++.h>
	#define rep(i, n) for (int i = 0; i < (n); i++)
	#define rep2(i, s, n) for (int i = (s); i < (n); i++)
	using namespace std;
	using ll = long long;
	using P = pair<int, int>;

	template <typename T>
	vector<T> smallest_prime_factors(T n) {
	vector<T> spf(n + 1);
	for (int i = 0; i <= n; i++) spf[i] = i;

	for (T i = 2; i * i <= n; i++) {
	// 素数だったら
	if (spf[i] == i) {
	for (T j = i * i; j <= n; j += i) {
	// iを持つ整数かつまだ素数が決まっていないなら
	if (spf[j] == j) {
	spf[j] = i;
	}
	}
	}
	}

	return spf;
	}

	template <typename T>
	vector<T> factolization(T x, vector<T> &spf) {
	vector<T> ret;
	while (x != 1) {
	ret.push_back(spf[x]);
	x /= spf[x];
	}
	sort(ret.begin(), ret.end());
	return ret;
	}

	int main() {
	constexpr int MAX = 1000000;

	auto spf = smallest_prime_factors(MAX);

	int n;
	cin >> n;

	int ans = 0;
	for (int i = 1; i < n; i++) {
	int count = 1;
	map<int, int> mp;
	auto result = factolization(n - i, spf);
	for (auto z : result) mp[z]++;
	for (auto itr = mp.begin(); itr != mp.end(); ++itr) {
	count *= itr->second + 1;
	}
	ans += count;
	}
	cout << ans << endl;
	return 0;
	}

view raw ABC179c.cpp hosted with

by GitHub

普通に約数列挙をした場合はTLEでしたが、このコードの場合は209msと1/10まで実行時間を短縮できます。

素因数分解はよく使うので、複数クエリの場合はこれを利用したいと思います。もちろん、一つの整数に対して素因数分解するならば、通常の方法の方が高速です。

また、spfを利用する場合は前処理が重いので最大値が1e12など巨大な数になる場合は使えません。せいぜいがintに収まるぐらいの整数のときです。

5.参考文献

AtCoder 版！マスター・オブ・整数 (素因数分解編)
複数のクエリの素因数分解を高速に処理する ←コードの大部分を参考にさせていただきました。
素因数分解のアルゴリズム

SPFを利用した素因数分解

1.ABC179C問題

2.約数列挙

3.SPFを利用したアルゴリズム

3-1.SPF

3-2.SPFの求め方

3-3.素因数分解

3-4.コード全体

4.結果

5.参考文献

関連記事