FN - Editorial

[djdolls]

PROBLEM LINK:

Practice
Contest

Author: Tom Chen
Tester: Jingbo Shang
Editorialist: Ajay K. Verma

DIFFICULTY:

MEDIUM-HARD

PREREQUISITES:

Quadratic residue, Tonelli-Shanks algorithm, Baby-step giant-step algorithm

PROBLEM:

Given a prime P of the form 5k + 1 or 5k - 1, and and integer C, find, if exists, the smallest non-negative integer n such that F_n = C (mod P), where F_n is the n-th fibonacci number.

QUICK EXPLANATION:

Based on the condition that (P mod 5) is either 1 or -1, we can deduce that 5 is a quadratic residue modulo P. Tonelli-Shanks algorithm can be used to find the integer x such that x² = 5 (mod P). Using the value of x, we can find the explicit formula for n-th fibonacci number modulo P.

// Golden ratio modulo P
y = (1 + x) / 2
F_n = (yⁿ - (-1/y)ⁿ) / x

Since we know the values of F_n and x, the above equation can be used to find the value of yⁿ by solving a quadratic equation modulo P.

Now, we know the value of y and yⁿ, and want to calculate the value of n. This is the classical Discrete Log problem and can be solved efficiently using Baby-step giant-step algorithm.

EXPLANATION:

The value of the n-th fibonacci number can be calculated using the explicit formula:
F_n = (yⁿ - (-1/y)ⁿ) / √5,
where y is the golden ratio, i.e.,
y = (1 + √5) / 2

Let us see, how this formula looks for fibonacci number modulo a prime P.
First, we need to check if √5 exists modulo P, i.e., if there is an integer x such that x² = 5 (mod P).

Legendre Symbol and golden ratio computation:

In order to check whether for a given integer a, there exist an integer x such that x² = a (mod P), one needs to compute the legendre symbol (a | P). Here, we want to compute (5 | P).

Using quadratic reciprocity, we can write
(5 | P) = (-1)^{((5 - 1)/2 * (P - 1)/2)} * (P | 5) = (P | 5)
= (P mod 5 | 5)

Based on the constraints (P mod 5) is either 1 or 4, both of them are a quadratic residue modulo 5 (This is because 1² = 1 (mod 5), and 2² = 4 (mod 5)). Hence, (5 | P) = (P mod 5 | P) = 1, i.e., 5 is a quadratic residue modulo P.

This means that there exist an integer x such that x² = 5 (mod P).
We will discuss later how to compute such x. For the time being let us assume that we have an oracle who computes the square root modulo a prime for us.

Hence, the expresison for golden ratio can be written as
y = (1 + x)/2 mod P.

Since P is an odd prime, 2 * (P + 1)/2 = 1 (mod 1), i.e., 2^-1 = (P + 1)/2. Therefore,
y = (P + 1)/2 * (1 + x).

Note that, (1 + x)/2 * (-1 + x)/2 = (x² - 1)/4 = 1 mod P.
Hence, y^-1 mod P exists.

Explicit formula for fibonacci number and bound on the answer (if exists):

Now that we know the value of golden ratio y, and √5 = x, we can write
Fⁿ = (yⁿ - (-1/y)ⁿ) / x

Note that x and P are relatively prime (this is because gcd(x² = 5, P) = 1). Hence, x^-1 mod P must exist. Therefore, we can write the expression for n-th fibonacci number explicitly as

Fⁿ = x^-1 * (yⁿ - (-1/y)ⁿ).

According to Fermat’s little theorem (a^P = a (mod P)). Hence, both yⁿ and (-1/y)ⁿ are periodic, with a period of at most P.
This implies that F_n mod P is also periodic with period at most P. Hence, for a given C, if there exist an n such that F_n = C, then also once such n exist, which is <= P. Hence, we only look for candidate n which are smaller than P.

Reduction into discrete log problem:

Now, we want the n-th fibonacci number to be C. This gives us the following equation:
C = x^-1 * (yⁿ - (-1/y)ⁿ)
C * x = yⁿ - (-1/y)ⁿ

Now we should consider two cases: One when n is odd, and the other one when n is even. For each case we should solve the above equation and compute the value of n. If the parity of computed n matches with the case that we are considering, then we have a solution, otherwise we do not have an n at which fibonacci sequence takes value C modulo P.

