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 *) strongLpspQ[n_?EvenQ] := False; strongLpspQ[n_?PrimeQ] := False; strongLpspQ[n_Integer?Positive] := (lprpTestA[n][[2]] == 1); strongLpspList = {}; k = 3; While[k < 100000, (* there are 12 strong lpsp (5459, 5777, 10877, ..., 97439) up to 100000 *) If[strongLpspQ[k], Print[k]; AppendTo[strongLpspList, k] ]; k += 2 ]; strongLpspList