快速幂 – Qizy's Database

【日常小测】题目1

2017年3月1日2017年3月27日 Qizy 2 Comments

题目大意

求$\sum_{i=1}^n{f_i^k}$
其中$f_x$为第$x$项$Fibonacci$数
$n \le 1e18, k \le 1e5$

解题报告

之前鏼爷讲过二次剩余，于是看到这个模数就知道方向了
在暴力求出$\sqrt 5$的二次剩余后，我们可以推出答案长这样
$$\sum_{j=0}^{k}{(-1)^{k-j} \cdot C_k^j \sum_{i-1}^n{(A^jB^{k-j})^i}}$$
于是我们搞一搞组合数，逆元什么的这题就做完啦！

Code

#include<bits/stdc++.h>
#define LL long long
using namespace std;

const int MOD = 1000000009;
const int A = 691504013;
const int B = 308495997;
const int REV_SQRT_FIVE = 276601605;
const int N = 100009;

int k,vout,PA[N],PB[N];

inline int read() {
	char c=getchar(); int f=1,ret=0;
	while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
	return ret * f;
}

inline int Pow(int w, LL t) {
	int ret = 1;
	for (;t;t>>=1,w=(LL)w*w%MOD) 
		if (t & 1) ret = (LL)ret * w % MOD;
	return ret;
}

inline int Mul(int a, LL b) {
	int ret = 0;
	for (;b;b>>=1,(a<<=1)%=MOD)
		if (b & 1) (ret += a) %= MOD;
	return ret;
}

int main() {
	freopen("A.in","r",stdin); 
	freopen("A.out","w",stdout); 
	LL n; cin>>n; k = read();
	PA[0] = PB[0] = 1;
	for (int i=1;i<=k;i++) {
		PA[i] = (LL)PA[i-1] * A % MOD;
		PB[i] = (LL)PB[i-1] * B % MOD;
	}
	for (int i=0,c=1,v;i<=k;i++) {
		v = (LL)PA[i] * PB[k-i] % MOD; 
		if (v == 1) v = Mul(v, n);
		else v = ((LL)Pow(v, n+1) - v) * Pow(v-1, MOD-2) % MOD;
		if (k-i & 1) (vout -= (LL)c * v % MOD) %= MOD;
		else (vout += (LL)c * v % MOD) %= MOD;
		c = ((LL)c * (k - i) % MOD) * Pow(i+1, MOD-2) % MOD; 
	} 
	printf("%d\n",((LL)vout*Pow(REV_SQRT_FIVE,k)%MOD+MOD)%MOD);
	return 0;
}

用全概率公式可以推得$p(i) = – \frac{{1 – p[i – 1]}}{2}$
然后推一推通项公式可以得到$p(i) = \frac{{{{( – 1)}^i} + {2^{i – 1}}}}{{3 \cdot {2^{i – 1}}}}$
然后我们惊奇的发现${{{( – 1)}^i} + {2^{i – 1}}}$总可以被3整除
而${{{( – 1)}^i} + {2^{i – 1}}}$又是一个奇数，与2的整次幂肯定不能约分
于是搞一搞快速幂什么的就可以了

Code

#include<bits/stdc++.h>
#define LL long long
using namespace std;

const int MOD = 1000000007;
const int REV = 333333336;
const int N = 100000+9;

LL k,arr[N];

inline int pow(LL w, LL t) {
	int ret = 1; while (t) {
		if (t & 1) ret = (LL)ret*w % MOD;
		w = w*w % MOD; t >>= 1;
	} return ret;
}

inline LL read(){
	char c=getchar(); LL ret=0,f=1;
	while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
	return ret*f;
}

int main(){
	k = read(); int tag = 1;
	for (int i=1;i<=k;i++) tag = ((arr[i] = read()) == 1 && tag);
	if (tag) puts("0/1");
	else {
		LL numerator=REV,denominator=2; tag = -1;
		for (int i=1;i<=k;i++) if (arr[i]%2 == 0) tag = 1;
		for (int i=1;i<=k;i++) denominator = pow(denominator,arr[i]);
		(denominator *= pow(2LL,MOD-2LL)) %= MOD;
		(numerator *= denominator + tag) %= MOD;
		printf("%d/%d",(int)numerator,(int)denominator);
	}
	return 0;
}

