lucas.gp

lucas(p) =
{ my(u = 4, q = 1<<p - 1);
  for(k=3, p, u = (sqr(u)-2) % q);
  u == 0;
}