Counting Closed Walks in a Directed Bridgeless Graph with Constraints Short version: Given a directed bridgeless (possibly non-simple) graph $G$ and a positive integer $n$, let $W(G,n)$ be the number of "walks covers" of $G$ of length $n$, where a walk cover is defined as a walk through $G$ that starts and ends on the same vertex and passes through every vertex and every edge at least once. Walk covers that are rotations of each other are not counted as distinct.
For example, for the following graph

we have the following values:
$$\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
n & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15\\
\hline
W(G,n) & 0 & 0 & 0 & 1 & 1 & 1 & 2 & 3 & 3 & 5 & 6 & 9 & 11 & 16 & 21 \\\hline
\end{array}$$
Is there a way to compute $W(G,n)$ without actually finding all the walks with these constraints, possibly using the adjacency matrix?
Long version: Apologies in advance for any misuse of terminology. My background is in physics, and my suspicion is that my question has a simple answer using some basic techniques from graph theory, but my cursory search has not yielded anything.
I've been investigating juggling patterns, which can be related to closed walks in a directed bridgeless non-simple graph known as the juggling state graph, where the vertices represent various arrangements of the balls in the air, and the edges represent throws that connect these states. This graph is "minimally" non-simple, as it only has one self-loop, connecting the so-called "ground state" to itself. For more on the juggling state graph, see here.
I'm specifically interested in counting the composite juggling patterns, which are patterns that can be "factored" into smaller patterns. In terms of the graph, this means that the closed walk representing the pattern can be written as two or more shorter closed walks in sequence, where all the sub-walks start and end on the same vertex. Looking at the vertices and edges in these closed walks, we find various subgraphs. For example, here are the possible subgraphs that can represent composite juggling patterns of length $n=4$:

And here are the subgraphs for composite juggling patterns of length $n=5$:

Note that the number of edges need not equal $n$, as repeating edges is allowed. Once these possible patterns are found, it is computationally quite easy to search for these subgraphs within the larger juggling state graph. We can count the number of times each subgraph $g$ appears in the full juggling state graph $G$, i.e., the number of subisomorphisms of between $g$ and $G$, dividing by the order of the automorphism group of $g$ to deal with overcounting. Let us denote this number of subisomorphisms with the overcounting corrected as $S(g,G)$.
If we simply add up $S(g,G)$ for all possible $g$, we will not get the number of composite patterns. This is because the correspondence between a juggling pattern and a particular $g$ is not one-to-one. For an example, consider the following graph:

This is one of the possible subgraphs for $n = 3$, and in that case it can only correspond to one kind of juggling pattern, namely a 1-cycle followed by a 2-cycle, or, in terms of the edges, $(a,b,c)$, up to rotation of the pattern. However, this is also a possible graph for $n = 5$, and in that case it can represent two different kinds of juggling patterns: a 1-cycle followed by two 2-cycles; $(a,b,c,b,c)$, or three 1-cycles followed by one 2-cycle: $(a,a,a,b,c)$. These are distinct walks, and thus distinct juggling patterns, corresponding to the same subgraph. Using the notation above, we would write $W(g,3) = 1$ and $W(g,5)=2$.
For another example with a simple graph this time, consider the following:

For $n=4$, there is only one pattern: $(b,a,c,d)$, but for $n=6$, there are two: $(b,a,c,d,c,d)$ and $(b,a,b,a,c,d)$. So here we have $W(g,4) = 1$ and $W(g,6)=2$.
Thus the correct expression for the number of composite patterns of length $n$ is
$$N_\text{composite}(G,n) = \sum_g W(g,n) S(g,G),$$
where the sum is done over the possible subgraphs. To finally arrive at the question, I would like to be able to compute $W(g,n)$ without actually figuring out all the walks explicitly. Right now I have some code that works out all the walk covers and then counts them, but as $n$ increases, it becomes very slow for certain graphs.
One idea I had was to use the trace of the $n$th power of the adjacency matrix, since the $(i,j)$ entry of $A^n$ counts the number of walks from vertex $i$ to vertex $j$ of length $n$. However, this count includes walks that don't cover the whole graph. There might be a clever way to subtract those off to get $W(g,n)$, but I couldn't work it out.
Postscript: I would also like a better understanding of how $W(g,n)$ depends on various properties of $g$. Here's a table of $W(g,n)$ for some small graphs:

