线性代数二之矩阵加速DP——数学作业,Arc of Dream

数学作业

description

solution

d p dp dp状态转移方程, d p i = d p i − 1 ∗ 1 0 l e n i + i ( m o d M ) dp_{i}=dp_{i-1}*10^{len_i}+i\pmod M dpi=dpi110leni+i(modM)

n n n巨大,分段矩阵加速
[ f i − 1 i − 1 1 ] × [ 1 0 l e n 0 0 1 1 0 1 1 1 ] = [ f i i 1 ] \begin{bmatrix} f_{i-1}&i-1&1 \end{bmatrix} \times \begin{bmatrix} 10^{len}&0&0\\ 1&1&0\\ 1&1&1 \end{bmatrix}= \begin{bmatrix} f_{i}&i&1 \end{bmatrix} [fi1i11]×10len11011001=[fii1]

code

#include <cmath>
#include <cstdio>
#include <cstring>
#define int long long
int n, mod;
int mi[20];
struct matrix {
	int n, m;
	int c[3][3];
	matrix() {
		n = m = 0;
		memset( c, 0, sizeof( c ) );
	}
	matrix operator * ( matrix &t ) {
		matrix ans;
		ans.n = n, ans.m = t.m;
		for( int i = 0;i <= n;i ++ )
			for( int j = 0;j <= t.m;j ++ )
				for( int k = 0;k <= m;k ++ )
					ans.c[i][j] = ( ans.c[i][j] + c[i][k] * t.c[k][j] % mod ) % mod;
		return ans;
	}
}fast, ret;

matrix qkpow( matrix x, int y ) {
	if( ! y ) return x;
	matrix ans;
	ans.n = ans.m = 2;
	for( int i = 0;i <= 2;i ++ )
		ans.c[i][i] = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x;
		x = x * x;
		y >>= 1;
	}
	return ans;
}

signed main() {
	scanf( "%lld %lld", &n, &mod );
	mi[0] = 1;
	for( int i = 1;i <= 18;i ++ )
		mi[i] = mi[i - 1] * 10;
	ret.n = 0, ret.m = 2, ret.c[0][2] = 1;
	fast.n = fast.m = 2;
	for( int i = 1;;i ++ ) {
		memset( fast.c, 0, sizeof( fast.c ) );
		fast.c[0][0] = mi[i] % mod;
		fast.c[1][0] = fast.c[1][1] = fast.c[2][0] = fast.c[2][1] = fast.c[2][2] = 1;
		if( pow( 10, i ) <= n )
			fast = qkpow( fast, mi[i - 1] * 9 );
		else
			fast = qkpow( fast, n - mi[i - 1] + 1 );
		ret = ret * fast;
		if( mi[i] > n ) break;
	}
	printf( "%lld\n", ret.c[0][0] % mod );
	return 0;
}

Arc of Dream

description

solution

A o D ( n ) = A o D ( n − 1 ) + a n − 1 b n − 1 = A o D ( n − 1 ) + ( a n − 2 ∗ A X + A Y ) ( b n − 2 ∗ B X + B Y ) AoD(n)=AoD(n-1)+a_{n-1}b_{n-1}=AoD(n-1)+(a_{n-2}*AX+AY)(b_{n-2}*BX+BY) AoD(n)=AoD(n1)+an1bn1=AoD(n1)+(an2AX+AY)(bn2BX+BY)

= A o D ( n − 1 ) + a n − 2 b n − 2 A X B X + a n − 2 A X B Y + b n − 2 B X A Y + A Y B Y =AoD(n-1)+a_{n-2}b_{n-2}AXBX+a_{n-2}AXBY+b_{n-2}BXAY+AYBY =AoD(n1)+an2bn2AXBX+an2AXBY+bn2BXAY+AYBY

满满的矩阵味道,赤裸裸的系数暗示
在这里插入图片描述

[ A o D ( i ) a i − 1 b i − 1 a i − 1 b i − 1 A Y B Y A Y B Y ] × [ 1 0 0 0 0 0 0 A X B X A X B X 0 0 0 0 0 A X A X A X 0 0 0 0 B X B X 0 B X 0 0 0 1 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 1 ] ∣ ∣ [ A o D ( i + 1 ) a i b i a i b i A Y B Y A Y B Y ] \begin{bmatrix} AoD(i)&a_{i-1}b_{i-1}&a_{i-1}&b_{i-1}&AYBY&AY&BY \end{bmatrix} \\ \times \\ \begin{bmatrix} 1&0&0&0&0&0&0\\ AXBX&AXBX&0&0&0&0&0\\ AX&AX&AX&0&0&0&0\\ BX&BX&0&BX&0&0&0\\ 1&1&0&0&1&0&0\\ 1&1&1&0&0&1&0\\ 1&1&0&1&0&0&1\\ \end{bmatrix} \\ || \\ \begin{bmatrix} AoD(i+1)&a_ib_i&a_i&b_i&AYBY&AY&BY \end{bmatrix} [AoD(i)ai1bi1ai1bi1AYBYAYBY]×1AXBXAXBX1110AXBXAXBX11100AX0010000BX001000010000000100000001[AoD(i+1)aibiaibiAYBYAYBY]

code

#include <cstdio>
#include <cstring>
#define int long long
#define mod 1000000007
struct matrix {
	int n, m;
	int c[8][8];
	matrix() {
		n = m = 0;
		memset( c, 0, sizeof( c ) );
	}
	matrix operator * ( matrix &t ) {
		matrix ans;
		ans.n = n, ans.m = t.m;
		for( int i = 1;i <= n;i ++ )
			for( int j = 1;j <= t.m;j ++ )
				for( int k = 1;k <= m;k ++ )
					ans.c[i][j] = ( ans.c[i][j] + c[i][k] * t.c[k][j] ) % mod;
		return ans;
	}
}g, ret;

matrix qkpow( matrix x, int y ) {
	matrix ans;
	ans.n = x.n, ans.m = x.m;
	for( int i = 1;i <= ans.n;i ++ )
		ans.c[i][i] = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x;
		x = x * x;
		y >>= 1;
	}
	return ans;
}

signed main() {
	int n, A0, AX, AY, B0, BX, BY;
	while( ~ scanf( "%lld %lld %lld %lld %lld %lld %lld", &n, &A0, &AX, &AY, &B0, &BX, &BY ) ) {
		g.n = g.m = ret.m = 7, ret.n = 1;
		if( ! n ) {
			printf( "0\n" );
			continue;
		}
		if( n == 1 ) {
			printf( "%lld\n", A0 * B0 % mod );
			continue;
		}
		memset( g.c, 0, sizeof( g.c ) );
		memset( ret.c, 0, sizeof( ret.c ) );
		ret.c[1][1] = ret.c[1][2] = A0 * B0 % mod;
		ret.c[1][3] = A0, ret.c[1][4] = B0;
		ret.c[1][5] = AY * BY % mod;
		ret.c[1][6] = AY, ret.c[1][7] = BY;
		g.c[1][1] = g.c[5][1] = g.c[5][2] = g.c[6][3] = g.c[7][4] = g.c[5][5] = g.c[6][6] = g.c[7][7] = 1;
		g.c[2][1] = g.c[2][2] = AX * BX % mod;
		g.c[3][1] = g.c[3][2] = AX * BY % mod;
		g.c[4][1] = g.c[4][2] = AY * BX % mod;
		g.c[3][3] = AX, g.c[4][4] = BX;
		g = qkpow( g, n - 1 );
		ret = ret * g;
		printf( "%lld\n", ret.c[1][1] );
	}
	return 0;
}

HH去散步

之前已经写过博客了,不再赘述

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值