mcfx's blog - 矩阵快速幂

mcfx's blog - 矩阵快速幂 /category/jzk/ zh-CN Tue, 20 Dec 2016 09:03:00 +0800 Tue, 20 Dec 2016 09:03:00 +0800 BZOJ 2432: [Noi2011]兔农 /archives/156/ /archives/156/ Tue, 20 Dec 2016 09:03:00 +0800 admin

###Description
农夫栋栋近年收入不景气，正在他发愁如何能多赚点钱时，他听到隔壁的小朋友在讨论兔子繁殖的问题。
问题是这样的：第一个月初有一对刚出生的小兔子，经过两个月长大后，这对兔子从第三个月开始，每个月初生一对小兔子。新出生的小兔子生长两个月后又能每个月生出一对小兔子。问第n个月有多少只兔子？
聪明的你可能已经发现，第n个月的兔子数正好是第n个Fibonacci(斐波那契)数。栋栋不懂什么是Fibonacci数，但他也发现了规律：第i+2个月的兔子数等于第i个月的兔子数加上第i+1个月的兔子数。前几个月的兔子数依次为：
1 1 2 3 5 8 13 21 34 …
栋栋发现越到后面兔子数增长的越快，期待养兔子一定能赚大钱，于是栋栋在第一个月初买了一对小兔子开始饲养。
每天，栋栋都要给兔子们喂食，兔子们吃食时非常特别，总是每k对兔子围成一圈，最后剩下的不足k对的围成一圈，由于兔子特别害怕孤独，从第三个月开始，如果吃食时围成某一个圈的只有一对兔子，这对兔子就会很快死掉。
我们假设死去的总是刚出生的兔子，那么每个月的兔子数仍然是可以计算的。例如，当k=7时，前几个月的兔子数依次为：
1 1 2 3 5 7 12 19 31 49 80 …
给定n，你能帮助栋栋计算第n个月他有多少对兔子么？由于答案可能非常大，你只需要告诉栋栋第n个月的兔子对数除p的余数即可。

###Input
输入一行，包含三个正整数n, k, p。

###Output
输出一行，包含一个整数，表示栋栋第n个月的兔子对数除p的余数。

###Sample Input
6 7 100

###Sample Output
7

###Solution
按模k余数写出数列，则为
1,1,...,a,0,
a,a,...,b,0,
b,b,...,c,0,
...
这个数列最后有可能每一排首项循环，也有可能不再出现0
每次求a的逆元，则其在fib数列中首次出现位置就是a后面下一个0的位置，如果没有0特判

###Code
```c++
#include<bits/stdc++.h>

typedef long long ll;

#define fo0(i,n) for(int i=0,i##end=n;i<i##end;i++)
#define fo1(i,n) for(int i=1,i##end=n;i<=i##end;i++)
#define fo(i,a,b) for(int i=a,i##end=b;i<=i##end;i++)

void exgcd(int a,int b,int &x,int &y){if(b)exgcd(b,a%b,y,x),y-=(a/b)*x;else x=1,y=0;}

ll n;
int k,p,pos[1000000],f[6000010]={0,1,1},p2[1000000];

struct mat
{
	int a[3][3];
	inline void operator *=(mat x)
	{
		ll t[3][3];memset(t,0,sizeof(t));
		fo0(i,3)fo0(j,3)fo0(k,3)t[i][k]+=(ll)a[i][j]*x.a[j][k];
		fo0(i,3)fo0(j,3)a[i][j]=t[i][j]%p;
	}
	inline mat operator *(mat x)
	{
		mat r=*this;r*=x;return r;
	}
};

mat pow(mat a,ll p)
{
	mat r=(mat){1,0,0,0,1,0,0,0,1};
	for(;p;p>>=1,a*=a)if(p&1)r*=a;
	return r;
}

struct data
{
	int st,len;mat x;
}s[1000001];

inline mat get(ll n,int tt)
{
	mat r=(mat){1,0,0,0,1,0,0,0,1};
	if(!n)return r;
	for(;n>=s[tt].len;n-=s[tt++].len)r*=s[tt].x;
	return r*pow((mat){0,1,0,1,1,0,0,0,1},n);
}

int main()
{
	scanf("%lld%d%d",&n,&k,&p);
	for(int i=3;;i++)
	{
		f[i]=f[i-1]+f[i-2];
		if(f[i]>=k)f[i]-=k;
		if(!pos[f[i]])pos[f[i]]=i;
		if(f[i]==1&&f[i-1]==0)break;
	}
	int ns=1,sc=1;
	for(;!p2[ns];sc++)
	{
		s[sc].st=ns;
		p2[ns]=sc;
		int t1,t2;
		exgcd(k,ns,t1,t2);
		if(t2<0)t2+=k;
		t1=pos[t2];
		if((ll)ns*t2%k!=1||!t1)break;
		s[sc].len=t1;
		s[sc].x=pow((mat){0,1,0,1,1,0,0,0,1},t1-1)*(mat){0,1,0,1,1,0,0,-1,1};
		ns=(ll)ns*f[t1-1]%k;
	}
	mat bfm=(mat){1,0,0,0,1,0,0,0,1},lpm=(mat){1,0,0,0,1,0,0,0,1},ans=(mat){1,0,0,0,1,0,0,0,1};
	ll bf=0,lp=0;
	fo1(i,p2[ns]-1)bf+=s[i].len,bfm*=s[i].x;
	fo(i,p2[ns],sc-1)lp+=s[i].len,lpm*=s[i].x;
	if(!lp)lp=1,lpm=(mat){0,1,0,1,1,0,0,0,1};
	if(n<bf)
	{
		ans*=get(n,1);
	}
	else
	{
		ans*=bfm;n-=bf;
		ans*=pow(lpm,n/lp);
		ans*=get(n%lp,p2[ns]);
	}
	int a1=ans.a[0][1]+ans.a[2][1];
	while(a1>=p)a1-=p;while(a1<0)a1+=p;
	printf("%d\n",a1);
}
```

]]> 0 /archives/156/#comments /feed/archives/156/ BZOJ 1009: [HNOI2008]GT考试 /archives/21/ /archives/21/ Fri, 09 Sep 2016 20:13:00 +0800 admin