Here are some of the patterns I've observed:

*

*Obviously, $W(g,n)$ is zero if $n$ is less than the number of edges in the graph. However, the converse is not true.

*If $g$ is a cycle of length $m$, then $W(g,n) = 1$ for $n \equiv 0 \text{ mod } m$ and $W(g,n) = 0$ otherwise.

*If $g$ is non-simple, $W(g,n)$ is never zero once $n$ is greater than or equal to the number of edges.

*If $g$ consists of two cycles of distinct lengths $m_1$ and $m_2$, then $W(g,n)$ is the number of compositions of $n$ (partitions where order matters) using only $m_1$ and $m_2$. This seems to be true whether the two cycles share one or more vertices.

*If $g$ consists of two distinct cycles of the same length $m$, then $W(g,n) = 0$ for $n \not\equiv 0 \text{ mod } m$ or $n < 2m$ and for $n \equiv 0 \text{ mod } m$ and $n \geq 2m$, $W(g,n)$ matches OEIS A052823, which has to do with necklace polynomials with two colors. This likely generalizes to more than two cycles of the same length, but I haven't checked.

*Adding a self-loop to a graph removes all zeros for sufficiently large $n$, but the placement of the loop matters greatly. For example, look at the last four columns of the table above. Placing the self-loop on particular vertices can cause the number of walk covers to increase dramatically. (Computing the bottom right value of 2004 took my laptop 4 hours.)

 A: Let our directed graph be $G=(V,E)$ where $V$ is the set of vertices and $E$ is the set of directed edges. Assume it has $n$ vertices and $m$ directed edges.
