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


Думаю что никого не удивит “дефолтное” решение для вычисления N-го числа Фибоначчи:

static int fib(int n)
{
  return n > 1 ? fib(n - 1) + fib(n - 2) : n;
}


Также, есть достаточно “модная” форма записи этого решения которая использует Y-комбинатор:

Func<int, int> fib  = Y<int,int>(f => n => n > 1 ? f(n - 1) + f(n - 2) : n);


Где Y определен, к примеру, вот так:

static Func<A, R> Y<A, R>(Func<Func<A, R>, Func<A, R>> f)
{
  Func<A, R> g = null;
  g = f(a=>g(a));
  return g;
}


Для больших N, такой подход непродуктивен. Как посчитать число N быстрее? Для начала давайте вспомним, почему например x*x считается быстрее чем Math.Pow(x,2)? Потому что для целочисленной степени можно не только обойтись без рядов Тейлора, но также оптимизировать вычисление для больших степеней путем создания временных переменных. Например, x⁴ можно считать как int y = x * x; return y * y; – и чем больше степень, тем больше экономия.

К чему это я? К тому что число Фибоначчи можно рассчитать с помощью следующей формулы:




Теперь понятно зачем нам целочисленное возведение в степень? Для матриц можно делать то же самое, только сначала нужно понять как это вообще делается. На просторах интеренета я нашел возможно идеальный алгоритм, который оптимизирует целочисленное возведение в степень. Прошу заметить, что в примере ниже степень имеет тип short.

public static long IntPower(int x, short power)
{
  if (power == 0) return 1;
  if (power == 1) return x;
  int n = 15;
  while ((power <<= 1) >= 0) n--;
  long tmp = x;
  while (--n > 0)
    tmp = tmp * tmp *
         (((power <<= 1) < 0) ? x : 1);
  return tmp;
}


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

struct mtx2x2
{
  public int _11, _12, _21, _22;
  public static mtx2x2 operator*(mtx2x2 lhs, mtx2x2 rhs)
  {
    return new mtx2x2
             {
               _11 = lhs._11*rhs._11 + lhs._12*rhs._21,
               _12 = lhs._11*rhs._12 + lhs._12*rhs._22,
               _21 = lhs._21*rhs._11 + lhs._22*rhs._21,
               _22 = lhs._21*rhs._12 + lhs._22*rhs._22
             };
  }
}


После этого нужно было определить две константы которые нам пригодятся – ту матрицу которую мы будет возводить в степень и единичную матрицу:

private static readonly mtx2x2 fibMtx = new mtx2x2 {_11 = 1, _12 = 1, _21 = 1};
private static readonly mtx2x2 identity = new mtx2x2 {_11 = 1, _22 = 1};


Теперь мы можем переписать метод IntPower() для матриц 2×2:

public static mtx2x2 IntPower(mtx2x2 x, short power)
{
  if (power == 0) return identity;
  if (power == 1) return x;
  int n = 15;
  while ((power <<= 1) >= 0) n--;
  mtx2x2 tmp = x;
  while (--n > 0)
    tmp = (tmp * tmp) * (((power <<= 1) < 0) ? x : identity);
  return tmp;
}


И определить новый метод для вычисления числа Фибоначчи:

static int fibm(short n)
{
  return IntPower(fibMtx, (short)(n-1))._11;
}


Вот и все. Думаю что сравнивать производительность бессмысленно – у меня fib(40) считается 4 секунды, а fibm(40) выводится моментально. ■