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

#cpp #математика

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

Вопрос: как получить формулу для вычисления n-ого члена последовательности без нахождения
всех предыдущих членов? Очевидно, что рекурсивный подход для n=10000 не подходит. 
UPD: время выполнения программы для любых n не должно превышать 1 секунду. 
    

Ответы

Ответ 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;
}


Ответ 2

Смотрите, давайте составим аналог матричной формулы, которая давала асимптотику O(log
n) для стандартного Фибоначчи.

Итак, что у нас есть? a[i+1] = a[i] + a[i-1] + a[i-2]. Запишем это в матричном виде.

| ? ? ? |   | a[i]   |   | a[i+1] |
| ? ? ? | * | a[i-1] | = | a[i]   |
| ? ? ? |   | a[i-2] |   | a[i-1] |


(За неимением скобок матрицу заключил в вертикальные линии) Знаками вопроса обозначены
не определённые пока коэффициенты.

С последней строкой матрицы просто:

            | a[i]   |
| ? ? ? | * | a[i-1] | = a[i-1]
            | a[i-2] |


поэтому берём в качестве этой строки | 0 1 0 |. Точно так же средняя строка будет
просто | 1 0 0 |. С верхней строкой немного сложнее. Но если мы вспомним, что a[i+1]
= a[i] + a[i-1] + a[i-2], и запишем

            | a[i]   |
| ? ? ? | * | a[i-1] | = a[i] + a[i-1] + a[i-2]
            | a[i-2] |


сразу получим строку | 1 1 1 |.

Дальше немного математики. Умножение вектора

| a[i]   |
| a[i-1] |
| a[i-2] |


на нашу матрицу

| 1 1 1 |
| 1 0 0 |
| 0 1 0 |


сдвигает его на один индекс вверх. А умножение на эту же матрицу N раз сдвинет на
N индексов. Но умножить на матрицу N раз — это всё равно, что умножить один раз на
N-ую степень матрицы. А возводить в степень за логарифмическое время мы умеем — для
этого есть классический «русский крестьянский метод».

Реализуем. Пишу на C#, в нём есть встроенный целый тип неограниченной точности.

Для начала, вспомогательный класс, представляющий матрицу:

class Matrix3x3
{
    BigInteger[,] payload = new BigInteger[3,3];
    public BigInteger this[int i, int j]
    {
        get { return payload[i, j]; }
        set { payload[i, j] = value; }
    }

    public Matrix3x3 Multiply(Matrix3x3 other)
    {
        var result = new Matrix3x3();
        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                for (int k = 0; k < 3; k++)
                    result.payload[i, j] += payload[i, k] * other.payload[k, j];
        return result;
    }

    static public Matrix3x3 Unit()
    {
        var result = new Matrix3x3();
        for (int k = 0; k < 3; k++)
            result.payload[k, k] = 1;
        return result;
    }
}


Стандартнее не бывает.

Можно было бы ещё определить класс для вектора, сделайте это сами для красоты.

Итак, алгоритм:

// заполняем матрицу
var matrix = new Matrix3x3();
matrix[0, 0] = 1;
matrix[0, 1] = 1;
matrix[0, 2] = 1;
matrix[1, 0] = 1;
matrix[2, 1] = 1;

// возводим в степень
var matrixPower = Matrix3x3.Unit();
for (var remainingExponent = N; remainingExponent > 0; remainingExponent /= 2)
{
    if (remainingExponent % 2 != 0)
        matrixPower = matrixPower.Multiply(matrix);
    matrix = matrix.Multiply(matrix);
}

// умножаем вектор (1 1 1) на матрицу
// нам нужно в общем-то только первое число результата, так что остальные не считаем
BigInteger result = 0;
for (int j = 0; j < 3; j++)
    result += matrixPower[0, j] * 1;


Всё!



Кстати, это общий способ решения рекуррентных соотношений. (Ну, один из них.)


Ответ 3

С учетом вашего 

return f(n-1)+f(n-2)+f(n-3); 


хуже написать просто невозможно... У вас будет экспоненциальный рост, так что уже
при n=50 вы рискуете не дожить до получения результата :)

Например, для f(10) значение f(3) вычисляется 44 раза (вместо одного!), для f(20)
- 19513 раз, для f(30) - 8646064 раза. Именно вычисляется, как 3 вызова функции...

В лоб - куда короче:

unsigned long long F(int N)
{
    if (N < 3) return 1;

    unsigned long long f0 = 1, f1 = 1, f2 = 1;
    for(int i = 3; i <= N; ++i)
    {
        unsigned long long f = f0+f1+f2;
        f0 = f1;
        f1 = f2;
        f2 = f;
    }
    return f2;
}


Проблема только с переполнением: F(10000) - 2647-значное число...
Время счета этого числа "методом в лоб" (я использую MAPM) - порядка 0.1с.


Ответ 4

С точки зрения математики, явные формулы из рекуррентных получаются с помощью аппарата
производящих функций. 

Называются такие числа числами трибоначчи. Полный вывод слишком длинный, но явная
формула тут:

https://ru.m.wikipedia.org/wiki/Числа_трибоначчи

Прошу прощения за ссылку, с телефона сложно картинки вырезать. Исправлю чуть позже.


Ответ 5

Если известно максимальное N, достаточно взять и один раз все протабулировать.

А от переполнения поможет библиотека GMP или ее аналог.