2019-ACM-ICPC-南昌区网络赛-H. The Nth Item-特征根法求通项公式+二次剩余+欧拉降幂

2019-ACM-ICPC-南昌区网络赛-H. The Nth Item-特征根法求通项公式+二次剩余+欧拉降幂

---

【Problem Description】

​ 已知\(f(n)=3\cdot f(n-1)+2\cdot f(n-2),(n\ge 2)\)，求\(f(n)\pmod {998244353}\)。

【Solution】

​ 利用特征根法求得通项公式为\(a_n=\frac{\sqrt{17}}{17}\cdot\Bigg(\Big(\frac{3+\sqrt{17}}{2} \Big)^n-\Big(\frac{3-\sqrt{17}}{2} \Big)^n\Bigg)\)。\(\sqrt{17}\pmod {998244353}\)可以用二次剩余求得。然后就可以用快速幂在\(O(log(n))\)的时间复杂度内求得\(a_n\)。但是因为\(T\le 10^7\)，所以还需优化。

​ \(n\le 10^{18}\)，进行欧拉降幂后\(n\le 10^9\)，令\(k=\lfloor \sqrt{n}\rfloor\)，则：

\[n=k\cdot t+r\\ x^n=x^{k\cdot t+r}\Leftrightarrow x^n=x^{k\cdot t}+x^r(t,r\le k) \]

然后预处理出\(x^r,r\le k\)以及\(x^{k\cdot t},t\le k\)。则\(x^n\)就可以\(O(1)\)查询了。

---

【Code】

/*
 * @Author: Simon 
 * @Date: 2019-09-08 13:13:29 
 * @Last Modified by: Simon
 * @Last Modified time: 2019-09-17 17:44:39
 */
#include<bits/stdc++.h>
using namespace std;
typedef int Int;
#define int long long
#define INF 0x3f3f3f3f
#define maxn 100005
const int mod=998244353;
int inv17;
/*
 * Author: Simon
 * 功能: 求解x^2=n(mod p),即x=sqrt(n)(mod p)
 * 复杂度: O(sqrt(p))
 */

/*类似复数 单位元为w(复数的单位元为-1)*/
struct Complex {
    int x, y, w;
    Complex() {}
    Complex(int x, int y, int w) : x(x), y(y), w(w) {}
};
/*类复数乘法 */
Complex mul(Complex a, Complex b, int p) {
    Complex ans;
    ans.x = (a.x * b.x % p + a.y * b.y % p * a.w % p) % p;
    ans.y = (a.x * b.y % p + a.y * b.x % p) % p;
    ans.w = a.w;
    return ans;
}
/*类复数快速幂 */
Complex Complexfpow(Complex a, int b, int mod) {
    Complex ans = Complex(1, 0, a.w);
    while (b) {
        if (b & 1) ans = mul(ans, a, mod);
        a = mul(a, a, mod);
        b >>= 1;
    }
    return ans;
}
int fpow(int a, int b, int mod) {
    int ans = 1;
    a %= mod;
    while (b) {
        if (b & 1) (ans *= a) %= mod;
        (a *= a) %= mod;
        b >>= 1;
    }
    return ans;
}
/*求解x^2=n(mod p) */
int solve(int n, int p) {
    n %= p;
    if (n == 0) return 0;
    if (p == 2) return n;
    if (fpow(n, (p - 1) / 2, p) == p - 1) return -1; /*勒让德定理判断n不是p的二次剩余 */
    mt19937 rnd(time(0));
    int a, t, w;
    do {
        a = rnd() % p;
        t = a * a - n;
        w = (t % p + p) % p;                    /*构造w=a^2-n */
    } while (fpow(w, (p - 1) / 2, p) != p - 1); /*找到一个w不是p的二次剩余 */
    Complex ans = Complex(a, 1, w);
    ans = Complexfpow(ans, (p + 1) / 2, p); /*答案为(a+w)^{(p+1)/2} */
    return ans.x;
}
pair<int,int> bit1[maxn],bit2[maxn];
Int main(){
#ifndef ONLINE_JUDGE
    //freopen("input.in","r",stdin);
    //freopen("output.out","w",stdout);
#endif
    ios::sync_with_stdio(false);
    cin.tie(0);inv17=fpow(17,mod-2,mod);
    int x=solve(17,mod); //二次剩余
    int lim=ceil(sqrt(1e9));
    int s1=(3+x)%mod*fpow(2,mod-2,mod)%mod,
        s2=(3-x)%mod*fpow(2,mod-2,mod)%mod;
    bit1[0].first=bit1[0].second=bit2[0].first=bit2[0].second=1;
    for(int i=1;i<=lim;i++){ //预处理
        bit1[i].first=bit1[i-1].first*s1%mod;
        bit1[i].second=bit1[i-1].second*s2%mod;
        bit2[i].first=fpow(s1,i*lim,mod);
        bit2[i].second=fpow(s2,i*lim,mod);
    }
    int q,n;cin>>q>>n;
    int ans=0;
    while(q--){
        int tmp=0;
        if(n==0) tmp=0;
        else if(n==1) tmp=1;
        else{
            int t=n%(mod-1);
            int t2=t/lim,t1=t%lim;
            tmp=(bit1[t1].first*bit2[t2].first%mod-bit1[t1].second*bit2[t2].second%mod)%mod*x%mod*inv17%mod;
        }
        tmp=(tmp+mod)%mod;
        n=tmp*tmp^n; ans^=tmp;
        // cout<<n<<endl;
    }
    cout<<(ans+mod)%mod<<endl;
#ifndef ONLINE_JUDGE
    cout<<endl;system("pause");
#endif
    return 0;
}