Как стать автором
Поиск
Написать публикацию
Обновить

Комментарии 29

Класс! Отличный пример практического использования расширения кольца!

Всего голосов 18: ↑18 и ↓0+18

А вы можете добавить на график метод с матрицами? Я всегда представлял, что есть иерархия методов вычисления чисел Фибоначчи:

итерация (для большинства случаев) -> бине (для быстрой оценки больших чисел) -> матрицы (для быстрого и точного вычисления чисел)

Кажется, что код с матрицами будет более естественным и понятным для большинства. И в нём гораздо сложнее сделать ошибку.

Но сам по себе метод, конечно, очень полезен, спасибо.

Всего голосов 5: ↑5 и ↓0+5


Код 
class Mat2x2:
    __slots__ = ['a', 'b', 'c', 'd']

    def __init__(self, a, b, c, d):
        self.a = a
        self.b = b
        self.c = c
        self.d = d

    def __matmul__(x, y):
        ans = x.copy()
        ans @= y
        return ans

    def __imatmul__(x, y):
        x.a, x.b, x.c, x.d = x.a * y.a + x.b * y.c, x.a * y.b + x.b * y.d, x.c * y.a + x.d * y.c, x.c * y.b + x.d * y.d
        return x

    def __pow__(self, exp):
        cur = Mat2x2(1, 0, 0, 1)
        base = self.copy()
        while exp:
            if exp & 1:
                exp -= 1
                cur @= base
            else:
                exp >>= 1
                base @= base
        return cur

    def copy(self):
        return Mat2x2(self.a, self.b, self.c, self.d)

    def __repr__(self):
        return f'{self.__class__.__name__}({self.a}, {self.b}, {self.c}, {self.d})'


