JZOJ 6481 黎曼几何

首先考虑 DP，设 fn,gn 分别表示把 n 个硬币移动到下一个圈 / 下下一个圈的最小步数。

若要把所有硬币移动到下一个圈，首先要把面值最大的移动到下一个圈，然后再把剩下的移到这个圈。
为了将面值最大的移动到下一个圈，要先把剩下的先移到下下个圈，然后再把最大的移到下一个圈。
显然这样可以构建起 fn 的 DP 转移方程： fn=2gn−1+1

若要把所有硬币移动到下下个圈，首先要把面海最大的移动到下下个圈，然后再把剩下的移到这个圈。
为了将面值最大的移动到下下个圈，要先把剩下的先移到下下个圈，然后再把最大的移到下一个圈。
然后再将剩下的移回原来的圈，再将最大的移动到下下个圈，然后再把剩下的移到下下个圈。
总结起来就是 gn=2gn−1+fn−1+2

（我也知道很啰嗦不过这些我相信读者自己能推出来）

然后这里正解是直接矩阵快速幂 + 分块预处理……
我选择推通项（

设 F,G 为 f,g 的生成函数，则 F(x)=2xG(x)+11−x−1G(x)=2xG(x)+xF(x)+21−x−2

用第一条式子替换掉第二条式子中的 F(x)，有 G(x)=2xG(x)+x(2xG(x)+11−x−1)+21−x−2=2x2G(x)+2xG(x)+x+21−x−x−2(1−2x−2x2)G(x)=x1−x(x+2)G(x)=x(1−x)[1−(1+√3)x][1−(1−√3)x](x+2)G(x)=16(x+2)[1+√31−(1+√3)x+1−√31−(1−√3)x]gn=(1+√3)n(12+13√3)+(1−√3)n(12−13√3)−1

（省去了因式分解部分分式等一系列步骤）

然而，可惜的是！3 不是模 998244353 的二次剩余！
不过我们依然可以扩域，即把每个数表示为 a+b√3 的形式，然后分块预处理幂。

代码：

#include <cstdio>
#include <utility>
#define add(x,y) (x + y >= mod ? x + y - mod : x + y)
#define dec(x,y) (x < y ? x - y + mod : x - y)
using namespace std;
const int W = 1 << 20;
const int mod = 998244353;
const int inv[] = {499122177,332748118};
typedef pair<int,int> cp;
inline cp operator+(const cp &a,const cp &b)
{
    return cp(add(a.first,b.first),add(a.second,b.second));
}
inline cp operator-(const cp &a,const cp &b)
{
    return cp(dec(a.first,b.first),dec(a.second,b.second));
}
inline cp operator*(const cp &a,const cp &b)
{
    return cp(((long long)a.first * b.first % mod + 3LL * a.second % mod * b.second % mod) % mod,((long long)a.first * b.second % mod + (long long)a.second * b.first % mod) % mod);
}
cp pw1[2][W + 5],pw2[2][W + 5];
inline cp query1(long long n)
{
    return pw1[0][n & (W - 1)] * pw1[1][n >> 20];
}
inline cp query2(long long n)
{
    return pw2[0][n & (W - 1)] * pw2[1][n >> 20];
}
inline int query_g(long long n)
{
    if(n <= 0)
        return 0;
    int ret = (query1(n) * cp(inv[0],inv[1]) + query2(n) * cp(inv[0],mod - inv[1])).first;
    return ret = dec(ret,1);
}
inline int query_f(long long n)
{
    if(n <= 0)
        return 0;
    int ret = query_g(n - 1);
    ret = add(ret,ret);
    return ret = add(ret,1);
}
int T,ans1,ans2;
long long n;
int main()
{
    freopen("riemannian.in","r",stdin),freopen("riemannian.out","w",stdout);
    pw1[0][0] = pw1[1][0] = pw2[0][0] = pw2[1][0] = cp(1,0);
    for(register int i = 1;i <= W;++i)
        pw1[0][i] = pw1[0][i - 1] * cp(1,1),
        pw2[0][i] = pw2[0][i - 1] * cp(1,mod - 1);
    for(register int i = 1;i < W;++i)
        pw1[1][i] = pw1[1][i - 1] * pw1[0][W],
        pw2[1][i] = pw2[1][i - 1] * pw2[0][W];
    scanf("%d",&T);
    for(;T;--T)
        scanf("%lld",&n),ans1 ^= query_f(n),ans2 ^= query_g(n);
    printf("%d %d\n",ans1,ans2);
}