Here, we only consider the case when n is even, and leave the other one as an exercise for the reader. In case of even n
yⁿ - (1/y)ⁿ = C * x = z

This is a quadratic equation in yⁿ. Solving this equation we get
yⁿ = (z ± √(z² + 4))/2.

Once again we use our oracle to compute the square root of (z² + 4) if exists. We have two possible values of yⁿ, we consider both values and compute the value of n from it.

In other words, now our problem has been reduced to the following problem: Given integer y, u and prime P, find the smallest even number n such that
yⁿ = u mod P

If we remove the constraint that n is even, then this problem is known as discrete log problem. In the next section, we discuss how to solve the discrete log problem, and the following section shows how to extend the discrete log problem when the extra constraint of n being even is enforced.

Discrete log problem (Baby step giant step algorithm):

We know that if the equation (yⁿ = u) has a solution, it must have a solution smaller than P. If we have an integer Q whose value is approximately √P, then we can write the integer n as
n = Qa + b, where both a and b are smaller than Q.

yⁿ = y^{Qa + b}
u = y^{Qa + b}
u * (y^-Q)^a = y^b

We compute the values of y⁰, y¹, …, y^{Q - 1} and store them in a map. This can be done in O (√P) time. If the map contains u, then we have already found a solution of the above equation. Otherwise we compute w = y^-Q, and iterate through its powers.

For the a-th power of w we have w^a = y^-Qa
Next we compute u * w^a, and search for it in the map, if this value is found in the map, that means we have y^b = u * y^-Qa, i.e., y^{Qa + b} = u, which is a solution of the equation we are looking for.

If we have iterated through all values of a from 1 to Q, and did not find a solution, that means the equation has no solution.

Family of solution for yⁿ = u mod P:

Let us say that we have two solutions n1 and n2 of the above equation, that means
yⁿ¹ = yⁿ² = u mod P
y^{(n1 - n2)} = 1 mod P
i.e., (n1 - n2) is a multiple of the order of y, Ord(y).

Ord(y) is the smallest integer m such that y^m = 1 (mod P). Based on Fermat’s little theorem y^{P - 1} = 1, Hence Ord(y) must divide (P - 1).
We can iterate thorough divisors of (P - 1), and for each divisor d check if y^d = 1. The smallest such d will be Ord(y).

Since there are O (√P) divisors of (P - 1), and computation of (y^d mod P) can be done in O (lg P) time, we can compute Ord(y) in O (√P lg P) time.

Once we have computed Ord(y) = m, the family of solution of the above equation can be represented as
n = n1 + k * m
where n1 is the smallest solution found by Baby step giant step algorithm.

If n1 is even, then we already have a solution of the discrete log problem with even constraint. If n1 is odd, and m is even, then we can see that the family of solution will consist of odd integers only. In this case we have no even solution of the discrete log. On the other hand, if both n1 and m is odd, then the smallest even solution of the discrete log problem is (n1 + m).

Square root modulo a prime computation:

Finally we discuss the oracle that we have been using to compute the square root modulo a prime P. This is in fact the Tonelli Shanks algorithm. The algorithm is well explained here Tonelli Shanks algorithm, and requires finding a quadratic nonresidue modulo the prime. Since there are (P + 1)/2 quadratic non-residues, a randomized algorithm will find find a quadratic non-reside in O (1) steps with high probability.

The overall complexity of the Tonelli Shanks algorithm is O (lg² P).

Time Complexity:

O (√P lg P)

AUTHOR’S AND TESTER’S SOLUTIONS:

Author’s solution will be put up soon.
Tester’s solution can be found here.

[shashanktandon]

Can anyone now tell where my code fails and can explain me the concept what is wrong in my concept??? Link to my submission is http://www.codechef.com/viewsolution/2825396

[shashanktandon]

Can anyone now tell where my code fails and can explain me the concept what is wrong in my concept??? Link to my submission is http://www.codechef.com/viewsolution/2825396

[sy2006]

As the bsgs algo can be performed in O(√P) time as it was said above, shouldn’t the Time Complexity be O(√P)?

[nitinj]

HOLYYYYYYYYYYYYYYYYYYYYYY MOTHER OF CODER GOD! THAT QUESTION…

[kevinsogo]

Baby-step giant-step can be implemented this way using a map/associative array (taken from wikipedia):

