数论练习二之BSGS算法——随机数生成器,Matrix,Lunar New Year and a Recursive Sequence,Fermat‘s Last Theorem

2 篇文章 0 订阅

[SDOI2013] 随机数生成器

description

solution

肯定是非常想找一个通项公式来表示第 n n n个数的
在这里插入图片描述

依据形式,考虑化成等比数列

x i + 1 + k = a ( x i + k ) = a ⋅ x i + b + t ⇒ k = b a − 1 x_{i+1}+k=a(x_i+k)=a·x_i+b+t\Rightarrow k=\frac{b}{a-1} xi+1+k=a(xi+k)=axi+b+tk=a1b

⇒ x i + b a − 1 = a i − 1 ( x 1 + b a − 1 ) ( m o d p ) ⇒ a i − 1 ≡ ( x i + b ⋅ ( a − 1 ) − 1 ) ⋅ ( x 1 + b ⋅ ( a − 1 ) − 1 ) − 1 ( m o d p ) \Rightarrow x_i+\frac{b}{a-1}=a^{i-1}(x_1+\frac{b}{a-1})\pmod p\Rightarrow a^{i-1}\equiv (x_i+b·(a-1)^{-1})·(x_1+b·(a-1)^{-1})^{-1}\pmod p xi+a1b=ai1(x1+a1b)(modp)ai1(xi+b(a1)1)(x1+b(a1)1)1(modp)

x i x_i xi就是题目中的 t t t,所以式子只有 a i − 1 a^{i-1} ai1这个指数未知,bsgs求解最小指数( p p p是质数)

但要先特判 a = 0 / 1 a=0/1 a=0/1两种情况,显然此情况下等比数列式子是不成立的

  • a=0

    x i = b x_i=b xi=b

  • a=1

    x i ≡ x 1 + b ( i − 1 ) ( m o d p ) ⇔ ( t − x 1 ) b − 1 ≡ i − 1 ( m o d p ) x_i\equiv x_1+b(i-1)\pmod p\Leftrightarrow (t-x_1)b^{-1}\equiv i-1\pmod p xix1+b(i1)(modp)(tx1)b1i1(modp)

code

#include <map>
#include <cmath>
#include <cstdio>
using namespace std;
#define int long long
map < int, int > mp;

int bsgs( int a, int r, int p ) {
	if( a % p == 0 ) return -1;
	mp.clear();
	int m = ceil( sqrt( p * 1.0 ) );
	int IS = 1;
	for( int i = 1;i <= m;i ++ ) {
		IS = IS * a % p;
		mp[IS * r % p] = i;
	}
	int SI = 1;
	for( int i = 1;i <= m;i ++ ) {
		SI = SI * IS % p;
		if( mp.count( SI ) ) 
			return i * m - mp[SI] + 1;
	}
	return -1;
}

int T, p, a, b, x1, t;

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % p;
		x = x * x % p;
		y >>= 1;
	}
	return ans;
}

int inv( int x ) {
	return qkpow( x, p - 2 );
}

signed main() {
	scanf( "%lld", &T );
	while( T -- ) {
		scanf( "%lld %lld %lld %lld %lld", &p, &a, &b, &x1, &t );
		if( x1 == t ) {
			printf( "1\n" );
			continue;
		}
		if( a == 0 ) {
			if( x1 == t ) printf( "1\n" );
			else if( b == t ) printf( "2\n" );
			else printf( "-1\n" );
			continue;
		}
		if( a == 1 ) {
			if( b == 0 ) printf( "-1\n" );
			else printf( "%lld\n", ( ( t - x1 ) % p + p ) % p * inv( b ) % p + 1 );
			continue;
		}
		printf( "%lld\n", bsgs( a, ( t + b * inv( a - 1 ) % p ) * inv( x1 + b * inv( a - 1 ) % p ) % p, p ) );
	}	
	return 0;
}

[BZOJ4128]Matrix

description

solution

因为保证了 A A A矩阵有逆,所以可以选择求矩阵的逆

但其实还是选择两边同乘矩阵要好一些

BSGS模板,有优质代码是写了非常迅速的hash

但是偶直接map里面塞矩阵,重载一下map的小于等于比较挺不戳的!
在这里插入图片描述

code

