Prueba Lucas Lehmer

Question

He intentado escribir una función que compruebe si es primo utilizando el Test de Lucas Lehmer

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

pero no funciona

Answer 1

Como señala gammatester, debes tomar el residuo menos positivo en cada iteración.

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

Anexo : Como señala gammatester, esto sólo es fiable para valores pequeños de $p$ (por ejemplo $p < 32$ ) ya que s * s puede desbordarse. Hay código aportado por los usuarios en el Intercambio de archivos MATLAB para aritmética de enteros de precisión arbitraria, pero si desea una implementación seria, probablemente sea mejor utilizar otra cosa (por ejemplo, C con GNU GMP ).