yukicoder : No.613 Solitude by the window

知識ゲー？
https://yukicoder.me/problems/no/613

問題

a(0)=2
a(n+1) = a(n)*(a(n)+4)
とする。
素数Mに対しa(N) mod Mを答えよ。

解法

b(n)=a(n)+2とすると、b(0)=4、b(n+1)=b(n)^2-2となる。
Editorialによると、b(n)の一般項はとなる。らしい。

以後Editorialでは√3の処理をかなり頑張っているが、そこは手を抜いても解くことができる。
いや、自分も平方剰余の相互法則とかフロベニウスとか知らないので…。

整数を返す関数s(k)・t(k)とし、(2+√3)^k = s(k) + t(k)・√3と表せるとする。
よって(2+√3)^(k+1)=(2s(k)+3t(k)) + (s(k)+2t(k))√3より、s(k+1)=2s(k)+3t(k)、t(k+1)=s(k)+2t(k)と置ける。

さて、この(2+√3)^kだが、kに対し(M-1)周期になるかと思えばそうではなく(M+1)または(M-1)周期のいずれかになる(このあたりはEditorial参照。自分も深く理解していない)。
面倒なので(M^2-1)周期と考えればどちら周期でも問題ない。
(2-√3)^k=s(k)-t(k)√3であることを考えると、よりb(n)が求められる。

ll N;
ll mo;
map<int,int> mp;

const int MAT=2;
struct Mat { ll v[MAT][MAT]; };
Mat mulmat(Mat& a,Mat& b,int n=MAT) {
	ll mo2=4*mo*mo;
	int x,y,z; Mat r;
	FOR(x,n) FOR(y,n) r.v[x][y]=0;
	FOR(x,n) FOR(z,n) FOR(y,n) {
		r.v[x][y] += a.v[x][z]*b.v[z][y];
		if(r.v[x][y]>mo2) r.v[x][y] -= mo2;
	}
	FOR(x,n) FOR(y,n) r.v[x][y]%=mo;
	return r;
}

Mat powmat(ll p,Mat a,int n=MAT) {
	int i,x,y; Mat r;
	FOR(x,n) FOR(y,n) r.v[x][y]=0;
	FOR(i,n) r.v[i][i]=1;
	while(p) {
		if(p%2) r=mulmat(r,a,n);
		a=mulmat(a,a,n);
		p>>=1;
	}
	return r;
}

ll modpow(__int128_t a, ll n,ll mo) {
	
	__int128_t r=1;
	while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1;
	return r;
}

void solve() {
	int i,j,k,l,r,x,y; string s;
	
	cin>>N>>mo;
	if(mo==2) return _P("0\n");
	
	Mat A={ .v={{2,3},{1,2}}};
	Mat P=powmat(modpow(2,N,mo*mo-1),A);
	
	cout<<(P.v[0][0]*2+mo-2)%mo<<endl;
}

まとめ

いや、b(n)の一般項がわからない時点でどうしようもないのよ…。