###Description
阿申准备报名参加GT考试，准考证号为N位数X1X2....Xn(0<=Xi<=9),他不希望准考证号上出现不吉利的数字。他的不吉利数学A1A2...Am(0<=Ai<=9)有M位，不出现是指X1X2...Xn中没有恰好一段等于A1A2...Am. A1和X1可以为0

###Input
第一行输入N,M,K.接下来一行输入M位的数。 N<=10^9,M<=20,K<=1000

###Output
阿申想知道不出现不吉利数字的号码有多少种，输出模K取余的结果.

###Sample Input
4 3 100
111

###Sample Output
81

###Solution
我们把不吉利的数字用s表示，用f[i][j]表示n[i]匹配到s[j]的情况数，则可以用一个转移矩阵a通过f[i-1][j]求出f[i][j]
a[i][j]表示在s中匹配了i位时，在后面添数，有多少种情况匹配到j位，也就是指s的前i位后面添数后最长的公共前后缀长度为j的情况数
kmp预处理，然后对于每个i，枚举添加的数，再仿照kmp进行处理，即可求出转移矩阵
矩阵快速幂优化就可以过了，时间复杂度$$O(logN*M^3)$$

###Code
```c++
#include<cstdio>
 
struct _matrix
{
    int s[21][21],a,b;
    _matrix()
    {
        for(int i=0;i<21;i++)
            for(int j=0;j<21;j++)
                s[i][j]=0;
    }
}temp;
 
int mod;
 
void mul(_matrix &x,_matrix y)
{
    for(int i=0;i<x.a;i++)
        for(int j=0;j<y.b;j++)
        {
            temp.s[i][j]=0;
            for(int k=0;k<x.b;k++)
                temp.s[i][j]+=x.s[i][k]*y.s[k][j];
            temp.s[i][j]%=mod;
        }
    temp.a=x.a;
    temp.b=y.b;
    x=temp;
}
 
int m[21],n,f[21];
 
int main()
{
    int mm;
    scanf("%d%d%d",&n,&mm,&mod);
    char t[20];
    scanf("%s",t);
    for(int i=0;i<mm;i++)m[i+1]=t[i]-48;
    f[1]=0;
    for(int i=2;i<=mm;i++)
    {
        int t=i-1;
        while(t!=0&&m[f[t]+1]!=m[i])t=f[t];
        if(t==0)
        {
            if(m[1]==m[i])
                f[i]=1;
            else
                f[i]=0;
        }
        else
        {
            f[i]=f[t]+1;
        }
    }
    _matrix a;
    a.a=mm+1;
    a.b=mm+1;
    a.s[0][0]=9,a.s[0][1]=1;
    for(int i=1;i<mm;i++)
    {
        for(int j=0;j<10;j++)
        {
            if(m[i+1]==j)
            {
                a.s[i][i+1]++;
                continue;
            }
            int t=i;
            while(t!=0&&m[f[t]+1]!=j)t=f[t];
            if(t==0)
            {
                if(m[1]==j)
                    a.s[i][t+1]++;
                else
                    a.s[i][0]++;
            }
            else
                a.s[i][f[t]+1]++;
        }
    }
    _matrix b;
    b.a=1;
    b.b=mm+1;
    b.s[0][0]=9;
    b.s[0][1]=1;
    int ans=0,to=n-1,no=1;
    while(no<to)
    {
        if(to&no)mul(b,a);
        mul(a,a);
        no<<=1;
    }
    for(int i=0;i<mm;i++)ans+=b.s[0][i];
    printf("%d\n",ans%mod);
    return 0;
}
```

]]> 0 /archives/21/#comments /feed/archives/21/