Lucas Lehmer Test

Question

I have tried to write a function that test if it is prime using Lucas Lehmer Test

function y=ex2(p)
s=4;
for i=3:p
    s=s.^2-2;

end
if (mod(s,2.^p-1)==0)
    y=1;
else
    y=0;
end
end

but it does not work

I guess your ex2 overflows. Apply mod to s=s.^2-2; and your final test is wrong, see e.g. https://en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test — gammatester, Jul 08 '18 at 13:51
You're using this for numbers of the form $2^n - 1$ or for integers in general? — Robert Soupe, Jul 08 '18 at 18:26

parsiad · Accepted Answer · 2018-07-08T15:29:30.793

As gammatester points out, you should take the least positive residue at each iteration.

function is_prime = ex2(p)
    assert(floor(p) == p);
    assert(p > 2);

    p = uint64(p);
    s = uint64(4);
    M = 2^p - 1;

    for i = 1:p-2
        s = mod(s * s - 2, M);
    end

    is_prime = (s == 0);
end

Addendum: As gammatester points out, this is only reliable for small values of $p$ (e.g., $p < 32$) since s * s may overflow. There is some user-contributed code on the MATLAB File Exchange for arbitrary precision integer arithmetic, but if you want a serious implementation, it's probably best to use something else altogether (e.g., C with GNU GMP).

(+1) But you should add, that this is reliable only up to $p<32$ because otherwise s*s might overflow, e.g. the function gives a wrong result for $p=61$ with Octave. — gammatester, Jul 08 '18 at 15:12

Lucas Lehmer Test

1 Answers1