【BZOJ 4522】[Cqoi2016] 密钥破解

2016年9月6日2016年9月6日 Qizy 3 Comments

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4522
数据生成器：http://paste.ubuntu.com/23140343/

这题，我觉得唯一的难点在分解因数
因为忘了棒棒糖怎么写，于是考试的时候只能写了50分的暴力
话说CQOI怎么能考模板呢？

#include<bits/stdc++.h>
#define LL long long
#define abs(x) ((x)>0?(x):-(x))
using namespace std;

LL e,N,c,r,d,n; 

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

inline LL EX_GCD(LL a, LL b) {
	LL t1, t2;
	EX_GCD(t1,a,t2,b);
	return (t1%r + r) % r;
}

inline LL Mul(LL a, LL b) {
	LL ret = 0;
	while (b) {
		if (b&1) ret  = (ret + a) % N;
		a = a*2 % N; b >>= 1;
	}
	return ret;
}

inline LL pow(LL w, LL t) {
	LL ret = 1;
	while (t) {
		if (t & 1) ret = Mul(ret, w);
		w = Mul(w, w); t >>= 1;
	}
	return ret;
}

inline LL GCD(LL a, LL b){
	while (b) a = a % b, swap(a, b);
	return a;
}

inline void Decompose(){
	for (int seed=(rand()&32767)%(N-1)+1;!r;seed=(rand()&32767)%(N-1)+1) {
		LL w = rand() % (N-1) + 1, seg=2, tmp = -1, cnt = 0, gcd;
		while (w != tmp && !r) {
			gcd = GCD(abs(w-tmp), N);
			if (1 < gcd && gcd < N) r = Mul(N/gcd-1,gcd-1); 
			if (++cnt == seg) seg <<= 1, tmp = w;
			w = (Mul(w,w) + seed) % (N-1) + 1;
		}
	}		
}

int main(){
	srand(999971); 
	cin>>e>>N>>c; Decompose(); 
	d = EX_GCD(e,r); n = pow(c,d);
	cout<<d<<' '<<n;
	return 0;
}

【UOJ 241】破坏发射台

2016年8月31日2017年1月12日 Qizy 21 Comments

题目传送门：http://uoj.ac/problem/241
官方题解：http://c-sunshine.blog.uoj.ac/blog/2026

考试的时候，推这个题推了两个半小时QAQ
一直觉得马上就可以推出来辣！结果一直没有推出来
不过最后还是成功把奇数的式子推出来了（虽然和std不同，但是正确的）
至于偶数的式子嘛，感觉题解说的不是很清楚，我也就来哔哔两句：
不难发现，如果我们仍然按位来递推，不方便处理穿过发射台的情况
于是我们定义二元组（f，g）表示第i位与i+n/2位的颜色情况
于是我们发现这样的定义就可以按二元组为阶段来递推了，因为这货就无后效性了
于是搞一搞矩阵乘法即可
ps：为什么我做这题的时候老想着数位DP

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#define LL long long
#define abs(x) ((x)>0?(x):-(x))
#define S(x,y) ((x)*3+(y)+1)
using namespace std;

const int MOD = 998244353;
const int N = 100000;
const int M = 9;

int n,m,vout;

struct Matrix{
	int a[M][M];
	inline Matrix() {}
	inline Matrix(int v) {memset(a,0,sizeof(a));for(int i=1;i<M;i++) a[i][i]=v;}
	inline Matrix operator = (const Matrix &B) {memcpy(&(*this),&B,sizeof(a));}
	inline Matrix operator * (const Matrix &B) {
		Matrix ret(0);
		for (int i=1;i<M;i++) for (int j=1;j<M;j++) for (int k=1;k<M;k++) 
			ret.a[i][j] = ((LL)a[k][j]*B.a[i][k] + ret.a[i][j]) % MOD;
		return ret;
	}
	inline Matrix operator ^ (int t) {
		Matrix ret(1),w=*this;
		while (t) {
			if (t & 1) ret = ret*w;
			w = w*w; t >>= 1;	
		}
		return ret;
	}
};

inline int read(){
	char c=getchar(); int ret=0,f=1;
	while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0') ret=ret*10+c-'0',c=getchar();
	return ret*f;
}