#include <map>
#include <cmath>
#include <cstdio>
#include <cstring>
using namespace std;
int n, mod;
struct matrix {
	int c[71][71];
	void clear() {
		memset( c, 0, sizeof( c ) );
	}
	matrix() {
		clear();
	}
	void init() {
		for( int i = 1;i <= n;i ++ ) {
			for( int j = 1;j <= n;j ++ )
				c[i][j] = 0;
			c[i][i] = 1;
		}
	}
	matrix operator * ( const matrix &t ) const {
		matrix ans;
		for( int i = 1;i <= n;i ++ )
			for( int j = 1;j <= n;j ++ )
				for( int k = 1;k <= n;k ++ )
					ans.c[i][j] = ( ans.c[i][j] + c[i][k] * t.c[k][j] ) % mod;
		return ans;
	}
	friend bool operator == ( const matrix &s, const matrix &t ) {
		for( int i = 1;i <= n;i ++ )
			for( int j = 1;j <= n;j ++ )
				if( s.c[i][j] ^ t.c[i][j] ) return 0;
		return 1;
	}
	friend bool operator < ( const matrix &s, const matrix &t ) {
		for( int i = 1;i <= n;i ++ )
			for( int j = 1;j <= n;j ++ )
				if( s.c[i][j] ^ t.c[i][j] ) return s.c[i][j] < t.c[i][j];
		return 0;
	}
}A, B, g, h;
map < matrix, int > mp;

int main() {
	scanf( "%d %d", &n, &mod );
	for( int i = 1;i <= n;i ++ )
		for( int j = 1;j <= n;j ++ )
			scanf( "%d", &A.c[i][j] );
	for( int i = 1;i <= n;i ++ )
		for( int j = 1;j <= n;j ++ )
			scanf( "%d", &B.c[i][j] );
	int m = ceil( sqrt( mod * 1.0 ) );
	g.init();
	for( int i = 1;i <= m;i ++ ) {
		g = g * A;
		h = g * B;
		mp[h] = i;
	}
	h.init();
	for( int i = 1;i <= m;i ++ ) {
		h = h * g;
		if( mp.count( h ) ) return ! printf( "%d\n", i * m - mp[h] );
	}
	return 0;
}

CF1106F Lunar New Year and a Recursive Sequence

problem

solution

f i = ∏ j = 1 k f i − j b j ( m o d p ) f_{i}=\prod_{j=1}^kf_{i-j}^{b_j}\pmod p fi=j=1kfijbj(modp)

p p p是非常特殊的质数,其原根为 3 3 3,可以用原根改写上式,令 3 g i = f i 3^{g_i}=f_i 3gi=fi

g i = ∑ j = 1 k g i − j ⋅ b j ( m o d φ ( p ) ) g_i=\sum_{j=1}^kg_{i-j}·b_j\pmod {\varphi(p)} gi=j=1kgijbj(modφ(p))

此时的 g g g可以用矩阵加速,初始化 f i = 0 ( i < k ) , f k = 1 f_i=0(i<k),f_k=1 fi=0(i<k),fk=1,转移 n − k n-k nk次变成 [ f n − k + 1 , . . . , f n ] [f_{n-k+1},...,f_n] [fnk+1,...,fn]
[ h i − k + 1 , . . . , h i ⏟ k ] × [ 0 0 . . . 0 b k 1 0 . . . 0 b k − 1 0 1 . . . 0 b k − 2 . . . . . . . . . . . . . . . 0 0 . . . 0 b 2 0 0 . . . 1 b 1 ] = [ h i + k + 2 , . . . , h i + 1 ⏟ k ] \begin{bmatrix} \underbrace{h_{i-k+1},...,h_i}_{k} \end{bmatrix}\times \begin{bmatrix} 0&0&...&0&b_k\\ 1&0&...&0&b_{k-1}\\ 0&1&...&0&b_{k-2}\\ ...&...&...&...&...\\ 0&0&...&0&b_2\\ 0&0&...&1&b_1 \end{bmatrix}= \begin{bmatrix} \underbrace{h_{i+k+2},...,h_{i+1}}_{k} \end{bmatrix} [k hik+1,...,hi]×010...00001...00..................000...01bkbk1bk2...b2b1=[k hi+k+2,...,hi+1]
同时,有 f k t ≡ f n ( m o d p ) , 3 g k = f k f_k^t\equiv f_n\pmod p,3^{g_k}=f_k fktfn(modp),3gk=fk,这里的 t t t矩等于加速矩阵自乘后的 h k h_k hk

f k T = 3 g k ⋅ T ≡ f n = 3 g n ( m o d p ) f_k^T=3^{g_k·T}\equiv f_n=3^{g_n}\pmod p fkT=3gkTfn=3gn(modp)

对于 3 g n ≡ f n ( m o d p ) 3^{g_n}\equiv f_n\pmod p 3gnfn(modp) g n g_n gn未知,用BSGS求得

对于指数用exgcd计算, g k × t ≡ g n ( m o d p − 1 ) ⇔ t × g k − ( p − 1 ) × y = g n g_k\times t\equiv g_n\pmod {p-1}\Leftrightarrow t\times g_k-(p-1)\times y=g_n gk×tgn(modp1)t×gk(p1)×y=gn

未知数为 g k , y g_k,y gk,y,求出后 f k = 3 g k f_k=3^{g_k} fk=3gk

code

