Формула Бине без плавающей точки / forpes.ru

Главная
Формула Бине без плавающей точки

Формула Бине без плавающей точки +71

16.01.2022 11:15

catstail1954 27 8100 Источник

Красивая формула Бине для чисел Фибоначчи содержит иррациональность - квадратный корень из пяти. Это делает ее непригодной для точного вычисления больших чисел Фибоначчи. Это кажется вполне очевидным. Предлагаю способ, как избавиться от зловредного корня и сделать формулу Бине пригодной для точных вычислений.

Интересно? Читайте дальше

Как хорошо известно, числа Фибоначчи – это целочисленная последовательность, первые два члена которой равны единице, а каждый последующий равен сумме двух предыдущих. За 500 лет, прошедших с момента ввода этой последовательности в математический обиход, она основательно изучена. Открыто много интереснейших формул с участием чисел Фибоначчи. Но одной из “непреходящих” учебных задач является вычисление чисел Фибоначчи. Для этого придумано много способов: от прямой рекурсии, основанной на формуле:

${F}_{n}={F}_{n-1}+{F}_{n-2}$

до матричного метода, описанного, например, в книге Д.Кнута [1]. Большая часть этих подходов (кроме матричного метода Кнута) основаны на рекуррентных свойствах последовательности Фибоначчи и позволяют вычислить величину Fn в лучшем случае за время O(n). Матричный метод Кнута (использующий матричное возведение в степень) позволяет вычислить число Фибоначчи за логарифмическое время [2].

Особняком в этом ряду алгоритмов располагается формула Бине (известная еще Муавру) имеющая вид:

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

Сказанное означает, что вычисления не будут точными; в них вносится погрешность ограничения. Мне однажды попалась на глаза публикация [3] в которой использовалась формула Бине для вычисления очень большого числа Фибоначчи, но реализация предполагала использование плавающей арифметики сверхвысокой разрядности (с тем, чтобы нужное число полностью уместилось в мантиссу).

Мы пойдем совсем другим путем!

Для начала, рассмотрим множество чисел вида:

при целых a и b. Легко убедиться в том, что это множество алгебраически замкнуто относительно операций обычного сложения и умножения:

Более того, умножение и сложение будут коммутативными, что тоже легко проверяется непосредственно. Кроме того, ноль и единица хорошо "вписываются" в рассматриваемое множество:

Вполне натурально реализуется и вычитание подобных чисел:

Можно определить и деление (разумеется, лишь в случае, когда делитель отличен от нуля). Результат деления можно определить как корень уравнения:

$(a+b\sqrt{5})x=(c+d\sqrt{5})$

Пусть

$x=e+f\sqrt{5}$

Тогда предыдущее равество эквивалентно следующему:

$(a+b\sqrt{5})(e+f\sqrt{5})=(c+d\sqrt{5})$

Раскрывая произведение в левой части последнего равенства, получим систему линейных уравнений для определения неизвестных e и f:

$(a+b\sqrt{5})(e+f\sqrt{5})=(ae+5bf)+(af+be)\sqrt{5}$

Отсюда:

$\left.\begin{matrix} ae+5bf & = & c\\ be+af & = & d \end{matrix}\right\}$

Главный определитель этой системы равен:

$\begin{vmatrix} a & 5b\\ b & a \end{vmatrix} = {a}^{2}-5{b}^{2}$

Поскольку a и b здесь целые числа, то значение определителя всегда отлично от нуля, а значит, система имеет единственное решение и деление определяется корректно. Впрочем, мы увлеклись. Деление нам не понадобится.

Мы пришли к тому, что рассматриваемое множество с операциям сложения и умножения образует кольцо [4].

А теперь - самое главное! Зачем нам нужен корень из пяти? Никто не мешает реализовать арифметику на множестве пар (a,b), в которой сложение, вычитание и умножение будет описываться формулами:

Таким образом, можно “благополучно забыть” про корень из пяти и реализовать прямое вычисление по формуле Бине. Для того, чтобы значение оказалось целым (и корень из пяти сократился), нужно, чтобы числитель дроби формулы Бине представлял собой число вида:

$0+r\sqrt{5}$

которое мы отождествим с "обычным" числом

$r\sqrt{5}$

Деление этого "обычного" иррационального числа на корень из пяти и даст нам искомый целый результат. Естественно, что в действительности делить не требуется, достаточно вычислить (используя описанную выше арифметику пар) два бинома:

