Xài thẳng công thức
https://rextester.com/RYA77507
#include <boost/multiprecision/gmp.hpp>
#include <iostream>
template <class T, size_t SqrtIntegerValue>
struct SqrtSum {
T rp = T(0); // rational part
T ip = T(0); // irrational part
SqrtSum(T rp, T ip = T(0)) : rp(rp), ip(ip) {}
SqrtSum operator+(const SqrtSum& rhs) const
{
return {rp + rhs.rp, ip + rhs.ip};
}
SqrtSum& operator*=(const SqrtSum& rhs)
{
T temp = rp * rhs.rp + ip * rhs.ip * SqrtIntegerValue;
ip = rp * rhs.ip + ip * rhs.rp;
rp = temp;
return *this;
}
SqrtSum operator*(SqrtSum rhs) const { return rhs *= *this; }
};
template <class T, size_t V>
auto pow(SqrtSum<T, V> a, size_t p)
{
SqrtSum<T, V> res = T(1);
for (; p; p >>= 1, a *= a)
if (p & 1)
res *= a;
return res;
}
int main()
{
using rational_t = boost::multiprecision::mpq_rational;
auto fNth = [](auto u, auto v, auto a, auto b) {
return [=](size_t n) { return (a * pow(u, n) + b * pow(v, n)).rp; };
};
auto f = [&](size_t n) {
using S2 = SqrtSum<rational_t, 2>;
return fNth(S2(1, 1), S2(1, -1), S2({2, 4}, {1, 4}),
S2({2, 4}, {-1, 4}))(n);
};
for (size_t n = 0; n <= 25; ++n)
std::cout << n << "\t" << f(n) << "\n";
auto fibo = [&](size_t n) {
using S5 = SqrtSum<rational_t, 5>;
return fNth(S5({1, 2}, {1, 2}), S5({1, 2}, {-1, 2}), S5(0, {1, 5}),
S5(0, {-1, 5}))(n);
};
for (size_t n = 0; n <= 25; ++n)
std::cout << n << "\t" << fibo(n) << "\n";
}