Find a multiplier which multiplies given indices into a formation of bits

I have a set of indices and a set of bit masks (the size of the sets is the same, and between 6-4096 but much closer to the lower bound). I also have a single bit mask which I will call the superset. I need to find a factor, that if for any index I have, if I multiply the factor and the index and mask bitand with the superset, the result is the subset. The masks and factor will all be 64-bit wide. This example is 8-bit cause I can't write 64-bit masks here obviously

mask = 00111000         ((p * index) & mask) == subset
index=3 sub=00010000 => ((p * 3) & mask) == 00010000
index=5 sub=00100000 => ((p * 5) & mask) == 00100000

Really, my input is just a set {i0, i1, ... iN} of indices, {s0, s1, ... sN} of masks, and a superset mask.
Of course, all the subsets are subsets of the super-set. For now, for different indices, there can be the same subset results, but I will be happy to receive soloutions that act as if this was not a constraint.

My best try for now is just guessing random numbers for the factor, and checking to see for each indice-subset pair if it holds true. If not, I go to the next number. But for the smallest form of this problem, I had 13 such pairs, but since the brute force algorithm is O(2^n) with a very large constant factor, I can't find a factor which holds true even for 7 pairs!

What I thought about is looking at the resulting values from the multiplication. Let's look at this problem from an 8-bit form: (Here a,b,c,d, etc mean if the i'th bit is set in the result of the multiplication)

P * Index = 2^0*a + 2^1*b + 2^2*c + 2^3*d + 2^4*e + 2^5*f + 2^6*g + 2^7*h
P = (1a + 2b + 4c + 8d + 16e + 32f + 64g + 128h) / Index
Now we want P to be integer, so do modulo (let's say Index is 5 for example)
P = (1a + 2b + 4c + 3d + 1e + 2f + 4g + 3h) / 5
P = (1(a+e) + 2(b+f) + 4(c+g) + 3(d+h)) / 5

So now we can say that every bit-mask that its sum of set bits multiplied by the coefficients above modulo 5 is 0, it's a valid P * Index, then we can find P by a simple divison. Hope you get me.

But the problem is I don't know how to extend this to multiple index-result pairs, and I really don't know what to do now.

You can use an SMT solver like z3 for this. You can create variables for all the things you want calculated and then set constraints on them. Here I'm assuming that the values you can choose are p and the indices, that is still not clear from your question. Completely untested example:

void solve(const std::vector<uint64_t>& subs, uint64_t mask) {
    using namespace z3;
    context c;
    solver s(c);
    expr p = c.bv_const("p", 64);
    std::vector<expr> indices;
    for (int i = 0; i < subs.size(); ++i)
        indices.push_back(c.bv_const(("indices_" + std::to_string(i)).c_str(), 9));

    for (int i = 0; i < subs.size(); ++i)
        s.add((p * zext(indices[i], 64 - 9) & mask) == subs[i]);

    if (s.check() == sat) {
        model m = s.get_model();
        uint64_t p_result = m.eval(p).get_numeral_uint64();
        std::vector<uint16_t> indices_result;
        for (int i = 0; i < subs.size(); ++i)
           indices_result.push_back(m.eval(indices[i]).get_numeral_uint64());
        // do something with the result
    } else {
        // no solution
    }
}