inline int pow(int w, int t){
	int ret = 1;
	while (t) {
		if (t & 1) ret = (LL)ret*w % MOD;
		w = (LL)w*w % MOD; t >>= 1;
	}
	return ret;
} 

int main(){
	n = read(), m = read();
	if (n%2 == 0) {
		Matrix ori(0),tra(0); ori.a[S(1,2)][1]=1;
		if (m > 3) tra.a[S(0,0)][S(1,2)] = tra.a[S(0,0)][S(2,1)] = (LL)(m - 2) * (m - 3) % MOD;
		if (m > 3) tra.a[S(0,0)][S(1,0)] = tra.a[S(0,0)][S(0,1)] = (LL)(m - 3) * (m - 3) % MOD;
		if (m > 3) tra.a[S(0,0)][S(2,0)] = tra.a[S(0,0)][S(0,2)] = (LL)(m - 3) * (m - 3) % MOD;
		if (m > 3) tra.a[S(0,0)][S(0,0)] = ((LL)(m-4) * (m-4) + m - 3) % MOD;
		
		if (m > 2) tra.a[S(1,0)][S(2,1)] = tra.a[S(0,1)][S(1,2)] = m - 2;
		if (m > 2) tra.a[S(1,0)][S(0,1)] = tra.a[S(0,1)][S(1,0)] = m - 2;
		if (m > 2) tra.a[S(1,0)][S(0,2)] = tra.a[S(0,1)][S(2,0)] = m - 2;
		if (m > 3) tra.a[S(1,0)][S(2,0)] = tra.a[S(0,1)][S(0,2)] = m - 3;
		if (m > 3) tra.a[S(1,0)][S(0,0)] = tra.a[S(0,1)][S(0,0)] = m - 3;
		
		if (m > 2) tra.a[S(2,0)][S(1,2)] = tra.a[S(0,2)][S(2,1)] = m - 2;
		if (m > 2) tra.a[S(2,0)][S(0,2)] = tra.a[S(0,2)][S(2,0)] = m - 2;
		if (m > 2) tra.a[S(2,0)][S(0,1)] = tra.a[S(0,2)][S(1,0)] = m - 2;
		if (m > 3) tra.a[S(2,0)][S(1,0)] = tra.a[S(0,2)][S(0,1)] = m - 3;
		if (m > 3) tra.a[S(2,0)][S(0,0)] = tra.a[S(0,2)][S(0,0)] = m - 3;
		
		tra.a[S(1,2)][S(0,1)] = tra.a[S(2,1)][S(1,0)] = 1;
		tra.a[S(1,2)][S(2,0)] = tra.a[S(2,1)][S(0,2)] = 1;
		tra.a[S(1,2)][S(2,1)] = tra.a[S(2,1)][S(1,2)] = 1;
		tra.a[S(1,2)][S(0,0)] = tra.a[S(2,1)][S(0,0)] = 1;
		
		tra = tra^(n/2-1); ori = ori * tra;
		vout = ((LL)ori.a[S(0,2)][1] + ori.a[S(1,0)][1] + ori.a[S(1,2)][1] + ori.a[S(0,0)][1]) % MOD;
		vout = ((LL)vout * m % MOD) * (m-1) % MOD;
		printf("%d\n",vout);
	}
	else {
		if (n == 3) printf("%d\n",((LL)(pow(m-1,n-1)-(m-1))*m % MOD + MOD) % MOD);
		else {
			Matrix ori(0); ori.a[1][1] = 1; ori.a[2][1] = m-2;
			Matrix tra(0); tra.a[1][2] = 1; tra.a[2][1] = m-1; tra.a[2][2] = m-2; 
			tra = tra^(n-5); ori = ori * tra; ori.a[1][1] = (LL)ori.a[1][1]*(m-1)%MOD;
			vout = (((LL)(m-1)*(m-2)%MOD)*ori.a[2][1]%MOD + (LL)ori.a[1][1]*(m-1)%MOD) % MOD;
			vout = (LL)(pow(m-1,n-1)-vout)*m % MOD;
			printf("%d\n",(vout+MOD)%MOD);
		}
	}
	return 0;
}

