Fibonacci numbers (Editorial)

It is known that Fibonacci numbers satisfy the following identity:

GC D (F_{n}, F_{m}) = F_{GC D (n, m)}

GCD(F_n, F_m) = F_{GCD(n,m)}

This means that the problem reduces to computing the $k$ -th Fibonacci number modulo $1 0^{8}$ 10^8, where

k = GC D (n, m)

k = GCD(n, m)

Since $n, m \leq 1 0^{18}$ n, m \le 10^{18}, we have $k \leq 1 0^{18}$ k ≤ 10^{18}. Therefore, it is necessary to compute the value of $F_{k} m o d 1 0^{8}$ F_k~mod~10^8 in $O (l o g_{2} k)$ O(log_2k) time.

Theorem. To efficiently compute Fibonacci numbers, it is convenient to use matrix exponentiation. It is known that:

(1110)^{k} = (F_{k + 1} F_{k} F_{k} F_{k - 1})

\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^k = \begin{pmatrix} F_{k+1} & F_k \\ F_k & F_{k-1} \end{pmatrix}

Base case. For $k = 1$ k = 1, we have:

(1110) = (F_{2} F_{1} F_{1} F_{0})

\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} F_2 & F_1 \\ F_1 & F_0 \end{pmatrix}

which is true since $F_{0} = 0, F_{1} = 1, F_{2} = 1$ F_0 = 0, F_1 = 1, F_2 = 1.

Inductive step. Assume that the formula holds for $k$ . Then, for $k + 1$ k + 1:

(1110)^{k + 1} = (F_{k + 1} F_{k} F_{k} F_{k - 1}) (1110) = (F_{k + 1} + F_{k} F_{k} + F_{k - 1} F_{k + 1} F_{k}) = (F_{k + 2} F_{k + 1} F_{k + 1} F_{k})

\begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix}^{k+1} = \begin{pmatrix}F_{k+1} & F_k \\ F_k & F_{k-1}\end{pmatrix} \begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix} = \begin{pmatrix}F_{k+1}+F_k & F_{k+1} \\ F_k+F_{k-1} & F_k\end{pmatrix} = \begin{pmatrix}F_{k+2} & F_{k+1} \\ F_{k+1} & F_k\end{pmatrix}

Thus, the formula is valid for all $k \geq 1$ k \ge 1.

It remains to implement raising a matrix to the power $k$ using fast exponentiation in $O (l o g_{2} k)$ O(log_2k) time.

Example

Fibonacci numbers are defined as follows:

Let's compute the Fibonacci numbers using matrix exponentiation:

(1110)^{2} = (1110) \cdot (1110) = (2111) = (F_{3} F_{2} F_{2} F_{1})

\begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix}^2 = \begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix} \cdot \begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix} = \begin{pmatrix}2 & 1 \\ 1 & 1\end{pmatrix} = \begin{pmatrix}F_3 & F_2 \\ F_2 & F_1\end{pmatrix}

(1110)^{4} = (2111)^{2} = (2111) \cdot (2111) = (5332) = (F_{5} F_{4} F_{4} F_{3})

\begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix}^4 = \begin{pmatrix}2 & 1 \\ 1 & 1\end{pmatrix}^2 = \begin{pmatrix}2 & 1 \\ 1 & 1\end{pmatrix} \cdot \begin{pmatrix}2 & 1 \\ 1 & 1\end{pmatrix} = \begin{pmatrix}5 & 3 \\ 3 & 2\end{pmatrix} = \begin{pmatrix}F_5 & F_4 \\ F_4 & F_3\end{pmatrix}

(1110)^{8} = (5332)^{2} = (5332) \cdot (5332) = (34212113) = (F_{9} F_{8} F_{8} F_{7})

\begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix}^8 = \begin{pmatrix}5 & 3 \\ 3 & 2\end{pmatrix}^2 = \begin{pmatrix}5 & 3 \\ 3 & 2\end{pmatrix} \cdot \begin{pmatrix}5 & 3 \\ 3 & 2\end{pmatrix} = \begin{pmatrix}34 & 21 \\ 21 & 13\end{pmatrix} = \begin{pmatrix}F_9 & F_8 \\ F_8 & F_7\end{pmatrix}