We first solve the problem of calculating $C(G,k)$ where we define $C(G,k)$ as the number of different types of closed walks in $G$ under the equivalence relation where two closed walks are the same if they are cyclic permutations  of each other.
We can calculate this with mobius inversion but I am not completely sure how to explain it, the idea is that you can find how many closed walks are the same when they are rotated $d$ times. This number is equal to the number of closed walks of length $d$ and can be calculated with the adjacency matrix of $G$.
Once we are able to calculate $C(G,k)$ we can calculate the number that we want, which is $W(G,k)$ by doing inversion again, but this time over the lattice of subsets of $E$.
For each subset $E'$ of $E$ we can consider the spanning subgraph $G'$ of $G$ with those edges. Let $C(G',k)$ be the number of closed walks in that subgraph.
We can calculate $W(G,k)$ as
$$\sum\limits_{ E'\subseteq E} C(G',k) (-1)^{|E| - \left|E'\right|}$$
Here is some code with c++ and armadillo
#include <bits/stdc++.h>
#include <armadillo>

using namespace arma;
using namespace std;

const int MAX = 110;
int n,m,k;
int e[MAX][2]; //stores the endpoints of each edge
int a[MAX]; // bitmask for the subset of edges
int d; // the number of divisors of k
mat A; // the adjacency matrix of the graphs
mat D; // the weird matrix for the mobius inversion
vector <int> divs;

int next(){ // iterates through the subset of edges
        for(int i=m-1;i>=0;i--){
                if( a[i] == 0){
                        a[i] = 1;
                        for(int j=i+1;j<m;j++){
                                a[j]=0;
                        }
                        return 1;
                }
        }
        return 0;
}

void buildadj(){ // builds the adjacency matrix for current bitmask 
        for(int i=0;i<n;i++){
                for(int j=0;j<n;j++){
                        A(i,j)  = 0;
                }
        }
        for(int i=0;i<m;i++){
                if( a[i] ) A(e[i][0],e[i][1]) = 1;
        }
}

long long walksperiod(int d){ // gets the number of d-periodic walks for current bitmask
        mat B = powmat(A,d);
        for(int i=0;i<n;i++){
                for(int j=0;j<n;j++){
                        if( i != j) B(i,j) = 0;
                }
        }
        long long res = 0;
        for(int i=0;i<n;i++){
                res += B(i,i);
        }
        return res;
}

long long solve(){ // gets the value for the current bitmask
        // what we want to do here is find the number of walks of each period divisible by d
        // we do that with mobius inversion
        long long res = 0;
        mat B(d,1);
        for(int i=0;i<d;i++){
                B(i,0) = walksperiod(divs[i]);
        }
        B = inv(D)*B;
        for(int i=0;i<d;i++){
                cout << divs[i] << " " << B(i) << endl;
                res += B(i)/divs[i];
        }
        return res;
}

int main(){
        cin >> n >> m >> k; // insert n and m and k
        A = mat(n,n);
        for(int i=1;i<=k;i++){ // we get the divisors of k for the moebius inv
                if( k%i == 0){
                        divs.push_back(i);
                }
        }
        d = divs.size();
        D = mat(d,d);
        for(int i=0;i<d;i++){
                for(int j=0;j<d;j++){
                        if( divs[i]%divs[j] == 0) D(i,j) = 1;
                        else D(i,j) = 0;
                }
        }

        for(int i=0;i<m;i++){
                int a,b;
                cin >> a >> b;
                e[i][0] = a;
                e[i][1] = b;
        }
        int start = 1;
        long long res = 0;
        while( start || next() ){
                start = 0;
                buildadj();
                int mp = 0; //current number of edges
                for(int i=0;i<m;i++){
                        if( a[i] ) mp ++;
                }
                if( (mp%2) == (m%2) ) res += solve();
                else res -= solve();
        }
        cout << res << endl;

}

A: To put a bow on Onir's excellent answer, we can derive an expression for $C(g,k)$ in terms of the adjacency matrix of $g$. We first define three different functions, suppressing the $g$ argument for now:

*

*$C_\text{total}(k)$ is the number of closed walks of length $k$, allowing for both cyclic permutations and walks that just repeat a subwalk of length $d$ $k/d$ times, where $d \mid k$

*$C_\text{proper}(k)$ is the number of closed walks not counting cyclic permutations as distinct and also ignoring walks with repeating subwalks.

*$C_\text{unique}(k)$ is the number of closed walks not counting cyclic permutations as distinct but allowing for repeating subwalks. (This is what we previously called $C(g,k)$.)

We see that $C_\text{total}(k) = \text{tr}\left(A^k\right)$, where $A$ is the adjacency matrix of $g$. We see that the walks counted by $C_\text{total}(k)$ can be grouped by their "true" length, i.e., the length before they repeat, which must divide $k$. Each walk of true length $d$ is counted by $C_\text{total}(k)$ $d$ times since we are allowing cyclic permutations, so we find
$$C_\text{total}(k) = \sum_{d\mid k}d \cdot C_\text{proper}(d).$$
Möbius inversion then yields
$$C_\text{proper}(k) = \frac{1}{k} \sum_{d\mid k} 
\mu\left(\frac{k}{d}\right) C_\text{total}(d)
= \frac{1}{k} \sum_{d\mid k} \mu\left(\frac{k}{d}\right) 
\text{tr}\left(A^d\right),$$
where $\mu(n)$ is the Möbius mu function. We can also show
$$C_\text{unique}(k) = \sum_{d\mid k} C_\text{proper}(d)
= \sum_{d\mid k}
\frac{1}{d} \sum_{c\mid d} \mu\left(\frac{d}{c}\right) 
\text{tr}\left(A^c\right)
= \frac{1}{k}\sum_{d \mid k}\phi\left(\frac{k}{d}\right)\text{tr}\left(A^d\right),
$$
where $\phi(n)$ is the Euler totient function. The key property in proving the last equality is
$$\phi(k) = k \sum_{d\mid k} \frac{\mu(d)}{d}.$$
