How to further accelerate arithmetic with Fermat Pseudoprime and Fibonacci number

I've been working on this all night, and I have made this go pretty fast, compared to my first iteration of the program, but now I'm out of ideas.

I'm trying to write a program to test (by good old fashioned brute force) if there is a composite positive integer $n \equiv 3 \pmod {10}$ or $n \equiv 7 \pmod {10}$ where $n$ is a Fermat Pseudoprime, and $n$ divides the $n+1^{\rm th}$ Fibonacci number.

Here is what I have so far:

$HistoryLength = 5; startingIter = 0; iters = startingIter; listLength = 10000; maxIters = 20; nextBlock[start_, listlength_] := Block[ {a, i}, a = DeleteCases[ Table[(3 + 10 i), {i, start, start + listlength - 1}], _?PrimeQ]; a = Join[a, DeleteCases[ Table[(7 + 10 i), {i, start, start + listlength - 1}], _?PrimeQ]]; a = DeleteCases[a, Except[_?(Divisible[Fibonacci[# + 1], #] &)]]; a = DeleteCases[a, Except[_?(Divisible[2^(# - 1) - 1, #] &)]]; Print["Found ", Length[a], " candidates (3 + 10n) or (7 + 10n) between ", TraditionalForm[start], " and ", (start + listLength - 1), "."]; a ]; Monitor[ Do[ Print[ "Iteration: ", ToString[iters + 1], ". Currently testing around: ", ToString[7 + 10*(1 + iters*listLength)], " and a Fibbonaci number with : ", Length[IntegerDigits[Fibonacci[7 + 10*(1 + iters*listLength)]]], " digits: ", nextBlock[1 + iters*listLength, listLength] ], {iters, startingIter, maxIters}], ProgressIndicator[iters, {startingIter, maxIters}] ]  The strategy here is that each time nextBlock[] is called, it first creates a list of all composite integers in the desired range which are congruent to$3$or$7 \mod 10$. Next, it removes any which do not divide the appropriate Fibonacci number, and then removes any that are not Fermat Pseudo-prime. It's fast enough on my system to do 20 iterations in a reasonable amount of time. Basically I have no clue what to do next to speed this up. I've looked into Compile[] but the Fibonacci numbers in question have hundreds of thousands of digits, so I can't do that. If anyone is good with parallelization, I could use some guidance since the documentation is difficult to follow. I have an i7 CPU and access to another PC with Mathematica and 4 cores, so if I could figure that out I would probably speed things up a bunch. I haven't done anything with GPU acceleration in many years, is it possible to port GMP onto CUDA or OpenCL? • Since it's a brute force search anyway, you might as well use FindInstance. Could you verify that this is equivalent: FindInstance[ (Mod[n, 10] == 3 || Mod[n, 10] == 7) && Mod[Fibonacci[n + 1], n] == 0 && Mod[2^(n - 1) - 1, n] == 0 && Not[Element[n, Primes]] && n > 0, n, Integers, 1] It will probably sit there for billions of years or forever. Are you sure there is a single number that satisfies these constraints? If you want a faster solution you'll have to use a more intelligent algorithm than brute force. May 25, 2015 at 15:33 • It's an open problem. Such a number would be a counterexample to a certain probabilistic Primality test. Really, this is more an interesting exercise for me to get better with Mathematica. May 25, 2015 at 15:41 • Well, I just ran FindInstance for a good length of time and it failed to either prove no such number existed or find any number... so my approach doesn't work :( May 25, 2015 at 16:05 • It is known that any such number is greater than$2^{64}\$ May 25, 2015 at 16:08
• ahh, well that explains it : Fibonacci[2^64 + 1] gives Overflow[ ] I'm thinking you might be able to leverage the Fibonacci square test somehow to avoid calculating the Fib altogether, if either (5 (k+1)^2 + 4) or (5 (k+1)^2 - 4) is square, k+1 is a Fibonacci number... tricky May 25, 2015 at 16:24

1. Never use pattern matching unless you absolutely have to. Using Cases instead of Select can make a huge difference.
2. Vectorize operations. Use Range instead of Table if you can.
3. Test several things at once. And[test1, test2, test3] will abort when it can for maximum efficiency ("short circuit evaluation").

Taking this into consideration your code looks like this:

nextBlock2[start_, listlength_] := Module[{a},
a = 10 Range[start, start + listlength - 1];
a = Join[3 + a, 7 + a];
Select[a, And[
Not@*PrimeQ,
Divisible[Fibonacci[# + 1], #] &,
Divisible[2^(# - 1) - 1, #] &
]]
]


And it should be much faster.

• This is great, thanks for your help! The only other change I made was to switch the order of Fermat-Pseudoprime testing and fibonacci testing since the fibonacci numbers are so big, so I'm eliminating more with the first test. May 25, 2015 at 15:51
• One more thing: I'd use PowerMod[] instead of using Divisible[] on a big pseudoprime for testing Fermat's little theorem. May 25, 2015 at 16:10