N-e число обобщённых Фибоначчи за O(log N)

[vesper-bot]

А не проще найти N-е обычное число Фибоначчи и потом умножить его на С? Просто потому, что если ряд Фибоначчи порождается парой (0, С) это самое С можно вынести за скобки. А дальше взять обычный алгоритм из тех, что находят N-е число Фибоначчи за O(log N) и умножать результат.

За сам алгоритм, который тут есть, спасибо — мне пока не попадался.

[nklvv]

Можно и так, но этот способ мне показался более изящным.
Спасибо за отзыв!

[iboltaev]

двоичное возведение в степень делается проще, без разложения на степени двойки (точнее, разложение присутствует неявно):

def binpow(m: Matrix, pow: Int): Matrix = {
  assert(pow >= 0, "pow must be non-negative")
  if (pow == 0) Matrix.identity
  else if (pow == 1) m
  else {
    val p = binpow(m, pow / 2)
    if (pow % 2 == 0) p * p
    else m * p * p 
  }
}

[mayorovp]

Ох, а рекурсия-то тут нахрена?

Есть же нормальный алгоритм быстрого возведения в степень чего угодно:

1. пусть a возводимый в степень объект, b — единица относительно умножения, а n — показатель степени;
2. если n чётно, то переприсвоим (n ← n/2, a ← a*a)
3. если n нечётно, то переприсвоим (n ← n-1, b ← b*a)
4. повторяем шаги 2-3 пока n > 0
5. теперь результат в переменной b

[nklvv]

Действительно, ваш алгоритм эффективнее на больших числах:

Спасибо, учту

[vlanko]

А если возвести матрицу в степень через NumPy?

[daniel-daniels]

Так полагаю, речь в статье об асимптотической сложности в первую очередь, а не в разнице реализаций на чистом Питоне и на С (на всякий случай поясню, что под реализацией на С имею ввиду реализацию, собственно, NumPy \кеп).

[nklvv]

Верно, готовые библиотеки использовать в рамках статьи я не планировал

[AlbertEinsteinEpoch21]

Все гораздо проще как я знаю, используя то что для фибоначи:
F(2k) = F(k)[2F(k + 1) − F(k)]
F(2k+1) = F(k+1)^2 + F(k)^2
www.nayuki.io/page/fast-fibonacci-algorithms

def fibonacci(n):
if n < 0:
raise ValueError(«Negative arguments not implemented»)
return _fib(n)[0]

def _fib(n):
if n == 0:
return (0, 1)
else:
a, b = _fib(n // 2)
c = a * (b * 2 — a)
d = a * a + b * b
if n % 2 == 0:
return (c, d)
else:
return (d, c + d)

[nklvv]

Ваш алгоритм вычисляет N-е число Фибоначчи из классического ряда Фибоначчи. Если же заменить

return (0, 1)

на

return (0, C)

(где C — некая константа), на выходе получается некорректный результат.
По бенчмаркам ваш код действительно эффективен на больших числах.

[AlbertEinsteinEpoch21]

В приведённом в статье примере есть смещение на 1)
Если исправить так, то будет корректно.
def Fib(num, key):
fib_1 = 0
fib_2 = key
for dec in range(num-1):
fib_1, fib_2 = fib_2, fib_1+fib_2
return fib_2

Что бы работало в Вашей логике(начальные члены 0, key)
можно сделать на пример так:
def fib2(n,key=1):
r=[]
while n>1:
if n%2==0:r+=[0];n=n//2
else:r+=[1];n=(n-1)//2
a=1;b=1
for i in r[::-1]:
a,b=a * (b * 2 — a),a * a + b * b
if i:a,b=b,a+b
if key>1:return a*key
return a
При этом мы избавляемся от рекурсии что часто полезно.
Хотя здесь разницы почти нет, при n=1кк выигрыш в пользу рекурсии но не большой.
Мне нравятся матрицы) но в рамках статьи польза от их использования далеко не очевидна, но умение говорить на этом языке (языке матриц, и оптимизации операций с ними) не повредит.
Спасибо за статью, и комментарии.

[AlbertEinsteinEpoch21]

Для интереса уменьшил число умножений чисел Фибоначчи за счёт тождества Кассини,
таким образом на каждом шаге делаем всего 2 умножения)
def fib2(n,key=1):
r=[]
while n>0:
if n%2==0:n=n//2;c=0
else:n=(n-1)//2 ;c=1
r+=[[c,n%2]];
a=0;b=1
for i,j in r[::-1]:
c=a*b;d=a*a;
e=c+d;
if j:e+=-1
else:e+=1
a,b=c+c-d,d+e
if i:a,b=b,a+b
if key>1:return a*key
return a

[ShashkovS]

