From the algorithm you describe and the comment you made on your question I understand that you are interested in walks and not in simple paths and that the graph is unweighted.
The algorithm you describe has at least exponential complexity since you examine each possible subset of $A$ (it has $2^{|A|} = O(2^N)$ subsets). We can use a much faster polynomial time algorithm.
One such algorithm has $O(k(N + M))$ computational complexity. You will need $k + 1$ boolean matrices, $B_0, B_1, \dots, B_k$, of size $N$. $B_0$ will be initialized to false, excluding its $s$-th position which will be initiallized to true. That is:
$B_0[u] = 1$, iff $u = s$
You will execute $k$ steps. In the $i$-th step you will calculate $B_i$ using the following formula:
$B_i[u] = 1$, iff $\exists v: (v, u) \in E \wedge B_{i - 1}[v] = 1$
This way $B_i[u]$ will be true iff there is a walk of length $i$ from $s$ to $u$. When you have calculated $B_k$ you can find the answer by checking $B_k$'s positions that match to $A$'s nodes' ids. Iff at least one of these array positions make $B_k$ true, then the answer is positive. The computational complexity of this algorithm is $O(k(N + M))$ and its space complexity is $O(kN)$. You can drop the space complexity to $O(N + M)$ if you observe that only the previous and the current boolean matrix $B_{i - 1}, B_i, \forall i \in \{1, 2, \dots, k\}$ are needed during each step. Bellow you can see a C++ implementation of this algorithm.
#include <iostream>
#include <vector>
using namespace std;
#define MAXN 100005 // Equal to maximum possible N + 5
vector<int> E[MAXN]; // Adjacency lists
vector<int> A; // Nodes that belong in set A
bool PB[MAXN]; // boolean array of the previous step
bool B[MAXN]; // current boolean array
int main()
{
int N, M, S_A, s, k;
// N: number of nodes
// M: number of directed edges
// S_A: size of set A
// s: starting node
// k: length of walk to be done
cin >> N >> M >> S_A >> s >> k;
for (int i = 0; i < M; i++)
{
int u, v;
cin >> u >> v;
E[u].push_back(v); // there is an edge from u to v
}
for (int i = 0; i < S_A; i++)
{
int a;
cin >> a;
A.push_back(a); // a belongs to A
}
// Create B_0
for (int i = 0; i < N; i++) B[i] = false;
B[s] = true;
// Perform the k repetitions
for (int i = 1; i <= k; i++)
{
// Copy B (B_i) to PB (B_{i-1}) and initialize B
for (int j = 0; j < N; j++)
{
PB[j] = B[j];
B[j] = false;
}
for (int j = 0; j < N; j++) if (PB[j])
{
// There is a walk from s to j of length i-1...
for (int l = 0; l < E[j].size(); l++)
{
// ...and an edge from j to E[j][l]...
B[E[j][l]] = true;
// ...so there is a walk of length i from s to E[j][l]
}
}
}
// Now B = B_k
for (int i = 0; i < S_A; i++)
{
int a = A[i];
if (B[a])
{
// There is a walk of length k from s to a
cout << "YES\n";
return(0);
}
}
// There is no walk of length k from k to any node of A
cout << "NO\n";
return(0);
}
This algorithm is fast if $k$ is small but when $k$ becomes larger its computational complexity increases significantly. We will now describe an algorithm which is faster for large values of $k$.
You can use the fact that the $(i, j)$ cell of $Adj^p$ (where $Adj$ is the adjacency matrix of the graph) is equal to the number of walks from $i$ to $j$ in the graph with length equal to $p$ (easy to prove using mathematical induction). In your case you are interested in walks of length equal to $k$ so you have to calculate $Adj^k$. You can efficiently calculate $Adj^k$ in $\lceil log_2k \rceil$ steps using exponentiation by squaring. Each step of the above calculation requires squaring an $N \times N$ matrix. This can be done in $O(N^3)$ using the naive algorithm or in $O(N^{2.373})$ using the fastest known algorithm. The fact that you want to know if a path between two nodes exists and not the number of different paths allows you to avoid large numbers by making the integer matrix boolean after each multiplication. When you have $Adj^k$ (in boolean format) you can check if there is a path of length $k$ from $s$ to a vertex in $A$ by checking the $s$-th row of $Adj^k$ and specifically the columns regarding nodes of $A$. Formally the result is positive iff
$$\exists u \in A: Adj^k(s, u) = 1$$
The algorithm is the following:
- Read input and calculate adjacency matrix of the graph: $O(N^2)$
- Calculate $Adj^k$ using exponentiation by squaring, fast matrix multiplication and making matrix boolean after each multiplication: $O(N^{2.373}logk)$
- Check if there exists a cell $Adj^k(s, u) = 1$ where $u \in A$: $O(N)$
- Output the reslt: $O(1)$
The total computational complexity of this algorithm is $O(N^{2.373}logk)$ and it is better than the first algorithm for large values of $k$. Both of them also include possibility of self edges as you mention in your comment on your question.
You are able to combine the two algorithms and run the fastest one based on the values of the parameters and achive $O(\min \{k(N + M), N^{2.373}logk\})$ computational complexity.
