# 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?

• Math.SE question. Dec 17, 2014 at 3:09

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.

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.