N-е число Фибоначчи за O(log N)

Это абсурд. Жонглирование размытыми категориями типа «практически всегда» не делает O(log n) из O(n).

еще чуть хуже :)

Итак, нам нужно найти бесконечную последовательность, удовлетворяющую соотношению
x_{n+k} + a_{k} * x_{n+k} +… + a_0 * x_{n} = 0.
Заметим, что если {x_n} и {y_n} — две подходящих последовательности, то последовательность {x_n + y_n} и {a*x_n} — тоже подходящие, т.е. интересующие нас последовательности образуют линейное подпространство в пространстве всех последовательностей. Его размерность равна k+1, так как последовательность, с одной стороны, однозначно задается k числами {x_0, ..., x_k}, с другой — любым двум таким различным наборам чисел соответствуют различные последовательности.

Рассмотрим многочлен:
f(t) = t^{k+1} + a_k * t^{k} +… + a_0 = 0

Как известно из основной теоремы алгебры, у него ровно N корней с учетом кратности.
Пусть c — корень кратности s+1. Тогда для любых констант C0, ..., CS последовательность
x_k = (C0 + C1*k + C2*k^2 +… + CS*k^s) * c^k будет решением (т.е. произвольный многочлен P(k) степени не выше S, умноженный на c^k).
Докажем это, т.е. что
P(n+k+1)*c^{n+k+1} + a_k * P(n+k) * c^{n+k} +… + a_0 * P(n) * c^n = 0.
В силу линейности, достаточно рассматривать многочлены вида P(k) = k^i. Перейдем к такому случаю, и сократим на c^n:
(n+k+1)^i * c^{k+1} + a_k * (n+k)^i * c^{k} +… + a_0 * n^i = 0.
Разложим это по степеням n (и вынесем из каждой скобки c в максимально возможной степени):
n^i * (c^{k+1} + a_k * c^k +… + a_0) + i*n^{i-1} * c * ((k+1) * c^{k} + k * c^{k-1} +… + a_1) +… = 0
Из общей теории получаем, что, так как c — корень кратности s+1 многочлена f, то f© = f'© = f''© =… = f^(s) © = 0 (т.е. все производные, вплоть до s-й, равны нулю).
Запишем явно производный от f вплоть до s-й:
f^(0) = t^{k+1} + a_k * t^{k} +… + a_0
f^(1) = (k+1)*t^k + a_k * k * t^{k-1} +…
…

Замечаем, что эти производные в точности равны тому, что стоит в скобках в разложении, если в них вместо t подставить c. А так как c — это корень S-й степени, то все производные равны нулю, т.е. все скобки равны нулю, и вся сумма равна нулю, что и требовалось доказать.

Таким образом, каждый корень кратности S дает нам S независимых решений. Решения, соответствующие разным корням, тоже независимы, сумма всех кратностей равна N = > мы нашли все решения.

Для тех, кому лень было читать вывод — готовый алгоритм:
1. Ищем последовательность, удовлетворяющую уравнению x_{n+k} + a_{k} * x_{n+k} +… + a_0 * x_{n} = 0
2. Пишем характеристический многочлен: t^{k+1} + a_k * t^{k} +… + a_0
3. Находим корни характеристического многочлена, т.е. представляем его в виде (t-c1)^s1 * (t-c2)^s2 * ..., где все ci попарно различны.
4. Общий вид последовательность: x_n = (C1_0 + C1_1*n +… + C1_s1 * n^{s1-1})*c1^n + (C2_0 + C2_1*n +… + C2_s2 * n^{s2-1})*c2^n +…
5. Если у нас фиксировано начало — то подставляем в эту формулу n=0,1,...,k и находим все константы.

Например, для чисел Фибоначчи:
x_{n+2} — x_{n+1} — x_{n} = 0

t^2 — t — 1 = 0
D = 5
t = (1 +- sqrt(5)) / 2
x_n = A*((1+sqrt(5))/2)^n + B*((1-sqrt(5))/2)^n

Разным значениям A и B соответствуют разные значения первых двух членов.
x_0 = A+B
x_1 = A(1+sqrt(5))/2 + B(1-sqrt(5))/2
В стандартном случае x_0 = 0, x_1 = 1
A+B = 0
A*sqrt(5)/2 — B*sqrt(5)/2 = 1

A+B=0
A-B=2/sqrt(5)

A = 1/sqrt(5)
B = -1/sqrt(5)

(что соответствует формулам из вики:D)

В ветку призывается nsinreal, которому не лень переводить TeX-овские формулы в картинки.

Вообще, формулы в текстовом виде воспринимать ой как сложно.

Что же говорить о тех людях, у которых линейка читается один семестр и то для галочки и зачета. Кстати, а как у вас было, линейна алгебра и аналитическая геометрия — это один предмет?

Меня все-таки интересует полноценное объяснение с возможностью применения для всех линейно-рекуррентных функций. Ну т.е. чтобы был хотя-бы намек как такие формулы генерировать. Или проще тупо через матрицы считать?

Кстати, во что она откомпилируется на x86?

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

def pow(a,b):
  res = 1
  starta = a
  while b:
     if b % 2 == 1 :
        res *= starta
     a *= a
     b /= 2

    res *= a
    return res

У Вас неудачная реализация рекурсивного варианта.
Вот так:

Mat pow2(const Mat &a, int b) {
    if (b == 0) {
        return Mat();
    }
    Mat res = pow2(a, b / 2);
    if (b % 2) {
        res *= a;
    }
    return res;
}

оно работает за 1.42 против 1.7 с O3 и 1.26 против 3.27 с O2. (3 млн. шагов, матрицы 3х3)

Проще можно определить как вполне объективным понятием, если учитывать, что людей много и большинству из них проще что-то.