class R5:
    def __init__(self, a=0, b=0):
        self.a = a
        self.b = b

    def __repr__(self):
        return f'R5({self.a}, {self.b})'

    def __str__(self):
        if self.a and self.b:
            return f'({self.a}{self.b:+}√5)'
        elif self.b:
            return f'{self.b}√5'
        else:
            return f'{self.a}'

    def __add__(x, y):
        return R5(x.a + y.a, x.b + y.b)

    def __sub__(x, y):
        return R5(x.a - y.a, x.b - y.b)

    def __mul__(x, y):
        return R5(x.a * y.a + x.b * y.b * 5, x.a * y.b + x.b * y.a)

    def __pow__(x, power):
        if power == 0:
            return R5(1, 0)
        elif power % 2 == 1:
            return x * (x ** (power - 1))
        else:
            sq = x ** (power // 2)
            return sq * sq

    def __floordiv__(x, n):
        return R5(x.a // n, x.b // n)


def fib_ite(n):
    c, p = 0, 1
    for _ in range(n):
        c, p = c + p, c
    return c


def fib_bine(n):
    return (R5(1, 1) ** n - R5(1, -1) ** n).b // 2 ** n


def fib_mat(n, *, fib_mat=Mat2x2(0, 1, 1, 1)):
    if n == 0:
        return 0
    fib_pow = fib_mat ** (n - 1)
    return fib_pow.d


# Проверка корректности
for i in range(0, 100):
    # print(fib_bine(i), fib_ite(i), fib_mat(i))
    assert fib_bine(i) == fib_ite(i) == fib_mat(i)


# Оценка времени работы
def time_fib(funcs, from_i=0, till_i=300000, step=20011, repeats=10):
    from timeit import timeit
    setup = 'from __main__ import ' + ', '.join(funcs)
    times = []
    for i in range(from_i, till_i, step):
        tts = [timeit(stmt=f'{func}({i})', setup=setup, number=repeats) / repeats for func in funcs]
        times.append([i] + tts)
        print(*times[-1])
    return times


# Отрисовка
def plot_fib(funcs, times):
    from matplotlib import pyplot as plt
    ax = plt.gca()
    base = [t[0] for t in times]
    for i, func in enumerate(funcs, start=1):
        cur = [t[i] for t in times]
        ax.plot(base, cur, label=func)
        for x, v in zip(base, cur):
            ax.annotate(f'{v:0.3}', xy=(x, v))
    ax.set_title('Время на возведение в степень')
    ax.set_xlabel('Номер числа Фибоначчи')
    ax.set_ylabel('Время на вычисление в секундах')
    ax.legend()
    plt.show()


funcs = ['fib_ite', 'fib_bine', 'fib_mat']
times = time_fib(funcs, from_i=0, till_i=300000, step=20011, repeats=10)
plot_fib(funcs, times)

Всего голосов 2: ↑2 и ↓0+2

Я, конечно, ожидал всякого, но разница аж в 4 раза! В итоге посмотрел на результаты промежуточных вычислений: длина чисел по формуле Бине получается примерно в 2.5 раза больше, чем в матричном виде, то есть уходит тупо больше операций на работу с длинными числами.

Комментарий пока не оценивали0

А в чём принципиальное отличие от матричного метода Кнута? И там, и там узким местом становится возведение элемента в степень.

В метода Кнута используется умножение матриц 2×2 — это 8 умножений, 4 сложения.

Вы же делаете 4 умножения, 2 сложения, 1 умножение на однозначное число. Но вам надо сделать это два раза (для A и для B): получаем 8 умножений, 4 сложения плюс 2 умножения на однозначное число.

Получается, что асимптотика та же, а количество операций даже чуть больше, чем для метода Кнута.

Так что присоединяюсь к предыдущему комментатору и реквестирую график с матрицами.

Всего голосов 22: ↑22 и ↓0+22

Пару лет назад была статья с интересным методом, вот он на питоне:

def fast_doubling(n):
    if n <= 2:
        return 1

    f_n = fast_doubling(n // 2)
    f_n_1 = fast_doubling(n // 2 + 1)
    if n % 2:
        return f_n ** 2 + f_n_1 ** 2
    else:
        return f_n * (2 * f_n_1 - f_n)

На 1000-м элементе работает намного быстрее, чем метод отсюда (предполагаю, из-за меньшего размера O, как и в той статье в сравнении с матричным методом), на бОльших числах не смог проверить, тк код из этой статьи выдал RecursionError. Интересно бы было добавить к сравнению.

Всего голосов 2: ↑2 и ↓0+2

Да, матричный метод "естественнее" для вычисления чисел Фибоначчи. Но суть заметки - в "реабилитации" формулы Бине.

Всего голосов 6: ↑6 и ↓0+6

Очень круто получилось, а связь между методами проследить не удаётся? Ну там, взаимно-однозначное соответствие пар с матрицами или что-нибудь в этом духе.

А то одинаковая суть алгоритма (быстрое возведение в степень) намекает.

Всего голосов 2: ↑2 и ↓0+2

Спасибо за реабилитацию. Всегда обидно за красивые математические формулы, которыми нельзя красиво воспользоваться на практике. Но в данном случае можно.

Единственное, что в реализации мозолит глаза, так это деление на 2^n. Его логично заменить на битовый сдвиг вправо, что явно быстрее.

Кроме того можно заметить что при каждом умножении числителей A и B возникает как минимум 1 множитель равный двум "2":

(1+sqrt(5))^2 = 2 * (3 + 1 * sqrt(5))

(1+sqrt(5))^3 = 4 * (4 + 2 * sqrt(5))

...

который можно сразу сократить с 2-кой из знаменателя. Т.е. эту 2-ку можно сокращать при каждом умножении с помощь того же двоичного сдвига (также реализуемо при быстром умножении). Тогда в конце деление на 2^n не потребуется. Такой подход как минимум сократит объём памяти для хранения больших чисел, хотя вероятно в ущерб производительности, т.к. сдвигать придётся каждой операции, а не в конце

Всего голосов 1: ↑1 и ↓0+1

Для быстрого вычисления (за логарифм от n) можно еще последовательности Люка использовать. Полезные формулы — вычисление F(n+m) и L(n+m) через F(n), F(m) и L(n), L(m). Ну еще F(2n), L(2n) и способ, аналогичный быстрому возведению в степень.

Всего голосов 2: ↑2 и ↓0+2
Вот такое получилось на базе последовательностей Люка:
def fib(n):
    b = 1 << (n.bit_length() - 1)
    f = 1
    l = 1
    a = 2
    while b != 1:
        b = b >> 1
        f = f * l
        l = l ** 2 + a
        if n & b != 0:
            a = 2
            fn = (f + l) >> 1
            l = fn + f + f
            f = fn
        else:
            a = -2
    return f
Всего голосов 2: ↑2 и ↓0+2

Люк рулит!
Всего голосов 2: ↑2 и ↓0+2

Спасибо, хороший стиль изложения. Ощутил себя на школьном факультативе :)

Всего голосов 3: ↑3 и ↓0+3

Есть один непонятный момент почему следующее верно?

Для того, чтобы значение оказалось целым (и корень из пяти сократился), нужно, чтобы числитель дроби формулы Бине представлял собой число вида: 0 + r*5^(1/2)

Почему при вычитании A - B первый компонент становится равен нулю?

Всего голосов 1: ↑1 и ↓0+1

Предположим, после всех вычислений числитель принял вид a+b*sqrt(5). Здесь a и b могут быть только целыми (поскольку все исходные данные целые, а мы складываем, умножаем и вычитаем). Но если a+b*sqrt(5) разделить на sqrt(5) (обычным образом), то получится a/sqrt(5)+b. Величина a/sqrt(5) не может быть целой ни при каком b, кроме нуля. Примерно так.

Всего голосов 1: ↑1 и ↓0+1

Деление можно представить как умножение на обратный элемент.


Обратный элемент:


(a + b \sqrt 5) ^ {-1}

домножим числитель и знаменатель на a - b \sqrt 5. Получим в знаменателе a^2 - 5b.


(a + b \sqrt 5) ^ {-1} = a/(a^2 - 5b) - b/(a^2 - 5b)\sqrt 5

Т.е. остаёмся в группе.

Всего голосов 1: ↑1 и ↓0+1

При умножении (a + \sqrt{5} b) (a - \sqrt{5} b) получается a^2 - 5 b^2. А у Вас в результате везде написано b^2 без квадрата.

Всего голосов 2: ↑2 и ↓0+2

Да, спасибо. Описка

Комментарий пока не оценивали0

Действительно, опечатался. Спасибо за замечание. (К сожалению, в комментарии уже исправить не получается.)

Комментарий пока не оценивали0

А существует ли способ работы с корнями (и другими иррациональными или бесконечно-дробными) числами без их конвертации в ЧсПТ или округления до ближайшего рационального?

Всего голосов 1: ↑0 и ↓1-1

Символьная манипуляция — способ работы с корнями без их конвертации. Данная статья как раз и является примером реализации простейшего символьного движка.

Вот пример работы WolframAlpha
Всего голосов 2: ↑2 и ↓0+2
Можно было классик с арифметикой замутить, чтобы дальше красиво было:

class R5:
    def __init__(self, a=0, b=0):
        self.a = a
        self.b = b

    def __repr__(self):
        return f'R5({self.a}, {self.b})'

    def __str__(self):
        if self.a and self.b:
            return f'({self.a}{self.b:+}√5)'
        elif self.b:
            return f'{self.b}√5'
        else:
            return f'{self.a}'

    def __add__(x, y):
        return R5(x.a + y.a, x.b + y.b)

    def __sub__(x, y):
        return R5(x.a - y.a, x.b - y.b)

    def __mul__(x, y):
        return R5(x.a * y.a + x.b * y.b * 5, x.a * y.b + x.b * y.a)

    def __pow__(x, power):
        if power == 0:
            return R5(1, 0)
        elif power % 2 == 1:
            return x * (x ** (power - 1))
        else:
            sq = x ** (power // 2)
            return sq * sq

    def __floordiv__(x, n):
        return R5(x.a // n, x.b // n)


def fib_bine(n):
    return (R5(1, 1) ** n - R5(1, -1) ** n).b // 2 ** n


def fib_ite(n):
    c, p = 0, 1
    for _ in range(n):
        c, p = c + p, c
    return c


for i in range(1, 100):
    assert fib_bine(i) == fib_ite(i)
Всего голосов 4: ↑4 и ↓0+4

Именно!

Комментарий пока не оценивали0

То, что полученное кольцо пар (a, b) является полем (т.е. можно делить такие целочисленные пары на ненулевые целочисленные пары и получать сновва целочисленные пары) неверно. Действительно, если попробовать разделить пару (1, 0) (т.е. просто 1) на пару (0, 1) (т.е. \sqrt{5}), то мы должны получить пару (0, 1/5) = (0, 0.2) (т.е. \sqrt{5}/5 = 1/\sqrt{5}), что не является целочисленной парой!

В чём же ошибка в рассуждении? Всё просто: если система линейных уравнений (СЛУ) с целочисленными коэффициентами имеет ненулевой определитель, то это значит, что оно имеет решение в поле, но не обязательно в кольце (например, простое СЛУ 2x = 1 имеет определитель 2, но и единственное решение 1/2 = 0.5).

В данном случае, чтобы получить поле, нужно рассматривать пары рациональных чисел (а не целых).

P.S. Это не влияет на суть статьи: мы не пользовались тем, что это поле, а все моменты вне этих абстракций тщательно описанны. Но ошибка есть ошибка.

Всего голосов 4: ↑4 и ↓0+4

Я и чувствовал, что поле все-таки не получится... Спасибо за очень важную поправку!

Комментарий пока не оценивали0

Но если бы a и b были рациональными, а не целыми, то поле получилось бы (о чем, собственно, Вы и написали).

Комментарий пока не оценивали0

Хотел бы еще отметить интересные св-ва представленного поля. Оно является алгебраическим расширением поля Q рациональных чисел (если мы считаем, что a, b \in Q), а именно Q(\sqrt5) - число \sqrt5 является корнем многочлена x^2-5 с рациональными коэффициентами.

Комментарий пока не оценивали0

Если идёт речь о прямом (в том или ином смысле) вычислении по формуле Бине, то можно заметить, что член \left((1 - \sqrt{5})/2\right)^n считать нет смысла: число \left((1 - \sqrt{5})/2\right)^nпо жизни находится между -0.62и -0.61, а значит его степени будут (по модулю) всё меньше и меньше (при возведении числа между 0 и 1 в степень n оно начинает катострофически быстро стремится к нулю при росте ((1+\sqrt{5})/2)^n/\sqrt{5}, а затем сделать один из следующих вариантов:
1. Округлить до ближайшего целого. Поскольку, уже с нулевого члена (n=0) модуль второго слагаемого ((1-\sqrt{5})/2)^n/\sqrt{5}будет по модулю меньше 1/2, то банальное округление до ближайшего целого числа сделает своё дело.
2. Определить знак ((1-\sqrt{5})/2)^n/\sqrt{5} и округлить в соответствующую сторону до ближайщего целого. Действительно, знак это члена – +когда n делится на 2, и -, когда не делится. Следовательно, если мы знаем, что нужно добавить число с данным знаком и по модулю <1и получить целое число, то это означает округление в сторону данного знака до ближайшего целого.

Главный вопрос. Пусть мы посчитали пару (a, b) = a + \sqrt{5} b. Как же её поделить на \sqrt{5} и округлить?

На него можно ответить так. Ну, \sqrt{5} b / \sqrt{5} = b, которое можно просто вычесть и времено про него забыть (дабы на округление оно не влияет). Теперь же нам нужно найти такое число c, что a/\sqrt{5} \approx c. Тут два варианта: либо по-"глупому" написать ряд для 1/\sqrt{5} (благо, он классно сходится) или посчитать префикс цепной дроби того числа (или ещё что-нибудь) и приближать ими (но это нудно и нужно писать оценки), либо бинпоиском искать такое число c, чьё c^2 будет ближе всех к 5a^2.

С другой стороны есть другой формальный подход. Кольцо таких пар (формально, \mathbb{Z}[\sqrt{5}]) обладает ещё одной интересной операцией: операцией сопряжения, которая паре (числу) (a, b) сопоставит пару \overline{(a,b)} := (a, -b). Выглядит как странная операция, которая меняет второй аргумент, но не всё так просто. Во-первых она сохраняет операции кольца (и поля, если рассматривать пары рациональных чисел):

  1. \overline{(a, b) + (c, d)} = \overline{(a, b)} + \overline{(c, d)}. Действительно,

\overline{(a, b) + (c, d)} = \overline{(a + b, c + d)} = (a + b, -c - d) = (a, -b) + (c, -d) = \overline{(a, b)} + \overline{(c, d)}

  1. \overline{(a, b) \cdot (c, d)} = \overline{(a, b)} \cdot \overline{(c, d)}. Действительно,

\begin{multline}\overline{(a, b) \cdot (c, d)} = \overline{(ab - 5cd, ac + bd)} = (ab - 5cd, -ac - bd)\\ = (a, -b) \cdot (c, -d) = \overline{(a, b)} \cdot \overline{(c, d)}\end{multline}

А во-вторых, он позволяет упростить много счёта. Давайте немного обозначим: \alpha := 1 + \sqrt{5}, \beta := 1 - \sqrt{5}. Тогда нам надо посчитать (\alpha^n + \beta^n)/2^n/\sqrt{5}. Но заметим, что \overline{\beta} = \alpha, а тогда \beta^n = \overline{\alpha}^n = \overline{\alpha^n}, т.е. высчитав пару для \alpha^n мы можем буквально мгновенно получить пару для \beta^n: просто измените знак второго члена на противоположный. И счёт фактически сократится в два раза.

Всего голосов 4: ↑4 и ↓0+4
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Истории

Ваш аккаунт

Разделы

Информация

Услуги

Техническая поддержка
© 2006–2024, Habr