7
\$\begingroup\$

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?

\$\endgroup\$
  • \$\begingroup\$ Please fix your indentation. \$\endgroup\$ – 200_success Dec 17 '14 at 2:46
  • 1
    \$\begingroup\$ Math.SE question. \$\endgroup\$ – Brythan Dec 17 '14 at 3:09
2
\$\begingroup\$

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.

\$\endgroup\$
2
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.