这道题目是练习比赛中的一道题：
计算第N个既是三角形数又是平方数的数，对10006699取模。(1 ≤ N ≤ 10^18) 最多有10^5组测试数据。
我发现可以使用递推关系式Ti = 6Ti-1 - Ti-2 + 2 来轻松计算，其中T0 = 0 且 T1 = 1。
我正在使用矩阵指数运算来实现每个测试用例大约O(log N)的性能，但很明显太慢了，因为有10^5个测试用例。事实上，即使在限制条件仅为(1 ≤ N ≤ 10^6)的情况下，这段代码也太慢了，我可以只做O(N)的预处理和O(1)的查询。
我应该改变解决问题的方法还是只优化代码的一些部分呢？
#include <ios>
#include <iostream>
#include <vector>
#define MOD 10006699
/*
Transformation Matrix:
0 1 0 t[i] t[i+1]
-1 6 1 * t[i+1] = t[i+2]
0 0 1 2 2
*/
std::vector<std::vector<long long int> > multi(std::vector<std::vector<long long int> > a, std::vector<std::vector<long long int> > b)
{
std::vector<std::vector<long long int> > c(3, std::vector<long long int>(3));
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
for (int k = 0; k < 3; k++)
{
c[i][j] += (a[i][k] * b[k][j]) % MOD;
c[i][j] %= MOD;
}
}
}
return c;
}
std::vector<std::vector<long long int> > power(std::vector<std::vector<long long int> > vec, long long int p)
{
if (p == 1) return vec;
else if (p % 2 == 1) return multi(vec, power(vec, p-1));
else
{
std::vector<std::vector<long long int> > x = power(vec, p/2);
return multi(x, x);
}
}
int main()
{
std::ios_base::sync_with_stdio(false);
long long int n;
while (std::cin >> n)
{
if (n == 0) break;
else
{
std::vector<std::vector<long long int> > trans;
long long int ans;
trans.resize(3);
trans[0].push_back(0);
trans[0].push_back(1);
trans[0].push_back(0);
trans[1].push_back(-1);
trans[1].push_back(6);
trans[1].push_back(1);
trans[2].push_back(0);
trans[2].push_back(0);
trans[2].push_back(1);
trans = power(trans, n);
ans = (trans[0][1]%MOD + (2*trans[0][2])%MOD)%MOD;
if (ans < 0) ans += MOD;
std::cout << ans << std::endl;
}
}
}