If you are familiar with solving regular systems of equations, this is not a major step up. When using real numbers in a system of equations, you do elimination like this:
[a b; c d] -> [a b; 0 d-(b*c/a)] -> [a 0; 0 d-(b*c/a)] -> [1 0; 0 1]
Note: here I use MATLAB matrix notation for ease-of-entry.
The important realization to make is that all of these matrix operations (i.e. dividing, multiplying, adding, and subtracting) exist for any field, not just the real numbers. If you are unfamiliar with the term field, it just means a group of values which allow for multiplication, negation, inversion, addition, etc.
This takes us to solving systems of xor equations. You currently describe your systems as a bunch of 16-bit values xor'd together. However, the way I'd chose to represent it as a bunch of bits xor'd together, For instance, if your first equation was:
p[0] = a[1] ^ a[2]
I'd represent this as:
p[0][0] = a[1][0] ^ a[2][0]
p[0][1] = a[1][1] ^ a[2][1]
…
Where the second set of brackets denotes the bit
offset in the 16-bit value. So, each of your little equations would be equivalent to 16 equations.
Single bit xor operations on booleans form a field. In this field, we make the "addition" operator equivalent to xor. We can define the addition and multiplication tables as follows:
1 + 0 = 0 + 1 = 1; 1 + 1 = 0 + 0 = 0
1 * 0 = 0 * 1 = 0 * 0 = 0; 1 * 1 = 1
Division is only possible by 1 (because you cannot divide by zero), and so the division operator leaves an element unchanged.
With this, you should be able to form a matrix for your system of xor equations. This matrix will consist entirely of 1's and 0's. Then use the gauss-jordan elimination algorithm (it's really not too hard to implement) just as you would for ordinary real numbers. This will allow you to invert the matrix and find the solution.
I was personally so intrigued by this question that I wrote up a small C++ matrix implementation that allows you to supply any field you like. This may be a good starting point, or you may even wish to use my code altogether! Here's the source code on Github: XorSystem. I specifically suggest looking at the invert() method on ANMatrix.