methodC[n_Integer?Positive] :=
Module[
  { i, js, capQ = 1, capP = 3, capD },
  If[Floor[Sqrt[n]]^2 == n, Return[{ 0, 0, 0 }] ];
  For[i = 1, i <= 1000, i++,
    capD = capP^2 - 4;
    js = JacobiSymbol[capD, n];
    If[js == 0,
      (* if capD == n, (e.g., n = 5) try the next capD instead of quitting;
         this enables 5 to be classified as a lucas probable prime.
      *)
      If[Abs[capD] != n, Return[{ 0, 0, 0 }] ]
    ];
    If[js == -1, Return[{ capD, capP, capQ }] ];
    capP++
  ];
  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 }
];

extraStrongLprpTest[n_Integer?Positive] :=
Module[
  (* given n, this returns a list of three 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.
       if n is an extra strong lucas probable prime, then the third number is 1; otherwise, it is 0.
       therefore, if n is an odd prime, this will return { 1, 1, 1 }.
       if n < 3, is even, or is a square, then return { 0, 0, 0 }.
     this is based on "Frobenius Pseudoprimes" by Jon Grantham, in
     Mathematics of Computation, volume 70, number 234 (March, 2000), pp. 873-891.
     note: a lucas pseudoprime is a lucas probable prime that is composite.
  *)
  { lprp = 0, strongLprp = 0, extraStrongLprp = 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, 0 }]
  ];
  { capD, capP, capQ } = methodC[n];
  If[capD == 0, Return[{ 0, 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;

  (* if d = 1 and capP = n, then v1 = vk is 0 (mod n), so already, n is a strong and an extra strong lprp *)
  If[ (d == 1) && (Mod[vk, n] == 0),
    strongLprp = 1;
    extraStrongLprp = 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),  (* u(d) == 0 *)
      strongLprp = 1;
      If[ (Mod[vk - 2, n] == 0) || (Mod[vk + 2, n] == 0),
        extraStrongLprp = 1  (* u(d) == 0 and v(d) = 2 or -2 (mod n) *)
      ]
    ];

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

      If[i < nBits - 1,  (* see theorem 2.3 in the above reference *)
        extraStrongLprp = 1;
      ];
    ];

  ];  (* end For i loop *)

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

 { lprp, strongLprp, extraStrongLprp }

];  (* end of module extraStrongLprpTest *)


extraStrongLpspQ[n_?EvenQ] := False;
extraStrongLpspQ[n_?PrimeQ] := False;
extraStrongLpspQ[n_Integer?Positive] := (extraStrongLprpTest[n][[3]] == 1);

extraStrongLpspList = {};
k = 3;
While[k < 100000,  (* there are 12 extra strong lpsp (989, 3239, 5777, ..., 75077) up to 100000 *)
  If[extraStrongLpspQ[k],
    Print[k];
    AppendTo[extraStrongLpspList, k]
  ];
  k += 2
];
extraStrongLpspList