FN - Editorial

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 Fn = C (mod P), where Fn 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 x2 = 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
Fn = (yn - (-1/y)n) / x

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

Now, we know the value of y and yn, 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:
Fn = (yn - (-1/y)n) / √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 x2 = 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 x2 = 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 12 = 1 (mod 5), and 22 = 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 x2 = 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 = (x2 - 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
Fn = (yn - (-1/y)n) / x

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

Fn = x-1 * (yn - (-1/y)n).

According to Fermat’s little theorem (aP = a (mod P)). Hence, both yn and (-1/y)n are periodic, with a period of at most P.
This implies that Fn mod P is also periodic with period at most P. Hence, for a given C, if there exist an n such that Fn = 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 * (yn - (-1/y)n)
C * x = yn - (-1/y)n

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
yn - (1/y)n = C * x = z

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

Once again we use our oracle to compute the square root of (z2 + 4) if exists. We have two possible values of yn, 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
yn = 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 (yn = 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.

yn = yQa + b
u = yQa + b
u * (y-Q)a = yb

We compute the values of y0, y1, …, yQ - 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 wa = y-Qa
Next we compute u * wa, and search for it in the map, if this value is found in the map, that means we have yb = u * y-Qa, i.e., yQa + 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 yn = u mod P:

Let us say that we have two solutions n1 and n2 of the above equation, that means
yn1 = yn2 = 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 ym = 1 (mod P). Based on Fermat’s little theorem yP - 1 = 1, Hence Ord(y) must divide (P - 1).
We can iterate thorough divisors of (P - 1), and for each divisor d check if yd = 1. The smallest such d will be Ord(y).

Since there are O (√P) divisors of (P - 1), and computation of (yd 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 (lg2 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.

7 Likes

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 CodeChef: Practical coding for everyone

HOLYYYYYYYYYYYYYYYYYYYYYY MOTHER OF CODER GOD! THAT QUESTION…

1 Like

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.

7 Likes

My solution very similar but based on matrices:

alt text

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:
alt text 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…

5 Likes

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 CodeChef: Practical coding for everyone

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

1 Like

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.

1 Like

@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.

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.

For this test case:
3 1999996771

Tester’s code gives the answer -1, but n=4 is a correct answer. Why is that?