hdu 2481 Toy

好题！！！没话说……

用到的知识面很多，这题的难点在于公式的推导。

原始推导过程见：http://hi.baidu.com/spellbreaker/item/d8bb3bda5af30be6795daa93

这个过程有点小问题，改进后的见：http://blog.csdn.net/acm_cxlove/article/details/7868589

下面是自己敲的代码：

 

  1 /*************************************************************************
  2     > File Name: xh.cpp
  3     > Author: XINHUA
  4     > Mail: 525799145@qq.com 
  5     > Created Time: 2013/7/19 星期五 15:11:25 新华
  6  ************************************************************************/
  7 
  8 #include<iostream>
  9 #include<stdio.h>
 10 #include<algorithm>
 11 #include<iomanip>
 12 #include<cmath>
 13 #include<string>
 14 #include<vector>
 15 using namespace std;
 16 int prime[30001],m,n;
 17 bool f[30001]={0};
 18 __int64 mod;
 19 struct ma
 20 {
 21     __int64 a[2][2];
 22     void init()
 23     {
 24         a[0][0]=a[1][1]=1;
 25         a[0][1]=a[1][0]=0;
 26     }
 27 }e;
 28 //使用二分模拟乘法计算（a和b的范围太大）
 29 __int64 mulmod(__int64 a,__int64 b)
 30 {
 31     a%=mod;b%=mod;
 32     if(a<0) a+=mod;
 33     if(b<0) b+=mod;
 34     __int64 ans=0;
 35     while(b)
 36     {
 37         if(b&1)
 38         {
 39             ans+=a;
 40             if(ans>=mod) ans-=mod;
 41         }
 42         a=a<<1;
 43         if(a>=mod) a-=mod;
 44         b=b>>1;
 45     }
 46     return  ans%mod;
 47 }
 48 ma operator*(ma ab,ma b)
 49 {
 50     ma ans;
 51     int i,j,k;
 52     for(i=0;i<2;i++)
 53     {
 54         for(j=0;j<2;j++)
 55         {
 56             ans.a[i][j]=0;
 57             for(k=0;k<2;k++)
 58                 ans.a[i][j]=(ans.a[i][j]+mulmod(ab.a[i][k],b.a[k][j]))%mod;
 59         }
 60     }
 61     return ans;
 62 }
 63 ma operator^(ma ab,int b)
 64 {
 65     ma ans;
 66     ans.init();
 67     while(b)
 68     {
 69         if(b&1) ans=ans*ab;
 70         b>>=1;
 71         ab=ab*ab;
 72     }
 73     return ans;
 74 }
 75 
 76 void init_prime()
 77 {
 78     int i,j;
 79     m=0;
 80     for(i=2;i<=30000;i++)
 81     {
 82         if(f[i]==0) prime[m++]=i;
 83         for(j=0;j<m&&i*prime[j]<=30000;j++)
 84         {
 85             f[i*prime[j]]=1;
 86             if(i%prime[j]==0) break;
 87         }
 88     }
 89 }
 90 __int64 euler(__int64 n)
 91 {
 92     __int64 ans=1;
 93     for(int i=0;i<m&&prime[i]*prime[i]<=n;i++)
 94     {
 95         if(n%prime[i]==0)
 96         {
 97             ans*=prime[i]-1;
 98             n/=prime[i];
 99             while(n%prime[i]==0)
100             {
101                 n/=prime[i];
102                 ans=(ans*prime[i])%mod;
103             }
104         }
105     }
106     if(n>1) ans*=n-1;
107     return ans%mod;
108 }
109 __int64 get_T(int k)
110 {
111     if(k==1) return 1;
112     else if(k==2) return 5;
113     ma temp=e^(k-2);
114     __int64 f=3*temp.a[0][0]+temp.a[1][0];
115     __int64 g=2*(f-(3*temp.a[0][1]+temp.a[1][1])-1);
116     return (f+g)%mod;
117 }
118 int main()
119 {
120     init_prime();
121     e.a[0][0]=3;e.a[0][1]=1;e.a[1][0]=-1;e.a[1][1]=0;
122     while(scanf("%d%I64d",&n,&mod)!=EOF)
123     {
124         mod=(__int64)n*mod;
125         __int64 sum=0;
126         for(int i=1;i*i<=n;i++)
127         {
128             if(n%i==0)
129             {
130                 sum=(sum+mulmod(euler(i),get_T(n/i)))%mod;
131                 if(i*i!=n)
132                     sum=(sum+mulmod(euler(n/i),get_T(i)))%mod;
133             }
134         }
135         sum/=n;
136         printf("%I64d\n",sum%(mod/n));
137     }
138     return 0;
139 }
View Code

 

 

 

posted @ 2013-07-20 09:03  _随心所欲_  阅读(315)  评论(0)  编辑  收藏  举报