Вычисление длинных математических последовательностей

Приветствую. Возник вопрос, насчет вычисления последовательности схожей с последовательностью чисел Фибоначчи, только в которой каждый член равен сумме трех предыдущих. Последовательность начинается с трех единиц: 1, 1, 1, 3, 5, 9, .... Я хочу находить число этой последовательности по номеру. Номер лежит в пределах от 1 до 10000. Я изучил последовательность Фибоначчи и нашел несколько формул для вычисления n-ого члена, в частности, формулу Бине(через золотое сечение), и другие. Но я не знаю, как их применить к моей последовательности.
Вопрос: как получить формулу для вычисления n-ого члена последовательности без нахождения всех предыдущих членов? Очевидно, что рекурсивный подход для n=10000 не подходит. UPD: время выполнения программы для любых n не должно превышать 1 секунду.

Ответ

Данный рекурсивный алгоритм нужно переписать итеративно. В этом случае сложность стает O(n), а памяти нужно будет константно (ну почти - числа то увеличиваются).
Вот пример, как посчитать и вывести
#include #include #include #include
using namespace std;
mpz_class pp(int n) { if (n < 1 || n >= 10000000) return 0; if (n < 3) return 1; mpz_class a = 1; mpz_class b = 1; mpz_class c = 1; for (int i = 3; i < n; i++) { mpz_class t = a + b + c; c = b; b = a; a = t; } return a; }
int main() { auto t1 = std::chrono::steady_clock::now();
cout << pp(1000000) << endl;
auto t2 = std::chrono::steady_clock::now(); auto d_milli = std::chrono::duration_cast( t2 - t1 ).count(); cout << d_milli << endl; return 0; }
Вопрос в том, насколько это быстро? Для миллионного элемента на моей машине потребовалось 22секунды. Для 100000 элемента порядка 145мс. Для 10000 элемента это 2-4 мс.
И на последок. Как это скомпилировать (на линуксе):
g++ -std=c++11 ff.cpp -lgmpxx -lgmp
P.S. нужно не забыть поставить библиотеку gmp.
upd
решил переписать с ручной реализацией сложения больших чисел. Код сложения не самый оптимальный, но gmp проигрывает всего лишь 6-8 раз. Если переписать с sse, думаю, можно будет вытянуть больше. Но даже для такой простой реализации в лоб, думаю уже не плохо.
#include #include #include #include
#include #include #include
#define CLEN 8 #define BASE 100000000
class bi { public: bi() { } bi(std::string data) { for (int i = 0; i < data.length(); i+= CLEN) { int r = std::stoi(data.substr(i, CLEN)); m_digits.push_back(r); } } bi(const bi& data) { m_digits = data.m_digits; }
bi(int data) { while (data > 0) { m_digits.push_back(data % BASE); data = data / BASE; } } bi& operator=(bi data) { if (this != &data) { m_digits = data.m_digits; } return *this; } friend std::ostream& operator<< (std::ostream& stream, const bi& b);
friend bi operator+( const bi& lhs, const bi& rhs ); private: std::vector m_digits;
};
std::ostream& operator<< (std::ostream& stream, const bi& b) { for (int i = b.m_digits.size() - 1; i >= 0; i--) { if (i != b.m_digits.size() - 1) { stream << std::setfill('0') << std::setw(CLEN); } stream << b.m_digits[i]; } return stream; } bi operator+( const bi& lhs, const bi& rhs ) { bi result; size_t ls = lhs.m_digits.size(); size_t rs = rhs.m_digits.size(); size_t lm = std::min(ls, rs);
result.m_digits.reserve(std::max(ls, rs) + 1); int p = 0; for (size_t i = 0; i < lm; i++) { int c = lhs.m_digits[i] + rhs.m_digits[i] + p; result.m_digits.push_back(c % BASE); p = c / BASE; } if (ls > rs) { for (size_t i = lm; i < ls; i++) { int c = lhs.m_digits[i] + p;
result.m_digits.push_back(c % BASE); p = c / BASE; } } if (rs > ls) { for (size_t i = lm; i < rs; i++) { int c = rhs.m_digits[i] + p;
result.m_digits.push_back(c % BASE); p = c / BASE; } } if (p > 0) { result.m_digits.push_back(p); } return result; }
template T pp(int n) { if (n < 1 || n >= 10000000) return 0; if (n < 3) return 1; T a = 1; T b = 1; T c = 1; for (int i = 3; i < n; i++) { T t = a + b + c; c = b; b = a; a = t; } return a; }
template int test(int n) { auto t1 = std::chrono::steady_clock::now();
std::cout << pp(n) << std::endl;
auto t2 = std::chrono::steady_clock::now(); auto d_milli = std::chrono::duration_cast( t2 - t1 ).count(); std::cout << d_milli << std::endl; return d_milli; }
int main() { int n = 10000; int gmp = test(n); int my = test(n); std::cout << "gmp = " << gmp << " my = " << my << std::endl; std::cout << "my/gmp = " << my/gmp << std::endl; return 0; }