discrete_log(a, b, p): // finds the discrete log of b base a modulo p
    // initialization
    m = ceiling(√p)
    store = ((new map))
    a_pow = 1
    for j = 0..m-1:
        store[a_pow] = j
        a_pow = (a_pow·a) mod p
    a_pow = mod_inverse(a_pow, p)
    // queries
    y = b
    for i = 0..m-1:
        if store has key y:
            j = store[y]
            return i·m + j
        y = (y·a_pow) mod p
    return None

Now, when solving FN we could potentially use this function 4 times (two quadratic equations, two roots each), with the same values a and p, so we could save a few operations if we do the first part (initializing m, store and a_pow) only once.

Now with a hash_map we can guarantee an O(√p) running time, but I decided to use something else entirely… sorting.

One can remove the use of a map by simply sorting the key-value pairs and performing a binary search on it when accessing. This works because all accesses to the map are done after all insertions.

discrete_log(a, b, p): // finds the discrete log of b base a modulo p
    // initialization
    m = ceiling(√p)
    store_pairs = ((new list))
    a_pow = 1
    for j = 0..m-1:
        add the pair (a_pow, j) to store_pairs
        a_pow = (a_pow·a) mod p
    sort store_pairs
    a_pow = mod_inverse(a_pow, p)
    // queries
    y = b
    for i = 0..m-1:
        find y in store_pairs using binary search
        if y is found:
            j = ((pair of y in store_pairs))
            return i·m + j
        y = (y·a_pow) mod p
    return None

This makes an O(√p log p) running time because of sorting and binary search. To hack off the O(log p) part, we use some tricks. First, let’s also sort the query values from the second loop, and perform some kind of merge (similar to merge sort).

discrete_log(a, b, p): // finds the discrete log of b base a modulo p
    // initialization
    m = ceiling(√p)
    store_pairs = ((new list))
    a_pow = 1
    for j = 0..m-1:
        add the pair (a_pow, j) to store_pairs
        a_pow = (a_pow·a) mod p
    sort store_pairs
    a_pow = mod_inverse(a_pow, p)
    // collect query pairs
    query_pairs = ((new list))
    y = b
    for i = 0..m-1:
        add the pair (y, i) to query_pairs
        y = (y·a_pow) mod p
    sort query_pairs
    // merge part
    L = R = 0
    while L < store_pairs.length and R < query_pairs.length:
        if store_pairs[L].first < query_pairs[R].first:
            L++
        else if store_pairs[L].first > query_pairs[R].first:
            R++
        else: // now equal
            j = store_pairs[L].second
            i = query_pairs[R].second
            return i·m + j
    return None

This removes the need for binary search, but this is still O(√p log p) because of sorting. To improve that, we perform radix sort, by taking advantage of the fact that the first component is in the range [0,p). We use m = ceiling(√p) as our base so we will only need two passes.

sort(L, m): // L is a list of ordered pairs of length m
    // initialize buckets
    for i in 0..m-1:
        buckets[i] = ((new list))

    // first pass
    for (a,b) in L:
        add the pair (a,b) to buckets[a%m]
    M = ((new list))
    for i in 0..m-1:
        for (a,b) in buckets[i]:
            add the pair (a,b) to M
        empty the list buckets[i]

    // second pass
    for (a,b) in M:
        add the pair (a,b) to buckets[a/m]
    empty the list L
    for i in 0..m-1:
        for (a,b) in buckets[i]:
            add the pair (a,b) to L
        empty the list buckets[i]

    return L

Now discrete_log runs in O(√p) time!

Of course we can still preprocess store_pairs and a_pow beforehand to save a few operations.

[kevinsogo]

The time complexity is reported as O(√p lg p) because they need to calculate Ord(y).

I get around that by simply finding all discrete logs and testing them all for solution.

[iscsi]

My solution very similar but based on matrices:

so, we would like to find an n where F_n is C, and F_n-1 D (what we don’t know), and of course F_n+1 is D+C.
But if we use Cassini’s identity:
and use C,D,and C+D we’ve got two a quadratic equation mod p (if n is even or not). If we solve we’ve got some candidate D and if we have a candidate we can find it is a solution or not and the indice of the power with baby step giant step on matrices…

[sy2006]

@kevinsogo In this problem we only need to find the smallest solution not all of them so I think there is no need for calculating ord(y), isn’t it? Just bsgs with hash table is enough.

[kevinsogo]

Calculating Ord(y) isn’t that slow in practice because most numbers n have much less than √n factors anyway.

But yeah, it isn’t really needed.