How to find the fixed points of a simple mod function elegantly?

Let's use Z3 solver:

(declare-const x (_ BitVec 32))
(assert (= x (bvxor (bvmul x #x00000156) #xfca802c7)))
(check-sat)
(get-model)

The result is '#x08f2cab1' = 150129329.

Python code:

def f(x, n):
    return ((x*0x156)^0xfca802c7) % n


solns = [1]  # The one solution modulo 2, see text for explanation
n = 1
while n < 2**32:
    prev_n = n
    n = n * 2
    lifted_solns = []
    for soln in solns:
        if f(soln, n) == soln:
            lifted_solns.append(soln)
        if f(soln + prev_n, n) == soln + prev_n:
            lifted_solns.append(soln + prev_n)
    solns = lifted_solns

for soln in solns:
    print soln, "evaluates to ", f(soln, 2**32)

Output: 150129329 evaluates to 150129329

Idea behind the algorithm: We are trying to find x XOR 0xfca802c7 = x*0x156 modulo n, where in our case n=2^32. I wrote it this way because the right side is a simple modular multiplication that behaves nicely with the left side.

The main property we are going to use is that a solution to x XOR 0xfca802c7 = x*0x156 modulo 2^(i+1) reduces to a solution to x XOR 0xfca802c7 = x*0x156 modulo 2^i. Another way of saying that is that a solution to x XOR 0xfca802c7 = x*0x156 modulo 2^i translates to one or two solutions modulo 2^(i+1): those possibilities are either x and/or x+2^i (if we want to be more precise, we are only looking at integers between 0, ..., modulus size - 1 when we say "solution").

We can easily solve this for i=1: x XOR 0xfca802c7 = x*0x156 modulo 2^1 is the same as x XOR 1 = x*0 mod 2, which means x=1 is the only solution. From there we know that only 1 and 3 are the possible solutions modulo 2^2 = 4. So we only have two to try. It turns out that only one works. That's our current solution modulo 4. We can then lift that solution to the possibilities modulo 8. And so on. Eventually we get all such solutions.

Remark1: This code finds all solutions. In this case, there is only one, but for more general parameters there may be more than one.

Remark2: the running time is O(max[number of solutions, modulus size in bits]), assuming I have not made an error. So it is fast unless there are many, many fixed points. In this case, there seems to only be one.