HDU 6395 Sequence(分块矩阵快速幂)

最新推荐文章于 2022-09-06 10:35:25 发布

原创最新推荐文章于 2022-09-06 10:35:25 发布 · 313 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#hdu多校 #矩阵快速幂 #分块

矩阵快速幂专栏收录该内容

13 篇文章

订阅专栏

本文介绍了一种特定序列的计算方法，利用快速幂和分段思想解决复杂数学问题，并提供了两种实现方式：直接分段和使用二分查找确定区间。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Sequence

Time Limit: 4000/2000 MS (Java/Others) Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 1585 Accepted Submission(s): 568

Problem Description
Let us define a sequence as below
$\left\{\begin{matrix} F_1=A\\ F_2=B\\ F_n=C*F_{n-2}+D*F_{n-1}+\left \lfloor \frac{p}{n} \right \rfloor \end{matrix}\right.$
Your job is simple, for each task, you should output Fn module 109+7.

Input
The first line has only one integer T, indicates the number of tasks.

Then, for the next T lines, each line consists of 6 integers, A , B, C, D, P, n.

$1 \leq T \leq 20$
$0≤A,B,C,D≤10^9$
$1≤P,n≤10^9$

Sample Input
2
3 3 2 1 3 5
3 2 2 2 1 4

Sample Output
36
24

思路

如果没有 $⌊Pn⌋\left \lfloor \frac{P}{n} \right \rfloor$ 就是一个简单的快速幂，但是现在有了这个，那么我们会发现 $⌊Pi⌋\left \lfloor \frac{P}{i} \right \rfloor$ ，在 $P$ 固定的时候 $⌊Pi⌋\left \lfloor \frac{P}{i} \right \rfloor$ 的值是分段的，例如我们让P为30，下面是不同的i所对应的 $⌊Pi⌋\left \lfloor \frac{P}{i} \right \rfloor$ 的值
$[1][2][3][4][5][6][7][8,10][11,15][16,30]3015107654321\begin{matrix} [1]&[2]&[3]&[4]&[5]&[6]&[7]&[8,10]&[11,15]&[16,30]\\ 30 &15 & 10 & 7 & 6 &5 & 4 & 3 & 2 &1 \end{matrix}$
因为在每个区间内 $⌊Pi⌋\left \lfloor \frac{P}{i} \right \rfloor$ 的值是固定的所以我们可以构造出快速幂的公式
$(FnFn−1⌊Pn⌋)=(DC1100001)(Fn−1Fn−2⌊Pn⌋)=(DC1100001)n−i(FiFi−1⌊Pn⌋)\begin{pmatrix} F_n\\ F_{n-1}\\ \left \lfloor \frac{P}{n} \right \rfloor \end{pmatrix}=\begin{pmatrix} D & C &1 \\ 1&0 &0 \\ 0&0 &1 \end{pmatrix}\begin{pmatrix} F_{n-1}\\ F_{n-2}\\ \left \lfloor \frac{P}{n} \right \rfloor \end{pmatrix}=\begin{pmatrix} D & C &1 \\ 1&0 &0 \\ 0&0 &1 \end{pmatrix}^{n-i}\begin{pmatrix} F_{i}\\ F_{i-1}\\ \left \lfloor \frac{P}{n} \right \rfloor \end{pmatrix}$
那现在是如何找出每段 $⌊Pi⌋\left \lfloor \frac{P}{i} \right \rfloor$ 的开始和结尾，通过观察式子我们会发现每次的每个区间的末项就是 $⌊PPi⌋\left \lfloor \frac{P}{ \frac{P}{i}} \right \rfloor$ ， $⌊PPi⌋+1\left \lfloor \frac{P}{ \frac{P}{i}} \right \rfloor+1$ 就是下一项的开始，这个分段的思想在很多莫比乌斯反演的题目中常用来优化，既然有了每一段的首项和末项，我们每次记录一下每一段区间开始前的前两项用来求得这一段最后一个数，然后再来判断一下， $⌊PPi⌋\left \lfloor \frac{P}{ \frac{P}{i}} \right \rfloor$ 和 $m i n (n, p)$ 哪个小即可，当i大于p的时候 $⌊Pi⌋\left \lfloor \frac{P}{i} \right \rfloor$ 为0，那就是一个二维的矩阵快速幂了， $3 x 3$ 矩阵中 $2 x 2$ 的部分直接拿来用就可以了

#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#define MAX 3
using namespace std;
typedef struct {
long long m[MAX][MAX];
}Matrix;
Matrix P;
Matrix I={1,0,0,0,1,0,0,0,1};
const long long mod=1e9+7;
Matrix Matrixmul(Matrix a,Matrix b)
{
    int i,j,k;
    Matrix c;
    for(i=0;i<MAX;i++)
        for(j=0;j<MAX;j++)
        {
            c.m[i][j]=0;
            for(k=0;k<MAX;k++)
            {
                c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;
            }
            c.m[i][j]%=mod;
        }
        return c;
}
Matrix quickpow(long long n)
{
    Matrix m=P,b=I;
    while(n>0)
    {
        if(n%2==1)
            b=Matrixmul(b,m);
        n=n/2;
        m=Matrixmul(m,m);
    }
    return b;
}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        memset(P.m,0,sizeof(P.m));
        int a,b,c,d,n,p;
        scanf("%d%d%d%d%d%d",&a,&b,&c,&d,&p,&n);
        if(n==2)
        {
            printf("%d\n",a);
            continue;
        }
        else if(n==1)
        {
            printf("%d\n",b);
            continue;
        }
        P.m[0][0]=d%mod;
        P.m[0][1]=c%mod;
        P.m[0][2]=1;
        P.m[1][0]=1;
        P.m[2][2]=1;
        long long ans=0;
        long long ta=a%mod;
        long long tb=b%mod;
        int m=min(n,p);
        for(int i=3,last;i<=m;i=last+1)
        {
            last=min(p/(p/i),m);
            Matrix A=quickpow(last-i+1);
            ans=(A.m[0][0]*tb%mod+A.m[0][1]*ta%mod+A.m[0][2]*(p/i)%mod+mod)%mod;
            ta=(A.m[1][0]*tb%mod+A.m[1][1]*ta%mod+A.m[1][2]*(p/i)%mod+mod)%mod;
            tb=ans;
        }
        if(n>p)
        {
            Matrix A=quickpow(n-max(p,2));
            ans=(A.m[0][0]*tb%mod+A.m[0][1]*ta%mod+mod)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}

我们还可以用二分来寻找区间的首项和尾项

#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#define MAX 3
using namespace std;
typedef struct
{
    long long m[MAX][MAX];
} Matrix;
Matrix P;
Matrix I= {1,0,0,0,1,0,0,0,1};
const long long mod=1e9+7;
Matrix Matrixmul(Matrix a,Matrix b)
{
    int i,j,k;
    Matrix c;
    for(i=0; i<MAX; i++)
        for(j=0; j<MAX; j++)
        {
            c.m[i][j]=0;
            for(k=0; k<MAX; k++)
            {
                c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;
            }
            c.m[i][j]%=mod;
        }
    return c;
}
Matrix quickpow(long long n)
{
    Matrix m=P,b=I;
    while(n>0)
    {
        if(n%2==1)
            b=Matrixmul(b,m);
        n=n/2;
        m=Matrixmul(m,m);
    }
    return b;
}
int erfen(int n,int x)
{
    int l=1,r=n;
    int mid;
    while(l<=r)
    {
        mid=(l+r)/2;
        if(n/mid>=x)
            l=mid+1;
        else
            r=mid-1;
    }
    return r;
}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        memset(P.m,0,sizeof(P.m));
        int a,b,c,d,n,p;
        scanf("%d%d%d%d%d%d",&a,&b,&c,&d,&p,&n);
        if(n==2)
        {
            printf("%d\n",a);
            continue;
        }
        else if(n==1)
        {
            printf("%d\n",b);
            continue;
        }
        P.m[0][0]=d%mod;
        P.m[0][1]=c%mod;
        P.m[0][2]=1;
        P.m[1][0]=1;
        P.m[2][2]=1;
        long long ans=0;
        long long ta=a%mod;
        long long tb=b%mod;
        int m=min(n,p);
        for(int i=3,last; i<=m; i=last+1)
        {
            last=min(erfen(p,p/i),m);
            Matrix A=quickpow(last-i+1);
            ans=(A.m[0][0]*tb%mod+A.m[0][1]*ta%mod+A.m[0][2]*(p/i)%mod+mod)%mod;
            ta=(A.m[1][0]*tb%mod+A.m[1][1]*ta%mod+A.m[1][2]*(p/i)%mod+mod)%mod;
            tb=ans;
        }
        if(n>p)
        {
            Matrix A=quickpow(n-max(p,2));
            ans=(A.m[0][0]*tb%mod+A.m[0][1]*ta%mod+mod)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}