Algorithm implementation – operators overloading

Let us declare the constant $M O D$ , modulo which all computations will be performed.

#define MOD 100000000

The function gcd computes the greatest common divisor of two numbers $a$ and $b$ .

long long gcd(long long a, long long b)
{
  return (!b) ? a : gcd(b, a % b);
}

Let us declare a Matrix class and write its constructor.

class Matrix
{
public:
  long long a, b, c, d;
  Matrix(long long a = 1, long long b = 0, 
         long long c = 0, long long d = 1) 
  {
    this->a = a; this->b = b;
    this->c = c; this->d = d;
  }

C++

Let us overload the matrix multiplication operator. All computations are performed modulo $M O D = 1 0^{8}$ MOD = 10^8.

  Matrix operator* (const Matrix &x)
  {
    Matrix res;
    res.a = (a * x.a + b * x.c) % MOD;
    res.b = (a * x.b + b * x.d) % MOD;
    res.c = (c * x.a + d * x.c) % MOD;
    res.d = (c * x.b + d * x.d) % MOD;
    return res;
  }

C++

Next, overload the operator for raising a matrix to the power $n$ . The time complexity of the algorithm is $O (l o g_{2} n)$ O(log_2n).

  Matrix operator^ (long long n)
  {
    Matrix x(*this);
    if (n == 0) return Matrix();
    if (n & 1) return x * (x ^ (n - 1));
    return (x * x) ^ (n/2);
  }
};

C++

The function fib returns the $n$ -th Fibonacci number $F_{n}$ F_n modulo $1 0^{8}$ 10^8.

long long fib(long long n)
{
  Matrix res(1,1,1,0);
  res = res ^ n;
  return res.b;
}

The main part of the program. Read the input data, compute, and print the value of $F_{GC D (n, m)}$ F_{GCD(n, m)}.

while(scanf("%lld %lld",&n,&m) == 2)
{
  d = gcd(n,m);
  printf("%lld\n",fib(d));
}

Algorithm implementation – functions

Let us declare the constant $M O D$ , modulo which all computations will be performed.

#define MOD 100000000

The function gcd computes the greatest common divisor of two numbers $a$ and $b$ .

long long gcd(long long a, long long b)
{
  return (!b) ? a : gcd(b, a % b);
}

Declare the Matrix structure — a $2 \times 2$ 2 \times 2 matrix. By default, the identity matrix is created, which is convenient for exponentiation.

struct Matrix 
{
  long long a, b, c, d;
  Matrix(long long a_ = 1, long long b_ = 0,
         long long c_ = 0, long long d_ = 1)
         : a(a_), b(b_), c(c_), d(d_) {}
};

C++

The function multiply multiplies two $2 \times 2$ 2 \times 2 matrices using the standard formula.

Matrix multiply(Matrix &x, Matrix &y) 
{
  return Matrix(
    (x.a * y.a + x.b * y.c) % MOD,
    (x.a * y.b + x.b * y.d) % MOD,
    (x.c * y.a + x.d * y.c) % MOD,
    (x.c * y.b + x.d * y.d) % MOD
  );
}

C++

The function power implements binary exponentiation of the matrix $ba se$ to the power $e x p$ .

Matrix power(Matrix base, long long exp) 
{
  Matrix result;
  while (exp > 0)
  {
    if (exp & 1)
      result = multiply(result, base);
    base = multiply(base, base);
    exp >>= 1;
  }
  return result;
}

C++

The function fib computes the $n$ -th Fibonacci number modulo $1 0^{8}$ 10^8 using the matrix method.

long long fib(long long n) 
{
  if (n == 0) return 0;
  Matrix m(1, 1, 1, 0);
  Matrix res = power(m, n);
  return res.b; // F(n)
}

C++

The main part of the program. Read the input data, compute, and print the value of $F_{GC D (n, m)}$ F_{GCD(n, m)}.

while (scanf("%lld %lld", &n, &m) == 2) 
{
  d = gcd(n, m);
  printf("%lld\n", fib(d));
}

Algorithm implementation – memoization

The following identity is known for Fibonacci numbers:

F_{n + m} = F_{m} \cdot F_{n + 1} + F_{m - 1} \cdot F_{n}

F_{n+m} = F_m \cdot F_{n+1} + F_{m-1} \cdot F_n

From it, the following special cases directly follow:

if $m = n$ m = n, then $F_{2 n} = F_{n} \cdot F_{n + 1} + F_{n - 1} \cdot F_{n}$ F_{2n} = F_n \cdot F_{n+1} + F_{n-1} \cdot F_n.
if $m = n + 1$ m = n + 1, then $F_{2 n + 1} = F_{n + 1} \cdot F_{n + 1} + F_{n} \cdot F_{n}$ F_{2n+1} = F_{n+1} \cdot F_{n+1} + F_n \cdot F_n.

These formulas allow Fibonacci numbers to be computed using a "divide and conquer" approach, reducing the computation of $F_{n}$ F_n to values with indices approximately half as large.

The base cases are handled separately:

F_{0} = 0, F_{1} = F_{2} = 1

F_0 = 0, F_1 = F_2 = 1

With this approach, the recursion depth is proportional to $l o g_{2} n$ log_2n, and the time complexity of the algorithm is $O (l o g_{2} n)$ O(log_2n).

Declare the constant $M O D$ , modulo which all computations will be performed.

#define MOD 100000000

Declare a variable $F$ of type map to store (memoize) the Fibonacci numbers that have already been computed.

map<long long, long long> F;

The function gcd computes the greatest common divisor of two numbers $a$ and $b$ .

long long gcd(long long a, long long b)
{
  return (!b) ? a : gcd(b,a % b);
}

The function fib computes the $n$ -th Fibonacci number modulo $1 0^{8}$ 10^8.

long long fib(long long n)
{

The base cases are handled separately.

  if (n == 0) return 0;
  if (n == 1) return 1;
  if (n == 2) return 1;

If the value of $F_{n}$ F_n has already been computed and stored in $F$ , it is returned immediately without recomputation.

  if (F[n]) return F[n];

Next, we use the decomposition of the number $n$ into:

$n = 2 k + 1$ n = 2k + 1 — odd case: $F_{n} = F_{k}^{2} + F_{k + 1}^{2}$ F_n=F_k^2 + F_{k+1}^2.
$n = 2 k$ n = 2k — even case: $F_{n} = F_{k} \cdot F_{k + 1} + F_{k - 1} \cdot F_{k}$ F_n = F_k \cdot F_{k+1} + F_{k-1} \cdot F_k.

  long long k = n / 2;

Store and return the result.

  if (n % 2 == 1) // n = 2*k + 1
    return F[n] = (fib(k) * fib(k) + fib(k+1) * fib(k+1)) % MOD;
  else // n = 2*k
    return F[n] = (fib(k) * fib(k+1) + fib(k-1) * fib(k)) % MOD;
}

The main part of the program. Read the input data, compute, and print the value of $F_{GC D (n, m)}$ F_{GCD(n, m)}.

while(scanf("%lld %lld",&a,&b) == 2)
{
  d = gcd(a,b);
  printf("%lld\n",fib(d));
}

Java implementation

import java.util.*;

public class Main
{
  static HashMap<Long, Long> F = new HashMap<Long, Long>();
  static long MOD = 100000000;
  
  static long gcd(long a, long b)
  {
    return (b == 0) ? a : gcd(b,a % b);
  }

  static long fib(long n)
  {
    if (n == 0) return 0;
    if (n == 1) return 1;
    if (n == 2) return 1;
    if (F.containsKey(n)) return F.get(n);

    long k = n / 2;
    if (n % 2 == 1)
    {
      F.put(n, (fib(k) * fib(k) + fib(k+1) * fib(k+1)) % MOD);
      return F.get(n);
    }
    else
    {
      F.put(n, (fib(k) * fib(k+1) + fib(k-1) * fib(k)) % MOD);    	
      return F.get(n);
    }
  }
  
  public static void main(String[] args)
  {
    Scanner con = new Scanner(System.in);
    while(con.hasNext())
    {
      long a = con.nextLong();
      long b = con.nextLong();
      long d = gcd(a,b);
      System.out.println(fib(d));  
    }
    con.close();
  }
}

Java