1. 程式人生 > 實用技巧 >[NOI2012]隨機數生成器【矩陣快速冪】

[NOI2012]隨機數生成器【矩陣快速冪】

NOI2012 隨機數生成器

題目描述

棟棟最近迷上了隨機演算法,而隨機數是生成隨機演算法的基礎。棟棟準備使用線性同餘法(Linear Congruential Method)來生成一個隨機數列,這種方法需要設定四個非負整數引數 \(m,a,c,X_0\),按照下面的公式生成出一系列隨機數 \(\{X_n\}\)

\[X_{n+1}=(aX_n +c)\bmod m \]

其中\(mod\ m\) 表示前面的數除以 \(m\) 的餘數。從這個式子可以看出,這個序列的下一個數總是由上一個數生成的。

棟棟知道這樣產生的序列具有良好的隨機性,不過心急的他仍然想盡快知道 \(X_n\) 是多少。由於棟棟需要的隨機數是 \(0,1,\dots,g-1\)

之間的,他需要將 \(X_n\)​ 除以 \(g\) 取餘得到他想要的數,即 \(X_n \bmod g\),你只需要告訴棟棟他想要的數 \(X_n \bmod g\) 是多少就可以了。

輸入格式

一行 \(6\) 個用空格分割的整數 \(m,a,c,X_0,n\)\(g\),其中 \(a,c,X_0\) 是非負整數,\(m,n,g\) 是正整數。

輸出格式

輸出一個數,即 \(X_n \bmod g\)

輸入輸出樣例

輸入

11 8 7 1 5 3

輸出

2

說明/提示

計算得 \(X_n=X_5=8\),故\((X_n \bmod g) = (8 \bmod 3) = 2\)

對於 \(100\%\) 的資料,\(n,m,a,c,X_0\leq 10^{18}\)\(1\leq g\leq 10^8\)\(n,m\geq 1\)\(a,c,X_0\geq 0\)

題意

給出了一個迷惑式子,讓你算出來式子的第\(n\)項,然後\(mod\ g\)的結果

分析

看到這樣一個個的遞推式子,一個個用\(for\)迴圈來推肯定不行,所以很容易就會想到要用到矩陣快速冪來求。那麼我們現在的主要任務就是構造矩陣來進行乘法運算。

首先看到題目中給出的式子:

\[X_{n+1}=(aX_n+c)\ mod\ m \]

取模運算可以暫且先不看,因為對結果沒什麼影響,在矩陣乘法的時候進行取模就行了。所以轉化成如下式子:

\[X_{n+1}=aX_n+c \]

那麼我們就可以根據這個式子來構造矩陣。由矩陣的乘法運算為結果矩陣的\(i\)\(j\)列為前邊矩陣一個的第\(i\)行乘以另一個的第\(j\)列,所以我們可以得出如下的矩陣遞推式子:

\[\left[ \begin{matrix} X_{n-1}\\ c \end{matrix} \right]\times \left[ \begin{matrix} a & 1\\ 0 & 1 \end{matrix} \right] = \left[ \begin{matrix} X_{n}\\ c \end{matrix} \right] \]

這裡用\(X_{n-1}\)這一列分別乘以右邊矩陣的第一第二行,得到結果的矩陣,那麼我們就可以根據這個遞推式子來進行矩陣快速冪。
這裡乘法的運算過程如下:

\[X_{n-1}\times a+c\times 1 = X_n \]

\[X_{n-1}\times 0+c\times 1 = c \]

由此得到結果矩陣

這裡的矩陣做乘法的時候需要用到龜速乘,不然會爆\(long\ long\)

程式碼

#include<bits/stdc++.h>
using namespace std;
#define int long long
struct Node{//矩陣結構體
	int a[5][5];
};
int m,a,c,x0,g,n;
int ksj(int a,int b){//龜速乘
	int ans = 0;
	while(b){
		if(b & 1)ans = (ans + a) % m;
		a = (a + a) % m;
		b >>= 1;
	}
	return ans;
}
Node Mul(Node a,Node b,int c){//矩陣乘法,記得取模
	 Node ans;
	 memset(ans.a,0,sizeof(ans.a));
	 for(int i=1;i<=2;++i){
		 for(int j=1;j<=2;++j){
			 for(int k=1;k<=2;++k){
				 ans.a[i][j] = (ans.a[i][j] + ksj(a.a[i][k],b.a[k][j])%c)%c;
			 }
		 }
	 }
	 return ans;
}
Node ans;
void qpow(Node &ans,Node b,int c){//矩陣快速冪
	while(c){
		if(c & 1)ans = Mul(b,ans,m);
		b = Mul(b,b,m);
		c >>= 1;
	}
}
signed main(){
	Node bas;
	scanf("%lld%lld%lld%lld%lld%lld",&m,&a,&c,&x0,&n,&g);
	ans.a[1][1] = x0;//初始化矩陣
	ans.a[2][1] = c;
	bas.a[1][1] = a;
	bas.a[1][2] = 1;
	bas.a[2][1] = 0;
	bas.a[2][2] = 1;
	qpow(ans,bas,n);
	printf("%lld\n",ans.a[1][1] % g);//得答案
}