Fibonacci Fever (Çözgüdiň düşündirişi)

Since $n$ can be as large as 10^{18}, iterating over the Fibonacci numbers one by one is impossible. The matrix method does not work either: the sum of $k$-th powers would require a matrix of size about $k$, and $k$ can be as large as 10^5. Therefore, we use an explicit formula for the Fibonacci numbers — Binet's formula:

where

The numbers $\phi$ and $\psi$ are the roots of the equation x^2 = x + 1, and by Vieta's formulas \varphi\psi = -1.

Expanding the power. Raise Binet's formula to the power $k$ and expand it using the binomial theorem:

Denote

Sum over $i$ from $1$ to $n$ and change the order of summation:

Thus, the sum of $n$ terms is reduced to a sum of k + 1 terms.

Geometric progression. The inner sum is the sum of a geometric progression with common ratio t = t_j:

The case t = 1 really occurs and must be handled separately (otherwise a division by zero occurs). For example, if $k$ is divisible by $4$, then for j = k / 2 we get

The power t^n is computed by binary exponentiation in O(\log n) time, so each term is found quickly even for n = 10^{18}.

Computing the ratios iteratively. There is no need to raise $\phi$ and $\psi$ to a power anew for each $j$. Consecutive ratios are related by

so it is enough to compute t_0 = \varphi^k once, and then multiply the current value by $\psi /\phi$ at each iteration.

Modular computations. Floating-point arithmetic is not suitable here: we need the exact answer modulo p = 10^9 + 7. We will work with numbers of the form a + b\sqrt{5}, where $a$ and $b$ are residues modulo $p$. Addition and multiplication of such numbers follow the rules

Numbers of the form a + b\sqrt{5} form the field \mathbb{Z}_p[\sqrt{5}]. In this field, every nonzero element is invertible:

where the denominator a^2 - 5b^2 is nonzero modulo $p$. This guarantees that all divisions in the solution are valid. Division by $2$ in the formulas for $\phi$ and $\psi$ is replaced by multiplication by 2^{-1} \bmod p = (p + 1) / 2.

The inverse of 2. The inverse of $2$ modulo $p$ is a number $x$ such that 2x \equiv 1 \pmod{p}.

Let us check that x = (p + 1) / 2 works:

• p = 10^9 + 7 is an odd prime, so p + 1 is even and (p + 1) / 2 is an integer. It is less than $p$, so it is an ordinary residue modulo $p$.

• Multiply it by $2$: 2 \cdot (p + 1) / 2 = p + 1 \equiv 1 \pmod{p}.

Hence, (p + 1) / 2 is exactly 2^{-1} \bmod p.

The required sum is an integer, so in the final value a + b\sqrt{5} the coefficient $b$ is necessarily zero, and the answer is $a$.

Example

Consider the first example: $\mathbf{n}=4$, \boldsymbol{k = 1}. For k = 1 the formula becomes

Let us compute the powers of $\phi$ using the relation \varphi^2 = \varphi + 1:

Their sum equals \dfrac{15 + 7\sqrt{5}}{2}. The powers of $\psi$ are obtained by replacing \sqrt{5} with -\sqrt{5}, so their sum equals \dfrac{15 - 7\sqrt{5}}{2}. Therefore,

which matches the answer: 1 + 1 + 2 + 3 = 7.

Consider an example for \boldsymbol{k = 3}. Expand the cube of the difference:

The common ratios of the progressions are (using \varphi\psi = -1):

Then

Each of the four sums is computed using the geometric progression formula. Let us check this for n = 3. The powers of $\phi$ are

Then

The sums for $\psi$ are obtained by replacing \sqrt{5} with -\sqrt{5}:

Substituting the values found:

The result matches the direct computation:

Consider an example for \boldsymbol{k = 4} (a special case). Here

and

The middle term corresponds to the common ratio t_2 = 1: the geometric progression formula does not apply to it, and the sum 1 + 1 + \ldots + 1 equals $n$. This is exactly the case handled separately in the program. For example, for n = 3 we get:

Algorithm implementation

Declare a constant for the modulus $MOD$.

#define MOD 1000000007

We'll store the values of binomial coefficients in the array $Cnk$.

vector<long long> Cnk;

To implement arithmetic in the field Z_{MOD} [\sqrt{5}], define a structure $Num$. This structure represents numbers of the form x + y \sqrt{5}.

struct Num
{
  long long x, y;
  Num(long long x = 0, long long y = 0) : x(x% MOD), y(y% MOD) {}
};

The function add implements addition in the field Z_{MOD} [\sqrt{5}].