а потом произвести вычитание. Не представляет труда реализовать этот подход на любом языке программирования. Мы сделаем это на Питоне (привлекает неограниченная разрядность целых в этом языке).

def prod_pairs(a,b):
    return (a[0]*b[0]+5*a[1]*b[1],a[0]*b[1]+a[1]*b[0])
def sub_pairs(a,b):
    return (a[0]-b[0],a[1]-b[1])
def pow_pair(a,n):
    c=a
    for _ in range(n-1):
        c=prod_pairs(c,a)
    return c
def fib_bine(n):
    x1=pow_pair((1,1),n)
    x2=pow_pair((1,-1),n)
    z=sub_pairs(x1,x2)
    return z[1]//(2**n)

Комментарии излишни - все очень просто. Сразу же возникает вопрос, а можно ли ускорить этот код? Очевидно, что "узким местом" здесь являтся возведение пары в целую степень. Для ускорения этой операции имеется стандартный прием - быстрое возведение в степень (этим же приемом воспользовался и автор [2]). Идея ускорения состоит в том, что для вычисления xⁿ вычисляется цепочка x -> x² -> x⁴ ->…->x^2k до тех пор, пока 2^k<=n, а затем аналогичным образом вычисляется x^(n-2k).

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

def pow_pair(a,n):
if (n==1):
    return a
c=copy(a)
k=1
while k*2<=n:
    if k<=n:
       c=prod_pairs(c,c)
       k=k*2
    p=n-k
    if p>=1:
       tmp=pow_pair(a,p)
       return prod_pairs(tmp,c)
    else:
       return c

Использование этого приема позволяет вычислять числа Фибоначчи за время, близкое к логарифмическому по формуле Бине и без использования арифметики с плавающей точкой. Для сравнения производительности предлагаемого метода с методом, основанным на простой рекурсии, написан следующий простейший код:

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

И что же? Несмотря на очевидную простоту кода fib_ite, функция fib_bine показывает значительно лучшие результаты. Так, на компьютере автора четырехсоттысячное число Фибоначчи по описываемому алгоритму вычисляется примерно за 2 сек, а прямыми итерациями – за 27 сек. На прилагаемом рисунке показаны разультаты тестов:

По горизонтальной оси откладывается номер рассчитываемого числа Фибоначчи, по вертикальной – время в секундах.

Получается, что формула Бине вполне пригодна для быстрых и точных вычислений чисел Фибоначчи.

Спасибо, что дочитали до конца, а также искреннее спасибо авторам, на которые я ссылался в этой заметке:

1. Д.Кнут Искусство программирования на ЭВМ, т.1, Основные алгоритмы. – М: Вильямс, - 2017. - 720 C.

2. N-е число Фибоначчи за O(log N) https://habr.com/ru/post/148336/

3. Расчет миллионного числа Фибоначчи https://habr.com/ru/company/skillfactory/blog/555914/

4. С.Ленг Алгебра. М.: Наука, - 1965. - 431 C.

Комментарии (27)

samsergey
16.01.2022 11:26
#23946109
+18
Класс! Отличный пример практического использования расширения кольца!

h1pp0

16.01.2022 11:37

#23946147

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

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

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

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

ShashkovS

17.01.2022 12:32

#23950287

Код

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)

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

DistortNeo
16.01.2022 11:57
#23946199
+21
А в чём принципиальное отличие от матричного метода Кнута? И там, и там узким местом становится возведение элемента в степень.

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

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

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

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

