Modular arithmetic brute-force congruence finder

My full code is running too slow. I do profile.run to my project and found out this function consumes a lot of time.

The main problem is just to return an array containing all \$x\$ satisfying both of these equations with some large primes like \$23857201\$.

$$(x^3+v) \mod P = 0$$

$$(3xx + 3hx + h^2) \mod P = 0$$

def afind(h,v):
        p=[]
        for x in range(23857201):               
           if (x**3+v) % 23857201 == 0 and (3*x*x+3*h*x+h**2)*h%23857201==0:
              p.append(x)           
     return p

It took 21.69 sec for afind(1,940). I tried Cython, but it didn't help. How can I speed this up?

In case it is you who figured out an \$h\$ factor in the code snippet, add \$x^3\$ to both sides and rewrite the second equation as \$(x + h)^3 \equiv x^3 (\mod P)\$, so that with the first one in mind, the system comes down to

\$x^3 + v \equiv 0 (\mod P)\$

\$(x + h)^3 + v \equiv 0 (\mod P)\$

The programmatic approach is obvious:

-
find all solutions to \$x^3 + v \equiv 0 (\mod P)\$, and

-
select those which are exactly \$h\$ apart.

Of course the first step is a key (the second is trivial). I am not that good in a theory of Diophantine cubics to suggest any more.

PS: Math.SE is indeed a right place to ask.

StackExchange Code Review Q#73884, answer score: 2