Num add(const Num& a, const Num& b)
{
  return Num((a.x + b.x) % MOD, (a.y + b.y) % MOD);
}

The function mult implements multiplication in the field Z_{MOD} [\sqrt{5}].

Num mult(const Num& a, const Num& b)
{
  return Num(
    (a.x * b.x + 5 * a.y % MOD * b.y) % MOD,
    (a.x * b.y + a.y * b.x) % MOD
  );
}

The function pw implements binary exponentiation in the field Z_{MOD} [\sqrt{5}], where the exponent $exp$ is an integer.

Num pw(Num base, long long exp)
{
  Num res(1, 0);
  while (exp)
  {
    if (exp & 1) res = mult(res, base);
    base = mult(base, base);
    exp >>= 1;
  }
  return res;
}

The function pw implements fast binary exponentiation to compute a^b modulo $MOD$.

long long pw(long long a, long long b)
{
  long long res = 1;
  a %= MOD;
  while (b)
  {
    if (b & 1) res = (res * a) % MOD;
    a = a * a % MOD;
    b >>= 1;
  }
  return res;
}

The function inv computes the multiplicative inverse of an element $a$ in the field Z_{MOD} [\sqrt{5}].

If

then its inverse element is given by

Num inv(const Num& a)
{
  long long den = (a.x * a.x - 5 * a.y * a.y) % MOD;
  if (den < 0) den += MOD;

  long long invDen = pw(den, MOD - 2);

  return Num(
    a.x * invDen % MOD,
    (MOD - a.y) * invDen % MOD
  );
}

The function divide implements division in the field Z_{MOD} [\sqrt{5}]:

Num divide(const Num& a, const Num& b)
{
  return mult(a, inv(b));
}

The main part of the program. Read the input data.

scanf("%lld %lld", &n, &k);

Compute the binomial coefficients C_k^0, C_k^1, C_k^2, ..., C_k^k and store them in the array $Cnk$. Note that

Cnk.resize(k + 1);
Cnk[0] = 1;
for (i = 1; i <= k; i++)
{
  Cnk[i] = Cnk[i - 1] * (k - i + 1) % MOD;
  Cnk[i] = Cnk[i] * pw(i, MOD - 2) % MOD;
}

Compute inv2 = 2^{-1} \mod MOD = (MOD + 1) / 2.

long long inv2 = (MOD + 1) / 2;

Initialize the constants: phi = \phi = \frac{1 + \sqrt{5}}{2} and psi = \psi = \frac{1 - \sqrt{5}}{2}.

Num phi(inv2, inv2);
Num psi(inv2, MOD - inv2);

Compute ratio = \psi / \phi.

Num ratio = divide(psi, phi);

Initialize the constants $1$ and $-1$ in the field Z_{MOD} [\sqrt{5}].

Num one(1, 0);
Num minus_one(MOD - 1, 0);

The answer will be accumulated in the variable $res$.

Num res(0, 0);

Initialize t = \phi^k. Then, in the loop, we will multiply this value by $ratio$, obtaining

Num t = pw(phi, k);

In the variable $res$ we compute the sum

The required sum contains k + 1 terms. Therefore, we perform k + 1 iterations.

for (j = 0; j <= k; j++)
{
  Num term;

At this point t = \phi^{k-j} \psi^j. Compute the sum of the geometric progression:

However, the case t = 1 should be handled separately, because the formula for the sum of a geometric progression contains division by $0$. When t = 1 the sum takes the form: 1 + 1 + 1 + ... + 1 ($n$ times) and is equal to $n$.

  if (t.x == 1 && t.y == 0)
    term = Num(n % MOD, 0);
  else
  {
    Num tn = pw(t, n);
    Num num = add(tn, minus_one);     // t^n - 1
    term = mult(t, num);              // t * (t^n - 1)
    Num den = add(t, minus_one);      // t - 1
    term = divide(term, den);         // t * (t^n - 1) / (t – 1)
  }

Compute term = (-1)^j C_k^j \cdot term.

  long long coef = Cnk[j];
  if (j & 1) coef = MOD - coef;

  term.x = term.x * coef % MOD;
  term.y = term.y * coef % MOD;

Add the next term $term$ to the result $res$.

  res = add(res, term);

Multiply $t$ by $ratio$, obtaining t = \phi^{k-j-1} \psi^{j + 1}.

  t = mult(t, ratio);
}

It remains to divide $res$ by (\sqrt{5})^k.

Num sqrt5(0, 1);
Num sqrt5k = pw(sqrt5, k);
res = divide(res, sqrt5k);

Print the answer.

printf("%lld\n", res.x);