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, 
      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), "."];

            "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?

  • $\begingroup$ 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. $\endgroup$
    – Histograms
    Commented May 25, 2015 at 15:33
  • $\begingroup$ 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. $\endgroup$ Commented May 25, 2015 at 15:41
  • $\begingroup$ 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 :( $\endgroup$
    – Histograms
    Commented May 25, 2015 at 16:05
  • $\begingroup$ It is known that any such number is greater than $2^{64}$ $\endgroup$ Commented May 25, 2015 at 16:08
  • $\begingroup$ 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 $\endgroup$
    – Histograms
    Commented May 25, 2015 at 16:24

1 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[
     Divisible[Fibonacci[# + 1], #] &,
     Divisible[2^(# - 1) - 1, #] &

And it should be much faster.

  • $\begingroup$ 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. $\endgroup$ Commented May 25, 2015 at 15:51
  • 1
    $\begingroup$ One more thing: I'd use PowerMod[] instead of using Divisible[] on a big pseudoprime for testing Fermat's little theorem. $\endgroup$ Commented May 25, 2015 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.