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 pIt took 21.69 sec for afind(1,940). I tried Cython, but it didn't help. How can I speed this up?
Modular arithmetic brute-force congruence finder
python;performance;python 3.x;mathematics
Quick manual optimization can be made by avoiding large intermediates. This is possible by generating sequential cubes with two consecutive additions:cubes:  0     1     8    27    64   125diff:      1     7    19    37    61diff:         6    12    18    24This allows good use of % to keep the numbers small. Eg:def afind(h, v):    P = 23857201    x_cubed_plus_v_delta = 1    x_cubed_plus_v = v    for x in range(P):        x_cubed_plus_v %= P        if x_cubed_plus_v == 0:            if (3*x**2 + 3*h*x + h**2) * h % P == 0:                yield x        x_cubed_plus_v_delta += 6 * x        x_cubed_plus_v += x_cubed_plus_v_deltaprint(list(afind(1, 940)))I used yield for simplicity; it doesn't noticeably affect timings. Using P as a variable has a small but unimportant cost.Times:Interpreter  Before  After  Before/AfterCPython 2     32.6    3.5       9.3CPython 3     12.1    6.3       1.9PyPy 2        19.8    0.4      49.5 PyPy 3        19.5    0.4      48.8Unfortunately CPython 3 went from fastest to slowest :(.A Cython version would look like:from libc.stdint cimport uint64_tdef afind(h, v):    cdef uint64_t P, x_cubed_plus_v_delta, x_cubed_plus_v, x    P = 23857201    x_cubed_plus_v_delta = 1    x_cubed_plus_v = v    res = []    for x in range(P):        x_cubed_plus_v %= P        if x_cubed_plus_v == 0:            if (3*x**2 + 3*h*x + h**2) * h % P == 0:                res.append(x)        x_cubed_plus_v_delta += 6 * x        x_cubed_plus_v += x_cubed_plus_v_delta    return resYou could quickly use it with:import pyximportpyximport.install()import pxprint(px.afind(1, 940))(not the recommended method for distibution).This is fast, taking ~0.075s. It doesn't support PyPy. Note that if P is determined at runtime, it takes closer to ~0.27s, which means PyPy is effectively close to optimal.Just for fun, here's the absolute fastest I could get this to go:from libc.stdint cimport uint64_tdef afind(uint64_t h, v):    cdef uint64_t P, P3, x_cubed_plus_v_delta, x_cubed_plus_v, x    P = 23857201    P3 = P * 3    x_cubed_plus_v_delta = 1    x_cubed_plus_v = v    res = []    for x in range(P):        if x_cubed_plus_v >= P:            x_cubed_plus_v -= P        if x_cubed_plus_v == 0:            if (3*x**2 + 3*h*x + h**2) * h % P == 0:                res.append(x)        x_cubed_plus_v_delta += 6 * x        if x_cubed_plus_v_delta >= P3:            x_cubed_plus_v_delta -= P3        if x_cubed_plus_v_delta >= P:            x_cubed_plus_v_delta -= P        if x_cubed_plus_v_delta >= P:            x_cubed_plus_v_delta -= P        x_cubed_plus_v += x_cubed_plus_v_delta    return resThis just avoids modulo, replacing it with subtraction. To get significantly faster you'll do best looking for a clever algorithm.