#include <map>
#include <cmath>
#include <cstdio>
#include <cstring>
using namespace std;
#define int long long
#define mod 998244353
struct matrix {
	int n, m;
	int c[101][101];
	matrix() {
		memset( c, 0, sizeof( c ) );
	}
	matrix operator * ( matrix &t ) const {
		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 - 1 );
		return ans;
	}
}h, t;

matrix qkpow( matrix x, int y ) {
	matrix ans;
	ans.n = ans.m = x.n;
	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;
}

map < int, int > mp;

int bsgs( int a, int r ) {
	int m = ceil( sqrt( mod * 1.0 ) );
	int IS = 1;
	for( int i = 1;i <= m;i ++ ) {
		IS = IS * a % mod;
		mp[IS * r % mod] = i;
	}
	int MS = 1;
	for( int i = 1;i <= m;i ++ ) {
		MS = MS * IS % mod;
		if( mp.count( MS ) )
			return i * m - mp[MS];
	}
	return -1;
}

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

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

signed main() {
	int k, n, fn;
	scanf( "%lld", &k );
	h.n = h.m = t.m = k, t.n = t.c[1][k] = 1;
	for( int i = 2;i <= k;i ++ )
		h.c[i][i - 1] = 1;
	for( int i = k;i;i -- )
		scanf( "%lld", &h.c[i][k] );
	scanf( "%lld %lld", &n, &fn );
	h = qkpow( h, n - k );
	t = t * h;
	int g = bsgs( 3, fn );
	int x, y, d = exgcd( t.c[1][k], mod - 1, x, y );
	if( g % d ) return ! printf( "-1\n" );
	else {
		x = ( x * ( g / d ) % ( mod - 1 ) + ( mod - 1 ) ) % ( mod - 1 );
		printf( "%lld\n", qkpow( 3, x ) );
	}
	return 0;
}

Fermat’s Last Theorem

description

solution

找出 p p p的原根 g g g,令 x = g a , y = g b , z = g c x=g^a,y=g^b,z=g^c x=ga,y=gb,z=gc
g a n + g b n ≡ g c n ( m o d p ) g^{an}+g^{bn}\equiv g^{cn}\pmod p gan+gbngcn(modp)
r = g n ⇒ r a + r b ≡ r c ( m o d p ) r=g^n\Rightarrow r^a+r^b\equiv r^c\pmod p r=gnra+rbrc(modp)
假设 a ≤ b a\le b ab

  • a ≤ c ⇒ 1 + r b − a ≡ r c − a ( m o d p ) a\le c\Rightarrow 1+r^{b-a}\equiv r^{c-a}\pmod p ac1+rbarca(modp)
  • a > c ⇒ r a − c + r b − a ≡ 1 ( m o d p ) a>c\Rightarrow r^{a-c}+r^{b-a}\equiv 1\pmod p a>crac+rba1(modp)

枚举 r k r^k rk,出现循环的时候就break
出现 1 + r k 1 ≡ r k 2 / r k 1 + r k 2 ≡ 1 ( m o d p ) 1+r^{k_1}\equiv r^{k_2}/r^{k_1}+r^{k_2}\equiv 1\pmod p 1+rk1rk2/rk1+rk21(modp)就是有解了

code

#include <cstdio>
#define int long long
#define maxn 1000005
int d[maxn], vis[maxn], mi[maxn];

int qkpow( int x, int y, int mod ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

int root( int p ) {
	int cnt = 0, phi = p - 1;
	for( int i = 2;i * i <= phi;i ++ )
		if( phi % i == 0 ) {
			d[++ cnt] = i;
			if( i * i != phi ) d[++ cnt] = phi / i;
		}
	for( int g = 1, i;;g ++ ) {
		for( i = 1;i <= cnt;i ++ )
			if( qkpow( g, d[i], p ) == 1 )
				break;
		if( i > cnt ) return g;
	}
}

signed main() {
	int T, n, p;
	scanf( "%lld", &T );
	for( int Case = 1;Case <= T;Case ++ ) {
		scanf( "%lld %lld", &n, &p );
		int g = root( p ), r = qkpow( g, n, p );
		int t = 1, x = -1, y, z;
		for( int i = 0;i < p;i ++ ) {
			if( i ) t = t * r % p;
			if( vis[t] == Case ) break;
			vis[t] = Case, mi[t] = i;
			int MS = 1 + t;
			if( MS >= p ) MS = 0;
			int IS = 1 - t + p;
			if( IS >= p ) IS = 0;
			if( vis[MS] == Case ) {
				x = 1, y = qkpow( g, i, p ), z = qkpow( g, mi[MS], p );
				break;
			}
			if( vis[IS] == Case ) {
				x = qkpow( g, i, p ), y = qkpow( g, mi[IS], p ), z = 1;
				break;
			}
		}
		if( ~ x ) printf( "%lld %lld %lld\n", x, y, z );
		else printf( "-1\n" );
	}
	return 0;
}
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值