斐波那契

题目大意

F[0]=0，F[1]=1，F[n]=F[n-1]+F[n-2]。
现在要求gcd(aF[n]+bF[n+1],cF[n]+dF[n+1])。
a,b,c,d<=1000，n<=1e18。

做法

我们可以先辗转相除。
即(a,b,c,d)=(c,d,a%c,b-(a/c)d)。
辗转到最后会把c变成0。
问题变成求gcd(aF[n]+bF[n+1],cF[n+1])。
先特判a=0或c=0，非常简单。
我们先令g=gcd(aF[n]+bF[n+1],F[n+1])。
这个很好求，显然g=gcd(aF[n],F[n+1])=gcd(a,F[n+1])=gcd(a,F[n+1]%a)。
因为大家都知道gcd(F[n],F[n+1])=1，这个可以归纳证明。
矩阵乘法求F[n+1]%a即可。
接下来考虑gcd(aF[n]+bF[n+1],cF[n+1])=g*gcd((aF[n]+bF[n+1])/g,cF[n+1]/g)。
然后我们发现了显然gcd((aF[n]+bF[n+1])/g,F[n+1]/g)=1。
因此答案是g*gcd((aF[n]+bF[n+1])/g,c)=gcd(aF[n]+bF[n+1],gc)。
同样在%gc意义下做矩阵乘法求斐波那契数即可。

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int mod=998244353;
int v[2][2],dis[2][2],ans[2][2],o[2][2],sta[80];
int i,j,k,l,t,a,b,c,d,g,x,y,top,ca,num,mo;
ll n,m;
int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}
void work(ll m,int mo){
            top=0;
            while (m){
                sta[++top]=m%2;
                m/=2;
            }
            ans[0][0]=ans[1][1]=1;
            ans[0][1]=ans[1][0]=0;
            while (top){
                fo(i,0,1)
                    fo(j,0,1)
                        o[i][j]=0;
                fo(k,0,1){
                    /*fo(i,0,1)
                        fo(j,0,1)
                            (o[i][j]+=(ll)ans[i][k]*ans[k][j]%mo)%=mo;*/
                    (o[0][0]+=(ll)ans[0][k]*ans[k][0]%mo)%=mo;
                    (o[0][1]+=(ll)ans[0][k]*ans[k][1]%mo)%=mo;
                    (o[1][0]+=(ll)ans[1][k]*ans[k][0]%mo)%=mo;
                    (o[1][1]+=(ll)ans[1][k]*ans[k][1]%mo)%=mo;
                }
                fo(i,0,1)
                    fo(j,0,1)
                        ans[i][j]=o[i][j];
                if (sta[top]){
                    fo(i,0,1)
                        fo(j,0,1)
                            o[i][j]=0;
                    fo(k,0,1){
                        /*fo(i,0,1)
                            fo(j,0,1)
                                (o[i][j]+=(ll)ans[i][k]*dis[k][j]%mo)%=mo;*/
                        (o[0][0]+=(ll)ans[0][k]*dis[k][0]%mo)%=mo;
                        (o[0][1]+=(ll)ans[0][k]*dis[k][1]%mo)%=mo;
                        (o[1][0]+=(ll)ans[1][k]*dis[k][0]%mo)%=mo;
                        (o[1][1]+=(ll)ans[1][k]*dis[k][1]%mo)%=mo;
                    }
                    fo(i,0,1)
                        fo(j,0,1)
                            ans[i][j]=o[i][j];
                }
                top--;
            }
}
int main(){
    freopen("fib.in","r",stdin);freopen("fib.out","w",stdout);
    dis[0][1]=dis[1][0]=dis[1][1]=1;
    scanf("%d",&ca);
    while (ca--){
        scanf("%lld%d%d%d%d",&n,&a,&b,&c,&d);
        while (c){
            x=a%c;
            y=b-(a/c)*d;
            a=c;b=d;
            c=x;d=y;
        }
        c=d;
        if (c<0) c=-c;
        if (a<0){
            a=-a;
            b=-b;
        }
        if (c==0){
            work(n,mod);
            num=(ll)ans[1][0]*a%mod;
            work(n+1,mod);
            (num+=(ll)ans[1][0]*b%mod)%=mod;
            (num+=mod)%=mod;
            printf("%d\n",num);
            continue;
        }
        if (a==0){
            num=gcd(b,c);
            work(n+1,mod);
            num=(ll)num*ans[1][0]%mod;
            (num+=mod)%=mod;
            printf("%d\n",num);
            continue;
        }
        work(n+1,a);
        g=gcd(a,ans[1][0]);
        t=0;
        work(n,g*c);
        t=(ll)a*ans[1][0]%(g*c);
        work(n+1,g*c);
        (t+=(ll)b*ans[1][0]%(g*c))%=(g*c);
        (t+=(g*c))%=(g*c);
        num=gcd(t,g*c);
        num%=mod;
        printf("%d\n",num);
    }
}