「NOI2012」随机数生成器

知识点: 矩阵加速数列

原题面 Loj [Luogu](https://www.luogu.com.cn/problem/P2044

简述

对于数列 \(X\)，有：

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

给定参数 \(n, m, X_0, a, c, g\)，求 \(X_n\bmod g\)。
\(1<n,m\le 10^{18}\)，\(0\le a,c,X_0\le 10^{18}\)，\(1\le g\le 10^8\)。
1S，512MB。

分析

数据范围这么大，递推式都给了，考虑矩阵加速数列。
有：

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

直接做就可以了，注意 \(m\le 10^{18}\)，开 ull，矩阵乘法中要写龟速乘。

代码

//知识点：矩阵加速数列
/*
By:Luckyblock
*/
#include <algorithm>
#include <cstdio>
#include <ctype.h>
#include <cstring>
#include <iostream>
#define ull unsigned long long
const int kMaxm = 2;
//=============================================================
struct Matrix {
  ull a[3][3];
} ori, ans;
ull m, a, c, x0, n, g;
//=============================================================
inline ull read() {
  ull f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void GetMax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void GetMin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
ull QuickProd(ull x_, ull y_) {
  ull ret = 0;
  while (y_) {
    if (y_ & 1) ret = (ret + x_) % m;
    x_ = (x_ + x_) % m, y_ >>= 1ll;
  }
  return ret;
}
Matrix operator * (const Matrix &a_, const Matrix &b_) {
  Matrix c;
	memset(c.a, 0, sizeof(c.a));
	for (int k = 1; k <= kMaxm; ++ k) {
	  for (int i = 1; i <= kMaxm; ++ i) {
	    for (int j = 1; j <= kMaxm; ++ j) {
	      c.a[i][j] = (c.a[i][j] + QuickProd(a_.a[i][k], b_.a[k][j])) % m;
      } 
    }
  }
	return c;
}
void Build(Matrix &x_) {
  memset(x_.a, 0, sizeof (x_.a));
  for (int i = 1; i <= kMaxm; ++ i) {
    x_.a[i][i] = 1ll;
  }
}
Matrix QuickPow(Matrix x_, ull y_) {
  Matrix ret;
  Build(ret);
  while (y_) {
    if (y_ & 1) ret = ret * x_;
    x_ = x_ * x_, y_ >>= 1ll;
  } 
  return ret;
}
//=============================================================
int main() {
  m = read(), a = read(), c = read(); 
  x0 = read(), n = read(), g = read();
  ori.a[1][1] = a, ori.a[2][1] = c;
  ori.a[1][2] = 0, ori.a[2][2] = 1;
  ans.a[1][1] = x0, ans.a[1][2] = 1;
  ans = ans * QuickPow(ori, n);
  std :: cout << ans.a[1][1] % g;
  return 0;
}