Wrong calculation when using Cython

Question

I implemented Lucas-Lehmer primality test to check Mersenne prime in python. Then I use Cython to speed up the calculation.

Original Python code:

def lucas_lehmer(p):    
    if p == 2:
        return True

    s = 4
    M = (1 << p) - 1

    for i in range(p-2):    
        s = ((s * s) - 2) % M
        print("Processed: {}%".format(100*i//(p-2)))

    if s == 0:
        return True

    else:
        return False

Cython code:

cpdef lucas_lehmer(int p):    
    if p == 2:
        return True

    cdef unsigned long long int M
    M = (1 << p) - 1

    cdef unsigned long long int s
    s = 4

    cdef int i

    for i in range(p-2):    
        s = ((s * s) - 2) % M
        print("Processed: {}%".format(100*i//(p-2)))

    if s == 0:
        return True

    else:
        return False

Running the original Python code, it works correctly. But for Cython, it's only correct with p = 31 and lower, testing with p = 61 and bigger (all tested p values are values that 2^p-1 is prime), it returns False (not a prime number), except for p = 86243.

For some p like 97, even though 2^97-1 is not a prime number, the program actually return True (is a prime number), which is a contradiction.

Why does this happen? Without using cdef for variable M and s, the calculation will be correct, but the performance won't get any improved.

Why is there a `M` parameter in the Python code if it's overwritten immediately after? — Guimoute, Mar 22 '20 at 16:48
Surely that's just some overflow? Are you asking where exactly it happens? — Kelly Bundy, Mar 22 '20 at 16:55
`M = (1 << p) - 1` How do you think `M` is being stored? Perhaps you should try to define a type that Cython and C understands. There are examples on StackOverflow and elsewhere that show how to define 128-bit integers for Cython. — rickhg12hs, Mar 22 '20 at 17:11
@rickhg12hs I've tried from the following link: https://stackoverflow.com/questions/27582001/how-to-use-128-bit-integers-in-cython. When I compile, I get the error: `int128: undeclared identifier` — The 2nd, Mar 22 '20 at 17:32
128bit integer is only a gcc extension (do you use gcc?), and it would not help here because the problem happens when (1<<32) is executed. — ead, Mar 22 '20 at 19:15
However, for p>64 (or p>128) you will need to use libraries which support arbitrary size integer (or use Python-integer, which isn't a bad choice). — ead, Mar 22 '20 at 19:16

Alistair · Answer 1 · 2020-03-22T20:20:18.430

Running a few tests on your code I found that M was always equal to 1 so I defined p as a cdef and got the required result.

Not sure exactly what the issue is but it's something to do with that bit operation on p. p needs to be of the same type as M for it to make sense and if one is cdef and one is python int somehow it doesn't work?

cpdef lucas_lehmer(int py):    
    cdef p
    p = py
    if p == 2:
        return True

    cdef M
    M = (1 << p) - 1

    cdef s
    s = 4

    cdef int i

    for i in range(p-2):    
        s = ((s * s) - 2) % M
        print("Processed: {}%".format(100*i//(p-2)))

    if s == 0:
        return True

    else:
        return False

This code works fine for me but I don't see any performance improvement. It didn't faster than python interpreter. — The 2nd, Mar 24 '20 at 03:39

Wrong calculation when using Cython

1 Answers1