Hatena::Grouptopcoder

kuuso1の日記

2015-02-25

Ad Infinitum 10

| 00:54 | Ad Infinitum 10 - kuuso1の日記 を含むブックマーク はてなブックマーク - Ad Infinitum 10 - kuuso1の日記

F: Fibonacci LCM

  • N 個の整数 A[n] 1 ≤ n ≤ N が与えられる
  • F1 = 1, F2 = 1, F(n) = F(n-1) + F(n-2) ( n ≥ 2 ) で与えられるフィボナッチ数列について、LCM( F(A[1]), F(A[2]),…, F(A[N]) )をmod(1e9+7) で求めよ

(制約)1≤ N ≤100, 1≤ A[n] ≤ 10^9

・半日くらい考えて、ググってみたらこのページをみつけてファッ?!ってなった。

 ポイントは GCM( F(n), F(m) ) = F( GCM(n,m) ) らしい。ほんとかよ。

・あとはLCMを求めるだけだが、A[n]が大きすぎるので直接 F(A[n])とかは求められないので、例によって行列累乗でmod(1e9+7)で求めることを考える

・N=2 の場合は 上の式から LCM( F(n), F(m) ) = F(n)*F(m)/GCM( F(n), F(m) ) = F(n)*F(m)/F( GCM(n,m) ) で求まるが、N ≥ 3 だと 途中でmodをとってしまうと順次計算では求められない。

・なので、包除原理を使う。

f:id:kuuso1:20150226005103p:image

・配るDPで高次のgcdを計算していくが、基本的には互いに素な因子がどんどん消えていくので、項数が増えるのは同じ因子を持つ数が多いケース。

・そこで、因子に対して個数を指数で管理すれば、GCDの組み合わせの項数が増えない代わりに、指数が組み合わせ的に大きくなる。

・が、指数累乗するのである程度大きくても問題ないうえ、mod(1e9+7)で非零のベキなので、フェルマーの小定理より指数はmod(1e9+6)で管理すればよい。

using System;
using System.Collections;
using System.Collections.Generic;
 
class TEST{
	static void Main(){
		Sol mySol =new Sol();
		mySol.Solve();
	}
}

class Sol{
	public void Solve(){
		
		mod=(long)1e9+7;
		
		Dictionary<long,long>[] D=new Dictionary<long,long>[N];
		for(int i=0;i<N;i++)D[i]=new Dictionary<long,long>();
		
		for(int i=0;i<N;i++){
			for(int j=N-1;j>0;j--){
				foreach(var k in D[j-1].Keys){
					if(gcd(k,A[i])==1)continue;
					if(!D[j].ContainsKey(gcd(k,A[i])))D[j].Add(gcd(k,A[i]),0);
					D[j][gcd(k,A[i])]+=D[j-1][k];
					D[j][gcd(k,A[i])]%=(mod-1);
				}
			}
			if(A[i]==1)continue;
			if(!D[0].ContainsKey(A[i]))D[0].Add(A[i],0);
			D[0][A[i]]+=1;
		}
		
		long Ans=1;
		for(int i=0;i<N;i++){
			Ans*=fib(A[i]);
			Ans%=mod;
		}

		int div=1;
		for(int i=1;i<N;i++){
			foreach(var k in D[i].Keys){
				long t=fib(k);
				if(div==1){
					t=modComb.ModInv(t);
				}
				Ans*=modComb.ModPow(t,D[i][k]);
				Ans%=mod;
			}
			div^=1;
		}
		
		Console.WriteLine(Ans);
		
	}
	int N;
	long[] A;
	public Sol(){
		N=ri();
		A=new long[N];
		for(int i=0;i<N;i++){
			A[i]=rl();
		}
	}
	long mod;
	
	long gcd(long a,long b){
		return a==0?b:gcd(b%a,a);
	}
	
	
	long fib(long x){
		long a=1;
		long b=1;
		long n=x;
		Mat[] MM=new Mat[40];
		MM[0]=new Mat(1,0,0,1);
		MM[1]=new Mat(1,1,1,0);
		for(int i=2;i<40;i++){
			MM[i]=(MM[i-1]*MM[i-1])%mod;
		}
		if(n==1){
			return a;
		}
		if(n==2){
			return b;
		}
		
		Mat X=new Mat(1,0,0,1);
		n-=2;
		for(int i=1;i<40;i++){
			if(((n>>(i-1))&1)>0){
				X=MM[i]*X;
				X%=mod;
			}
		}
		long ans=((X.a*b)%mod+(X.b*a)%mod)%mod;
		return ans;
	}
	
	class Mat{
		public long a,b,c,d;
		public Mat(long x,long y,long z,long w){
			a=x;b=y;c=z;d=w;
		}
		
		public static Mat operator*(Mat x,Mat y){
			return new Mat(x.a*y.a+x.b*y.c,
					x.a*y.b+x.b*y.d,
					x.c*y.a+x.d*y.c,
					x.c*y.b+x.d*y.d);
		}
		public static Mat operator%(Mat x,long m){
			return new Mat(x.a%m,x.b%m,x.c%m,x.d%m);
		}
		public override String ToString(){
			return (a+" "+b+" "+c+" "+d).ToString();
		}
	}

	static int ri(){return int.Parse(Console.ReadLine());}
	static long rl(){return long.Parse(Console.ReadLine());}
}

class modComb{
	public static long mod=(long)(1e9+7);

	public static long ModPow(long n,long k){
		if(n>mod)n%=mod;
		if(k==0)return 1;
		if(n==0)return 0;
		List<long> L=new List<long>();
		
		int idx=0;
		while(mod>(0x1<<idx))idx++;
		long[] Table=new long[idx];
		Table[0]=n;
		for(int i=1;i<idx;i++){
			Table[i]=Table[i-1]*Table[i-1];
			Table[i]%=mod;
		}
		
		long ret=1;
		for(int i=0;i<idx;i++){
			if(( k&(0x1<<i)) >0){
				ret*=Table[i];
				ret%=mod;
			}
		}
		
		return ret;
		
	}
	public static long ModInv(long n){
		//(mod-2)乗を返す
		return ModPow(n,mod-2);
	}
}

・終了30分前にAC。久しぶりにACでガッツポーズでた。