好题！！！没话说……
用到的知识面很多，这题的难点在于公式的推导。
原始推导过程见：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 ",sum%(mod/n)); 137 } 138 return 0; 139 }