В питоне можно переопределять естественные операции умножения и возведения в степень, чтобы не писать MultiplyMatrix(a, b), а писать коротко a @ b

Для этого нужно определить __matmul__ и __pow__.

А ещё для вычислений с боооольшими числами помогает библиотека gpmy2.

Чтобы считать время выполнения, можно использовать библиотеку timeit, и запускать код десяток раз. Тогда графики будут гладкими, а не такими скачущими.

Код

from timeit import timeit
from gmpy2 import mpz
from matplotlib import pyplot as plt


class Mat2x2:
    __slots__ = ['a', 'b', 'c', 'd', 'use_mpz']

    def __init__(self, a, b, c, d, use_mpz=False):
        self.a = mpz(a) if use_mpz else a
        self.b = mpz(b) if use_mpz else b
        self.c = mpz(c) if use_mpz else c
        self.d = mpz(d) if use_mpz else d
        self.use_mpz = use_mpz

    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, self.use_mpz)
        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, self.use_mpz)

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


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


def matrix_fib_mpz(n, *, fib_mat=Mat2x2(0, 1, 1, 1, use_mpz=True)):
    fib_pow = fib_mat ** n
    return fib_pow.d


def fib(num, key=1):
    fib_1 = 0
    fib_2 = key
    for dec in range(num):
        fib_1, fib_2 = fib_2, fib_1 + fib_2
    return fib_2


setup = 'from __main__ import fib, matrix_fib, matrix_fib_mpz'
number = 10
from_i = 0
till_i = 500000
step = 51237
fib_base = []
run_times_fib = []
run_times_matrix_fib = []
run_times_matrix_fib_mpz = []
for i in range(from_i, till_i, step):
    assert fib(i) == matrix_fib(i)
    t1 = timeit(stmt=f'fib({i})', setup=setup, number=number) / number
    t2 = timeit(stmt=f'matrix_fib({i})', setup=setup, number=number) / number
    t3 = timeit(stmt=f'matrix_fib_mpz({i})', setup=setup, number=number) / number
    fib_base.append(i)
    run_times_fib.append(t1)
    run_times_matrix_fib.append(t2)
    run_times_matrix_fib_mpz.append(t3)
    print(i, t1, t2, t3)

ax = plt.gca()
ax.plot(fib_base, run_times_fib, label='simple fib')
ax.plot(fib_base, run_times_matrix_fib, label='matrix fib')
ax.plot(fib_base, run_times_matrix_fib_mpz, label='matrix fib with mpz')
for x, v1, v2, v3 in zip(fib_base, run_times_fib, run_times_matrix_fib, run_times_matrix_fib_mpz):
    ax.annotate(f'{v1:0.3}', xy=(x, v1))
    ax.annotate(f'{v2:0.3}', xy=(x, v2))
    ax.annotate(f'{v3:0.3}', xy=(x, v3))
ax.set_title('Время на возведение в степень')
ax.set_xlabel('Номер числа Фибоначчи')
ax.set_ylabel('Время на вычисление в секундах')
ax.legend()
plt.show()

[nklvv]

Да, с переопределением привычных операторов выглядит опрятнее.
Библиотеки numpy, gpmy2 и т.д. принципиально не хотел применять, хотелось реализовать на чистом Python.
За инфу про timeit спасибо, действительно бенчмарки выходят более гладкими

[Pavgran]

А ещё можно возводить не матрицу в степень, а характерестический многочлен. Например, матрица Q = [[0,1], [1,1]] имеет характерестический многочлен x^2-x-1, и, значит, Q^2-Q-E = 0.
Пользуясь этим, возводим x в нужную степень по модулю x^2-x-1, получаем многочлен вида ax+b, подставляем туда Q и получаем ответ. Перемножение двух матриц 2*2 — это 8 умножений и 6 сложений, перемножение двух двучленов по модулю x^2-x-1 — это 4 умножения и 3 сложения.
Формула умножения: (ax+b)(cx+d)=acx^2+(ad+bc)x+bd=ac(x+1)+(ad+bc)x+bd=(ac+ad+bc)x+(bd+ac).
В случае большего количества членов в рекуррентном соотношении находить степень матрицы с помощью характеристического многочлена становится ещё выгоднее.

[Pavgran]

В добавок к этому, если умножение медленнее, чем сложение (при увеличении размера чисел, например) полиномы можно перемножать с помощью техник быстрого умножения. (ax+b)(cx+d)=((a+b)(c+d)-bd)x+(ac+bd). Здесь 3 умножения и 4 сложения/вычитания.

Также для совсем больших чисел есть вариант воспользоваться китайской теоремой об остатках. Умножение при этом становится той же сложности, что и сложение, но идут потери на операции взятия остатков. Я не знаю, будет ли такой вариант быстрее. Его уже нужно аккуратно реализовывать и сравнивать.