tetelevm
16.01.2022 13:16
#23946495
+2
Пару лет назад была статья с интересным методом, вот он на питоне:
```
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. Интересно бы было добавить к сравнению.

catstail1954 Автор
16.01.2022 14:03
#23946649
+4
Да, матричный метод "естественнее" для вычисления чисел Фибоначчи. Но суть заметки - в "реабилитации" формулы Бине.
1. aamonster
  16.01.2022 17:18
  #23947359
  +2
  Очень круто получилось, а связь между методами проследить не удаётся? Ну там, взаимно-однозначное соответствие пар с матрицами или что-нибудь в этом духе.
  
  А то одинаковая суть алгоритма (быстрое возведение в степень) намекает.

Kohelet
16.01.2022 15:12
#23946907
+2
Для быстрого вычисления (за логарифм от n) можно еще последовательности Люка использовать. Полезные формулы — вычисление F(n+m) и L(n+m) через F(n), F(m) и L(n), L(m). Ну еще F(2n), L(2n) и способ, аналогичный быстрому возведению в степень.
1. stan_volodarsky
  16.01.2022 21:11
  #23948031
  Тут есть разбор методов, в том числе и последовательность Люка:
  https://stackoverflow.com/questions/14661633/finding-out-nth-fibonacci-number-for-very-large-n
1. paluke
  17.01.2022 12:17
  #23950191
  +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
  1. ShashkovS
    18.01.2022 11:55
    #23954939
    +2
    
    Люк рулит!

sci_nov
16.01.2022 20:15
#23947907
+2
Спасибо, хороший стиль изложения. Ощутил себя на школьном факультативе :)

ARad
17.01.2022 08:50
#23949261
+1
Есть один непонятный момент почему следующее верно?

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

Почему при вычитании A - B первый компонент становится равен нулю?
1. catstail1954 Автор
  17.01.2022 08:57
  #23949275
  Предположим, после всех вычислений числитель принял вид a+b*sqrt(5). Здесь a и b могут быть только целыми (поскольку все исходные данные целые, а мы складываем, умножаем и вычитаем). Но если a+b*sqrt(5) разделить на sqrt(5) (обычным образом), то получится a/sqrt(5)+b. Величина a/sqrt(5) не может быть целой ни при каком b, кроме нуля. Примерно так.

primetalk
17.01.2022 10:52
#23949711
+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. lounres
  17.01.2022 17:16
  #23951763
  +2
  При умножении (a + \sqrt{5} b) (a - \sqrt{5} b) получается a^2 - 5 b^2. А у Вас в результате везде написано b^2 без квадрата.
  1. catstail1954 Автор
    17.01.2022 17:17
    #23951765
    Да, спасибо. Описка
  1. primetalk
    17.01.2022 17:19
    #23951775
    Действительно, опечатался. Спасибо за замечание. (К сожалению, в комментарии уже исправить не получается.)

StjarnornasFred
17.01.2022 10:57
#23949747
-1
А существует ли способ работы с корнями (и другими иррациональными или бесконечно-дробными) числами без их конвертации в ЧсПТ или округления до ближайшего рационального?
1. Soukhinov
  17.01.2022 12:19
  #23950209
  +2
  Символьная манипуляция — способ работы с корнями без их конвертации. Данная статья как раз и является примером реализации простейшего символьного движка.
  
  Вот пример работы WolframAlpha

ShashkovS

17.01.2022 11:51

#23950057

Можно было классик с арифметикой замутить, чтобы дальше красиво было:

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)

catstail1954 Автор
17.01.2022 17:17
#23951771
Именно!

lounres
17.01.2022 17:32
#23951861
+3
То, что полученное кольцо пар (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. Это не влияет на суть статьи: мы не пользовались тем, что это поле, а все моменты вне этих абстракций тщательно описанны. Но ошибка есть ошибка.
1. catstail1954 Автор
  18.01.2022 11:16
  #23954635
  Я и чувствовал, что поле все-таки не получится... Спасибо за очень важную поправку!
1. catstail1954 Автор
  18.01.2022 11:21
  #23954661
  Но если бы a и b были рациональными, а не целыми, то поле получилось бы (о чем, собственно, Вы и написали).

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

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

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

С другой стороны есть другой формальный подход. Кольцо таких пар (формально, $\mathbb{Z}[\sqrt{5}]$ ) обладает ещё одной интересной операцией: операцией сопряжения, которая паре (числу) сопоставит пару $\overline{(a,b)} := (a, -b)$ . Выглядит как странная операция, которая меняет второй аргумент, но не всё так просто. Во-первых она сохраняет операции кольца (и поля, если рассматривать пары рациональных чисел):
1. $\overline{(a, b) + (c, d)} = \overline{(a, b)} + \overline{(c, d)}$ . Действительно,
1. $\overline{(a, b) \cdot (c, d)} = \overline{(a, b)} \cdot \overline{(c, d)}$ . Действительно,
А во-вторых, он позволяет упростить много счёта. Давайте немного обозначим: $\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$ : просто измените знак второго члена на противоположный. И счёт фактически сократится в два раза.