Modular arithmetic brute-force congruence finder

Question

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?

Answer 1

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   125
diff:      1     7    19    37    61
diff:         6    12    18    24

This 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_delta

print(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/After
CPython 2     32.6    3.5       9.3
CPython 3     12.1    6.3       1.9
PyPy 2        19.8    0.4      49.5 
PyPy 3        19.5    0.4      48.8

Unfortunately CPython 3 went from fastest to slowest :(.

A Cython version would look like:

from libc.stdint cimport uint64_t

def 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 res

You could quickly use it with:

import pyximport
pyximport.install()

import px

print(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_t

def 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 res

This just avoids modulo, replacing it with subtraction. To get significantly faster you'll do best looking for a clever algorithm.

Answer 2

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.