—– UPD 2017.1.12 —–
今天听毛爷爷讲这个题
居然感觉自己只会做奇数情况 QwQ
其实重点是，这货偶数情况的转移似乎非常神奇的样子
定义状态的方式，非常值得借鉴

【BZOJ 2242】[SDOI2011] 计算器

2016年8月10日2016年8月27日 Qizy 236 Comments

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=2242
数据生成器：http://paste.ubuntu.com/22875643/
数据传送门：http://pan.baidu.com/s/1dE71QHr

这个题，做了一个下午+一个晚上，最后才发现是把”cannot”中间多加了一个‘a’
MD，老子好想杀人！ (╯‵□′)╯︵┻━┻

题目本身很简单，就是快速幂+拓展欧几里得+BSGS
但网上的程序，多少有问题，反正我能搜出来的，我都可以叉掉
我这份代码问题应该比较少吧：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#define LL long long
using namespace std;

const int MAXN = 100000;
const double EPS = 1e-5;

namespace Hash{
	const int MOD = 99971;
	int head[MAXN],nxt[MAXN],val[MAXN],data[MAXN],T;
	
	inline void init(){
		T = 0;
		memset(head,0,sizeof(head));
	}
	
	inline void insert(int w, int d){
		nxt[++T] = head[w%MOD];
		val[T] = w; data[T] = d;
		head[w%MOD] = T;
	}	
	
	inline int find(int w){
		for (int i=head[w%MOD];i;i=nxt[i])
			if (val[i] == w) return data[i];
		return -1;
	}
};

inline int read(){
	char c=getchar(); int ret=0,f=1;
	while (c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0'){ret=ret*10+c-'0';c=getchar();}
	return ret*f;
}

inline int quick_pow(int a,int b,int c){
	int ret = 1; while (b) {
		if (b & 1) ret = ((LL)ret * a) % c;
		a = (LL)a*a%c; b >>= 1;
	} return ret;
}

int GCD(int a, int b){return b?GCD(b,a%b):a;}

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

inline int EX_GCD(int a, int b, int c) {
	int ret,tmp,gcd=GCD(a,b); 
	if (c%gcd) return -1;
	else {
		EX_GCD(a/gcd,b/gcd,ret,tmp);
		ret = (((LL)c/gcd*ret) % b + b) % b;
		return ret;
	}
}

inline void solve(int a, int b, int c) {
	int ret = EX_GCD(a,c,b); 
	if (~ret) printf("%d\n",ret);
	else printf("Orz, I cannot find x!\n");
}

inline void BSGS(int a, int b, int c) {
	Hash::init(); a %= c; b %= c;
	if (!a && b > 1) {printf("Orz, I cannot find x!\n"); return;}
	
	int lim = int(ceil(sqrt(c))+EPS);
	for (int i=0,w=1;i<=lim;i++) {
		if (w == b) {printf("%d\n",i); return;}
		Hash::insert(w,i);
		w = ((LL)w*a) % c;
	}
	
	int w_ori = quick_pow(a,lim,c), rev_ori = quick_pow(w_ori,c-2,c);
	for (int i=1,w=w_ori,rev=rev_ori;i<=lim;i++) {
		int tmp = Hash::find(((LL)rev*b) % c);
		if (tmp > -1) {printf("%d\n",i*lim+tmp); return;}
		w = ((LL)w*w_ori) % c;
		rev = ((LL)rev*rev_ori) % c;
	}
	printf("Orz, I cannot find x!\n");
}

int main(){
	int T=read(),ty=read(); while (T--) {
		int a=read(), b=read(), c=read();
		if (ty == 1) printf("%d\n",quick_pow(a,b,c));
		else if (ty == 2) solve(a,b,c);
		else BSGS(a,b,c);
	}
	return 0;
}

【算法笔记】快速幂与快速乘

2016年6月27日2016年6月27日 Qizy 1 Comment

quick_pow_and_quick_mul

PS:貌似也有人喜欢把快速乘叫成慢速乘？~~比如我~~

Tag: 快速幂

【日常小测】题目1

相关链接

题目大意

解题报告

Code

【Codeforces 696C】PLEASE

相关链接

解题报告

Code

【BZOJ 4522】[Cqoi2016] 密钥破解

【UOJ 241】破坏发射台

【BZOJ 2242】[SDOI2011] 计算器

【算法笔记】快速幂与快速乘