Let Bus \(j\) be the first bus Takahashi will take after he starts his travel in City \(Y_i\) at time \(X_i\). If there is no such bus, or \(Z_i \leq S_j\), then obviously Takahashi is in City \(Y_i\) at time \(Z_i\). On the other hand, if \(S_j \lt Z_i \leq T_j\), then Takahashi is on Bus \(j\) at time \(Z_i\).
Hereinafter, we will consider the other cases.

Now, let us consider binary-searching for the positive integer \(k\) satisfying the following conditions:

• Suppose that Takahashi has taken \(k\) buses, starting from Bus \(j\). Then, the time \(k'\) that Takahashi gets off the last bus is less than or equal to \(Z\).

Once we have obtained such \(k\), we can easily find which bus Takahashi is on, or which city Takahashi is in, with a simple case divisions.

In order to solve the subproblem, what information we need is “which is the last bus he took.” This information can be obtained with a technique called “doubling” in an \(O(\log M)\) time each, after a precalculation of \(O(M \log M)\) time.

Specifically, we can compute it as follows.

First, for each \(i=1,2,\ldots,M\) and \(j=0,1,2,\ldots,L\), where \(\lceil \log_2 M \rceil = L\), we precompute the values \(\text{r}[i][j]\) defined by:

• Starting from Bus \(i\), Takahashi has transferred exactly \(2^j\) times (i.e. he is on his \((2^j+1)\)-th bus). Then, Takahashi is on the Bus \(\text{r}[i][j]\) (among the \(M\) buses).

Such values \(\text{r}[i][j]\) can be derived from the values \(\text{r}[1][j-1],\text{r}[2][j-1],\ldots,\text{r}[M][j-1]\) if \(1 \leq j\), so if we first compute \(\text{r}[i][0]\) for every \(i\ (1 \leq i \leq M)\), and then compute them in the increasing order of \(j\), we can find all the values in a total of \(O(ML)=O(M \log M)\) time.

Next, we will actually tackle the subproblem. In order to obtain “the index of the last bus Takahashi took when he has taken exactly \(k\) buses starting from Bus \(j\),” we update a variable \(now\) initialized with \(j\) with the following steps for each \(i=0,1,2,\ldots,L\):

• If the \(2^i\)’s place of \(k-1\) in binary notation is \(1\), then replace \(now\) with \(\text{r}[now][i]\).

The final value of \(now\) obtained in the steps above is what we wanted. The computation complexity is \(O(L)=O(\log M)\).

Therefore, we could solve the problem. The total time complexity is \(O(M \log M + Q \log^2 M)\).

Sample code (Python)

N,M,Q = map(int,input().split())
bus = [[] for i in range(N+1)]
A = [0]*M
B = [0]*M
S = [0]*M
T = [0]*M
for i in range(M):
    A[i],B[i],S[i],T[i] = map(int,input().split())
    bus[A[i]].append([S[i],i])
for i in range(N+1):
    bus[i].sort()
import bisect
dp = [[0]*20 for i in range(M)]
for i in range(M):
    nex = bisect.bisect_left(bus[B[i]],[T[i],0])
    if nex == len(bus[B[i]]):
        dp[i][0] = i
    else:
        dp[i][0] = bus[B[i]][nex][1]
for i in range(1,20):
    for j in range(M):
        dp[j][i] = dp[dp[j][i-1]][i-1]
for i in range(Q):
    X,Y,Z = map(int,input().split())
    now = bisect.bisect_left(bus[Y],[X,0])
    if now == len(bus[Y]) or Z <= bus[Y][now][0]:
        print(Y)
        continue
    now = bus[Y][now][1]
    if Z <= T[now]:
        print(*[A[now],B[now]])
        continue
    left = 0
    right = M
    while left+1 < right:
        mid = (left+right)//2
        now2 = now
        val = mid
        idx = 0
        while val:
            if val % 2:
                now2 = dp[now2][idx]
            val //= 2
            idx += 1
        if T[now2] < Z:
            left = mid
        else:
            right = mid
    val = left
    idx = 0
    while val:
        if val % 2:
            now = dp[now][idx]
        val //= 2
        idx += 1
    now2 = bisect.bisect_left(bus[B[now]],[T[now],0])
    if now2 == len(bus[B[now]]) or Z <= bus[B[now]][now2][0]:
        print(B[now])
        continue
    now = bus[B[now]][now2][1]
    print(*[A[now],B[now]])

Note that we can do the binary search during doubling. Sometimes such optimization is required, so it is worth remembering.

Sample code (C++)

#include<bits/stdc++.h>
using namespace std;

int main(){
	int N,M,Q; cin >> N >> M >> Q;
	vector<int> A(M),B(M),S(M),T(M);
	vector<vector<pair<int,int>>> bus(N+1);
	for(int i=0; i<M; i++){
		cin >> A[i] >> B[i] >> S[i] >> T[i];
		bus[A[i]].push_back(make_pair(S[i],i));
	}
	for(auto &p: bus) sort(p.begin(),p.end());
	vector<vector<int>> table(20,vector<int>(M));
	for(int i=0; i<M; i++){
		auto itr = lower_bound(bus[B[i]].begin(),bus[B[i]].end(),make_pair(T[i],-1));
		if(itr == bus[B[i]].end()) table[0][i] = i;
		else table[0][i] = (*itr).second;
	}
	for(int i=1; i<20; i++){
		for(int j=0; j<M; j++){
			table[i][j] = table[i-1][table[i-1][j]];
		}
	}
	while(Q--){
		int X,Y,Z; cin >> X >> Y >> Z;
		auto itr = lower_bound(bus[Y].begin(),bus[Y].end(),make_pair(X,-1));
		if(itr == bus[Y].end()){
			cout << Y << endl;
			continue;
		}
		int now = (*itr).second;
		if(Z <= S[now]){
			cout << Y << endl;
			continue;
		}
		if(Z <= T[now]){
			cout << A[now] << " " << B[now] << endl;
			continue;
		}
		for(int i=19; i>=0; i--){
			int next = table[i][now];
			if(T[next] < Z) now = next;
		}
		itr = lower_bound(bus[B[now]].begin(),bus[B[now]].end(),make_pair(T[now],-1));
		if(itr == bus[B[now]].end()){
			cout << B[now] << endl;
			continue;
		}
		int next = (*itr).second;
		if(Z <= S[next]){
			cout << B[now] << endl;
			continue;
		}
		cout << A[next] << " " << B[next] << endl;
	}
}

Otherwise, one can solve the problem by making use of the property that the following graph will form a rooted tree: each vertex represents a bus and each edge from a bus to another bus means that Takahashi will take the latter bus right after taking the former bus. The time complexity is again \(O((M+Q) \log M)\).

Sample code (C++)

#include<bits/stdc++.h>
using namespace std;

int main(){
    int N,M,Q; cin >> N >> M >> Q;
    vector<int> A(M),B(M),S(M),T(M);
    vector<vector<pair<int,int>>> bus(N+1);
    for(int i=0; i<M; i++){
        cin >> A[i] >> B[i] >> S[i] >> T[i];
        bus[A[i]].push_back(make_pair(S[i],i));
    }
    for(auto &p: bus) sort(p.begin(),p.end());
    vector<vector<int>> edge(M);
    vector<bool> isroot(M);
    for(int i=0; i<M; i++){
        auto itr = lower_bound(bus[B[i]].begin(),bus[B[i]].end(),make_pair(T[i],-1));
        if(itr == bus[B[i]].end()) isroot[i] = true;
        else edge[(*itr).second].push_back(i);
    }
    vector<int> X(Q),Y(Q),Z(Q);
    vector<vector<int>> ans(Q),mem(M);
    for(int i=0; i<Q; i++){
        cin >> X[i] >> Y[i] >> Z[i];
        auto itr = lower_bound(bus[Y[i]].begin(),bus[Y[i]].end(),make_pair(X[i],-1));
        if(itr == bus[Y[i]].end() || Z[i] <= (*itr).first) ans[i] = {Y[i]};
        else mem[(*itr).second].push_back(i);
    }
    vector<int> time,city;
    function<void(int)> dfs = [&](int now){
        time.push_back(T[now]);
        time.push_back(S[now]);
        city.push_back(B[now]);
        city.push_back(A[now]);
        for(auto i: mem[now]){
            int idx = lower_bound(time.rbegin(),time.rend(),Z[i])-time.rbegin();
            if(idx == time.size()) ans[i] = {*city.begin()};
            else if(idx%2 == 0) ans[i] = {*(city.rbegin()+idx)};
            else ans[i] = {*(city.rbegin()+idx-1),*(city.rbegin()+idx)};
        }
        for(auto next: edge[now]){
            dfs(next);
        }
        time.pop_back();
        time.pop_back();
        city.pop_back();
        city.pop_back();
    };
    for(int i=0; i<M; i++){
        if(!isroot[i]) continue;
        dfs(i);
    }
    for(auto p: ans){
        for(auto q: p) cout << q << " ";
        cout << endl;
    }
}