题意：
Given n (1 <= n <= 1018), You should solve for
where
得到1E9+7的循环节是222222224,222222224的循环节是183120.
然后三次矩阵快速幂就行了。
需要注意每次都要判断那一层的n是否为0和1。
暴力解法：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 |
/* *********************************************** Author :111qqz Created Time :Mon 31 Oct 2016 02:47:01 PM CST File Name :code/hdu/4291.cpp ************************************************ */ #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> #include <vector> #include <queue> #include <set> #include <map> #include <string> #include <cmath> #include <cstdlib> #include <ctime> #define fst first #define sec second #define lson l,m,rt<<1 #define rson m+1,r,rt<<1|1 #define ms(a,x) memset(a,x,sizeof(a)) typedef long long LL; #define pi pair < int ,int > #define MP make_pair using namespace std; const double eps = 1E-8; const int dx4[4]={1,0,0,-1}; const int dy4[4]={0,-1,1,0}; const int inf = 0x3f3f3f3f; const LL loop0 = 1E9+7; const LL loop1=222222224,loop2=183120; //暴力得到每一层的循环节。 /* LL find_loop(LL mod) { LL a,b,c; a = 0 ; b = 1 ; for ( int i = 2 ; ; i++) { c = (a+3*b%mod)% mod; a = b; b = c; // cout<<"a:"<<a<<" b:"<<b<<endl; if (a==0&&b==1) { return i-1; // i的时候得到第i项，判断循环的时候是判断g[i-1]==g[0]&&g[i]==g[1]，所以循环节长度是i-1; } } }*/ LL n; struct Mat { LL mat[2][2]; void clear() { ms(mat,0); } void out() { cout<<mat[0][0]<<" "<<mat[0][1]<<endl; cout<<mat[1][0]<<" "<<mat[1][1]<<endl; cout<<endl; } }M,M1; Mat mul(Mat a,Mat b,LL MOD) { Mat c; c.clear(); for ( int i = 0 ; i < 2 ; i++) for ( int j = 0 ; j < 2 ; j++) for ( int k = 0 ; k < 2 ; k++) c.mat[i][j] = (c.mat[i][j] + a.mat[i][k]*b.mat[k][j]%MOD)%MOD; return c; } Mat power(Mat a,LL b,LL MOD) { Mat res; res.clear(); for ( int i = 0 ;i < 2 ; i++) res.mat[i][i] = 1; while (b>0) { if (b&1) res = mul(res,a,MOD); b = b>>1LL; a = mul(a,a,MOD); } return res; } void Mat_init() { M.clear(); M.mat[0][1] = 1; M.mat[1][0] = 1; M.mat[1][1] = 3; M1.clear(); M1.mat[1][0] = 1; } int main() { #ifndef ONLINE_JUDGE freopen("code/in.txt","r",stdin); #endif /* loop1 = find_loop(loop0); printf("loop1:%lld\n",loop1); loop2 = find_loop(loop1); printf("loop2:%lld\n",loop2); */ Mat_init(); while (scanf("%lld\n",&n)!=EOF) { if (n==0) { puts("0"); continue; } if (n==1) { puts("1"); continue; } LL cur = n; Mat ans; ans = power(M,cur-1,loop2); ans = mul(ans,M1,loop2); cur = ans.mat[1][0]; if (cur!=0&&cur!=1) //三次快速幂，每次都要记得判断初始项。 { ans = power(M,cur-1,loop1); ans = mul(ans,M1,loop1); cur = ans.mat[1][0]; } if (cur!=0&&cur!=1) { ans = power(M,cur-1,loop0); ans = mul(ans,M1,loop0); cur = ans.mat[1][0]; } printf("%lld\n",cur); } #ifndef ONLINE_JUDGE fclose(stdin); #endif return 0; } |
再来一个比较一般的做法：
参考递推式循环节
Acdreamer的博客_广义Fibonacci数列找循环节
摘重点：
今天早上起来后，看了下代码，为什么要判断5是不是p的模二次剩余呢，为什么是5呢
想了想，5对于斐波那契数列来讲，不就是x^2=x+1的delta么，那么这题的递推式是x^2=3*x+1，delta=3*3+4=13，然后我就把勒让德符号判断二次剩余那里改成13，然后对应的暴力出13及13以内的素数对应的循环节，交了一发，AC了
所以综上所述：
是模
的二次剩余时，枚举
的因子
的因子
对于第二种非剩余的情况，理论上是枚举(p+1)(p-1)的因子，实际上常常只是枚举2*(p+1)的因子
对于c=5的情况，是有定理：
不过对于c不等于5的情况。。。该结论是否一定成立呢…感觉不是很好证，求指教。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 |
/* *********************************************** Author :111qqz Created Time :Mon 31 Oct 2016 08:22:17 PM CST File Name :code/hdu/4291_2.cpp ************************************************ */ #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> #include <vector> #include <queue> #include <set> #include <map> #include <string> #include <cmath> #include <cstdlib> #include <ctime> #define fst first #define sec second #define lson l,m,rt<<1 #define rson m+1,r,rt<<1|1 #define ms(a,x) memset(a,x,sizeof(a)) typedef long long LL; #define pi pair < int ,int > #define MP make_pair using namespace std; const double eps = 1E-8; const int dx4[4]={1,0,0,-1}; const int dy4[4]={0,-1,1,0}; const int inf = 0x3f3f3f3f; struct Mat { LL mat[2][2]; void clear() { ms(mat,0); } }M,M1; Mat mul (Mat a,Mat b,LL mod) { Mat c; c.clear(); for ( int i = 0 ; i < 2 ; i++) for ( int j = 0 ; j < 2 ; j++) for ( int k = 0 ; k < 2 ; k++) c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]%mod)%mod; return c; } Mat mat_ksm(Mat a,LL b,LL mod) { Mat res; res.clear(); for ( int i = 0 ; i < 2 ; i++) res.mat[i][i] = 1; while (b>0) { if (b&1) res = mul(res,a,mod); b = b >> 1LL; a = mul(a,a,mod); } return res; } LL gcd(LL a,LL b) { return b?gcd(b,a%b):a; } const int N = 1E6+7; bool prime[N]; int p[N]; void isprime() //一个普通的筛 { int cnt = 0 ; ms(prime,true); for ( int i = 2 ; i < N ; i++) { if (prime[i]) { p[cnt++] = i ; for ( int j = i+i ; j < N ; j+=i) prime[j] = false; } } } LL ksm( LL a,LL b,LL mod) { LL res = 1; while (b>0) { if (b&1) res = (res * a) % mod; b = b >> 1LL; a = a * a % mod; } return res; } LL legendre(LL a,LL p) //勒让德符号,判断二次剩余 { if (ksm(a,(p-1)>>1,p)==1) return 1; return -1; } LL pri[N],num[N];//分解质因数的底数和指数。 int cnt; //质因子的个数 void solve(LL n,LL pri[],LL num[]) { cnt = 0 ; for ( int i = 0 ; p[i] * p[i] <= n ; i++) { if (n%p[i]==0) { int Num = 0 ; pri[cnt] = p[i]; while (n%p[i]==0) { Num++; n/=p[i]; } num[cnt] = Num; cnt++; } } if (n>1) { pri[cnt] = n; num[cnt] = 1; cnt++; } } LL fac[N]; int cnt2; //n的因子的个数 void get_fac(LL n)//得到n的所有因子 { cnt2 = 0 ; for (int i = 1 ; i*i <= n ; i++) { if (n%i==0) { if (i*i!=n) { fac[cnt2++] = i ; fac[cnt2++] = n/i; } else fac[cnt2++] = i; } } } int get_loop( LL p) //暴力得到不大于13的素数的循环节。 { LL a,b,c; a = 0 ; b = 1 ; for ( int i = 2; ; i++) { c = (a+3*b%p)%p; a = b; b = c; if (a==0&&b==1) return i-1; } } /* 2->3 3->2 5->12 7->16 11->8 13->52 */ const LL LOOP[10]={3,2,12,16,8,52}; LL ask_loop(int id) { return LOOP[id]; } LL find_loop(LL n) { solve(n,pri,num); LL ans = 1; for ( int i = 0 ; i < cnt ; i++) { LL rec = 1; if (pri[i]<=13) rec = ask_loop(i); else { if (legendre(13,pri[i])==1) get_fac(pri[i]-1); else get_fac((pri[i]+1)*(3-1)); //为什么可以假设循环节不大于2*(p+1)??? sort(fac,fac+cnt2); for ( int j = 0 ; j < cnt2 ; j++) //挨个验证因子 { Mat tmp = mat_ksm(M,fac[j],pri[i]); //下标从0开始，验证fac[j]为循环节，应该看fib[0]==fib[fac[j]]和fib[1]==fib[fac[j]+1]是否成立 tmp = mul(tmp,M1,pri[i]); if (tmp.mat[0][0]==0&&tmp.mat[1][0]==1) { rec = fac[j]; break; } } } for ( int j = 1 ; j < num[i] ; j++) rec *=pri[i]; ans = ans / gcd(ans,rec) * rec; } return ans; } void init() { M.clear(); M.mat[0][1] = M.mat[1][0] = 1; M.mat[1][1] = 3; M1.clear(); M1.mat[1][0] = 1; } LL n; LL loop0 = 1E9+7; LL loop1,loop2; int main() { #ifndef ONLINE_JUDGE freopen("code/in.txt","r",stdin); #endif /* printf("%d\n",get_loop(2)); printf("%d\n",get_loop(3)); printf("%d\n",get_loop(5)); printf("%d\n",get_loop(7)); printf("%d\n",get_loop(11)); printf("%d\n",get_loop(13)); */ init(); isprime(); while (~scanf("%lld\n",&n)) { if (n==0||n==1) { printf("%lld\n",n); continue; } LL loop1 = find_loop(loop0); LL loop2 = find_loop(loop1); // printf("loop1:%lld loop2:%lld\n",loop1,loop2); LL cur = n; Mat ans = mat_ksm(M,cur-1,loop2); ans = mul(ans,M1,loop2); cur = ans.mat[1][0]; if (cur!=0&&cur!=1) { Mat ans = mat_ksm(M,cur-1,loop1); ans = mul(ans,M1,loop1); cur = ans.mat[1][0]; } if (cur!=0&&cur!=1) { Mat ans = mat_ksm(M,cur-1,loop0); ans = mul(ans,M1,loop0); cur = ans.mat[1][0]; } printf("%lld\n",cur); } #ifndef ONLINE_JUDGE fclose(stdin); #endif return 0; } |
说点什么
您将是第一位评论人！