How to further accelerate arithmetic with Fermat Pseudoprime and Fibonacci number

Question

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?

Answer 1

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.