methodA[n_Integer?Positive] :=
Module[
  { i, js, capD = 5 },
  If[Floor[Sqrt[n]]^2 == n, Return[{ 0, 0, 0 }] ];
  For[i = 1, i <= 1000, i++,
    js = JacobiSymbol[capD, n];
    If[js == 0,
      (* if capD == (+,-)n, (e.g., n = 5 or 11) try the next capD instead of quitting;
         this small modification of Selfridge's method A
         enables 5 and 11 to be classified as lucas probable primes.
      *)
      If[Abs[capD] != n, Return[{ 0, 0, 0 }] ]
    ];
    If[js == -1, Return[{ capD, 1, (1 - capD) / 4 }] ];
    If[capD < 0, capD = -capD + 2 , capD = -(capD + 2)]
  ];
  Print["n = ", n, ": could not find D, P, and Q"]
];

getBitsAndS[m_Integer?Positive] :=
Module[
(* express m as d * 2^s, where d is odd, and get the bits of m *)
  { d, s, bits },
  s = IntegerExponent[m, 2];
  d = m/(2^s);
  bits = IntegerDigits[m, 2];
  { d , s, bits }
];

lprpTestA[n_Integer?Positive] :=
Module[
  (* given n, this returns a list of two numbers:
       if n is a lucas probable prime, then the first number is 1; otherwise, it is 0.
       if n is a strong lucas probable prime, then the second number is 1; otherwise, it is 0.
       therefore, if n is an odd prime, this will return { 1, 1 }.
       if n < 3, is even, or is a square, then return { 0, 0 }.
     this is based on "Lucas Pseudoprimes" by Robert Baillie and S. S. Wagstaff, Jr., in
     Mathematics of Computation, volume 35, number 152 (October, 1980), pp. 1391-1417.
     this uses Selfridge's method A on page 1401 of that article to obtain D, P, and Q.
     note: a lucas pseudoprime is a lucas probable prime that is composite.
  *)
  { lprp = 0, strongLprp = 0,
    capD, capP, capQ,
    d, s, bits, nBits, i,
    uk, vk, qk, uNew, vNew, qPowerNew, uTemp, vTemp
  },

  If[ (n < 3) || EvenQ[n],
    Return[{ 0, 0 }]
  ];
  { capD, capP, capQ } = methodA[n];
  If[capD == 0, Return[{ 0, 0 }] ];

  (* d is the odd part of n + 1 *)
  { d, s, bits } = getBitsAndS[ n + 1 ];
  nBits = Length[bits];

  (* start at the left of the bits array.
     example: if n = 43, then n + 1 = (binary)101100; compute u1, u2, u4, u5, u10, u11, u22, u44 *)
  (* initialize u1, v1, and Q^1 *)
  uk = 1;
  vk = capP;
  qk = capQ;

  (* d = 1 and if capP = n, then v1 = v(d) is 0 (mod n), so already, n is a strong lprp *)
  If[ (d == 1) && (Mod[vk, n] == 0), strongLprp = 1 ];

  For[i = 2, i <= nBits, i++,
    uNew = uk * vk;      (* double the subscript *)
    vNew = vk * vk - 2 * qk;
    qPowerNew = qk^2;
    uk = Mod[uNew, n];
    vk = Mod[vNew, n];
    qk = Mod[qPowerNew, n];
    If[bits[[i]] == 1,  (* increase the subscript by 1 *)
      uTemp = capP * uk + vk;
      If[ ! EvenQ[uTemp], uTemp += n];  (* make uTemp be even *)
      uNew = uTemp / 2;
      vTemp = capD * uk + capP * vk;
      If[ ! EvenQ[vTemp], vTemp += n];
      vNew = vTemp / 2;
      qPowerNew = qk * capQ;
      uk = Mod[uNew, n];
      vk = Mod[vNew, n];
      qk = Mod[qPowerNew, n];
    ];

    If[ (i == nBits - s) && (uk == 0), strongLprp = 1 ];  (* u(d) == 0 *)
    If[ (i >= nBits - s) && (i < nBits) && (vk == 0), strongLprp = 1 ]; (* v(d*2^i) == 0 *)

  ];  (* end For i loop *)

  If[uk == 0, lprp = 1];

 { lprp, strongLprp }

];  (* end of module lprpTestA *)


lpspQ[n_?EvenQ] := False;
lpspQ[n_?PrimeQ] := False;
lpspQ[n_Integer?Positive] := (lprpTestA[n][[1]] == 1);

lpspList = {};
k = 3;
While[k < 100000,  (* there are 57 lpsp (323, 377, 1159, ..., 97439) up to 100000 *)
  If[lpspQ[k],
    Print[k];
    AppendTo[lpspList, k]
  ];
  k += 2
];
lpspList