torus711 のアレ

主に競技プログラミングの問題について書きます.PC 以外だと数式が表示されないかもしれないです

GCDLCM2

問題概要

 正整数の配列が与えられる.この配列に対し,以下の操作を任意回行う.

  1. 2 つの要素 $x, y$ を選ぶ
  2. $x, y$ を削除する
  3. $\mathit{ GCD }( x, y ), \mathit{ LCM }( x, y )$ (それぞれ最大公約数と最小公倍数)を追加する

 配列の要素の和として達成できるものの最大数を $S$ として,$S \bmod{ 10^9+7 }$ を求めよ.

 入力となる配列は,生成のためのアルゴリズムとそれに与えるパラメータとして与えられる.完成した配列の要素数は $10^5$ 以下,各要素は $10^7$ 以下となる.

解法

 まずは一回の操作で何が起こるか考えていきます.$x = y$ のとき $x = y = \mathit{ GCD }( x, y ) = \mathit{ LCM( x, y ) }$ となって変化しないので,以降 $x \neq y$ として,ついでに面倒なので $x < y$ とします(一般性は損なわない).また,手順 3 で加えられる要素を $x' = \mathit{ GCD }( x, y ), y' = \mathit{ LCM }( x, y )$ とします.$\mathit{ LCM }( x, y ) = \frac{ xy }{ \mathit{ GCD }( x, y )}$ なので,変化はそれぞれ \[
x' = x \frac{ \mathit{ GCD }( x, y ) }{ x } \\
y' = y \frac{ x }{ \mathit{ GCD }( x, y ) }
\] と書け,$\frac{ x }{ \mathit{ GCD }( x, y ) }$ を $y$ に掛けて,$x$ を割るという操作になっています.$\frac{ x }{ \mathit{ GCD }( x, y ) } > 1$ であることに注意すると,$y \mapsto y'$ による値の増分が,$x \mapsto x'$ による減少分を上回ります.このことから,有効な操作が可能ならばできるだけやった方が得になると分かります.
 $\mathit{ LCM }$ を取る操作は,2 数を素因数分解した各素因数について,指数の大きい方を掛け合わせて数を作る操作です.可能な限り操作をした時に作られる最大の要素は,配列中の全ての要素について,各素因数の最大の指数を採用してできる値となります.この操作を要素数分やって,作られる各値の和(の剰余)を計算すれば問題を解くことができます.
 従って,与えられた配列の各要素を素因数分解し,素数ごとに指数を求めてソートした列を用意すれば実装できます.

 計算量についてです.与えられた配列の要素数を $N$ とし,そこに含まれる最大の要素を $M $ とします.
 まず,各要素を素因数分解するために $O( N \sqrt M )$ 時間かかります.一つの要素が含む素因数の種類(=出てくる指数の数)は $O( \log M )$ と言えて,分類後のソートに $O( N \log M \log( N \log M ) )$ 時間かかります.大きい指数を順番に取り出すときに使い終わった素数を残していると $O( \frac{ NM }{ \log_{ \mathrm{ e } } M } )$ (素数定理より)となってやばいですが,使い終わった(指数が $0$ しか残ってない)素数を消していけば $O( N \log M )$ にできます.
 結局,全体では $O( N \sqrt M + N \log M \log( N \log M ) )$ となります.

コード

using LL = long long;
using VI = vector< int >;
using VVI = vector< vector< int > >;
using PII = pair< int, int >;

#define REP2( i, n ) REP3( i, 0, n )
#define REP3( i, m, n ) for ( int i = ( int )( m ); i < ( int )( n ); ++i )
#define GET_REP( a, b, c, F, ... ) F
#define REP( ... ) GET_REP( __VA_ARGS__, REP3, REP2 )( __VA_ARGS__ )
#define ALL( c ) begin( c ), end( c )

#define SZ( v ) ( (int)( v ).size() )

#define PB push_back

#define fst first
#define snd second

constexpr int MOD = 1000000007;

// 素数かどうかの表を計算
// Eratosthenes の篩 : O( N log log N )
vector<bool> eratosthenes( const int N );

// N 以下の素数のリストを生成する
// include : eratosthenes
vector<int> getPrimes( const int N );

// 素因数分解 O( √N )
vector< long long > primeFactorization( long long N );

// a^x を mod で求める
// 反復二乗法
// O( log x )
long long mod_pow( long long a, long long x, long long mod );

class GCDLCM2
{
public:
	int getMaximalSum( vector <int> start, vector <int> d, vector <int> cnt )
	{
		const int L = SZ( start );
		const int N = accumulate( ALL( cnt ), 0 );

		VI primes = getPrimes( 10000000 );

		VVI exps( SZ( primes ) );
		REP( i, L )
		{
			REP( j, cnt[i] )
			{
				const int a = start[i] + j * d[i];
				const auto factors = primeFactorization( a );

				map< int, int > counts;
				for_each( ALL( factors ), [&]( const LL p ){ ++counts[p]; } );
				for_each( ALL( counts ), [&]( const PII &c ){ exps[ lower_bound( ALL( primes ), c.fst ) - begin( primes ) ].PB( c.snd ); } );
			}
		}
		for_each( ALL( exps ), []( VI &row ){ sort( ALL( row ) ); } );

		LL res = 0;
		REP( i, N )
		{
			LL a = 1;
			VI nprimes;
			VVI nexps;

			REP( j, SZ( primes ) )
			{
				if ( exps[j].empty() )
				{
					continue;
				}

				( a *= mod_pow( primes[j], exps[j].back(), MOD ) ) %= MOD;
				exps[j].pop_back();

				nprimes.PB( primes[j] );
				nexps.PB( move( exps[j] ) );
			}
			( res += a ) %= MOD;

			swap( primes, nprimes );
			swap( exps, nexps );
		}

		return res;
	}
};