Counting adjacency matrices Here is a question that has come up in the context of a problem that involves counting partially ordered sets.  
For an adjacency matrix $A$, let $p$ be the sum of elements in the strict upper triangle (upper triangle minus the diagonal) of $A$, and $q$ be the sum of elements in the strict upper triangle of $A^2$. For fixed values of $p$ and $q$, is it possible to compute the cardinality of the set of all such matrices $A$? 
If yes, how does one go about it? My guess is that the problem may be easier to tackle if we demand that $A$ has some extra symmetry, but I have not been able to arrive at any definite conclusion. Though I am interested in the generic case (without any added symmetries), solutions for any special cases will also be helpful. So will be any references that deal with similar problems.
 A: EDIT (4/21/20): new link re: function inserted at end

Source: my tweets, with minor errors removed (https://twitter.com/krzhang/status/1252529588049072128)
Let's assume $A$ is symmetric and 0 on the diagonal. (disclaimer: I'm guessing this is not what you meant by "symmetry" because of poset context, but it may still be helpful) This means we are really working with unlooped undirected graphs, where
- p is how many edges the graph has, and 
- q is now many 2-paths (disregarding order) that are not loops. 
Now, cool observation: 2-paths that are not loops can be identified with their middle point and 2 neighbors. So for each vertex $i$ of deg. $d_i$, it contributes $d_i(d_i-1)/2$ 2-paths that are not loops.
So our problem becomes: "How many ways are there to split $p$ into nonnegative integers $d_1 + ... + d_n,$ such that $\sum d_i(d_i - 1)/2 = q$?"
Some manipulation gives $\sum d_i^2 = 2q + p$, so this problem really reduces to
"Given the first and 2nd power sums of $d_1 ... d_n$, how many sets of nonnegative $d_n$ are there?" or the quite-beautiful probabilistic form:
"How many nonnegative integral distributions are there of a fixed mean and variance?"
There're number theory constraints here, so I guess this is hard (which means original problem is even harder). However, computationally this isn't bad. Here's a solution:


*

*construct a 3-d infinite array so that $P[x][y][z]$ = "the number of ways to solve this with $x$ numbers such that their 1-power sum (sum) is $y$ and their 2-power sum is $z$

*use dynamic programming to build this layer by layer by $x$. So compute everything with $x = 1$ first, then reduce each problem with $x+1$ to those with $x$ by summing over different values for the first element.


This gives an $O(n^2p(2q+p))$ algorithm.
Link (h/t Boris Alexeev): This last function is explored at https://mathworld.wolfram.com/SumofSquaresFunction.html as "sum of squares function." As I predicted, it seems number-theoretical and thus needs to be written as sums of modular functions for larger $n$ and $k$. Therefore, finding a closed form seems super hard.
