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

Answer 1

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).