Fast Trial Division

Question

I'm going to write a custom Trial Division primality test. I know that PrimeQ will first try trial division and then switches to PowerMod.

TrialFactorFreeQ[N_, Max_] := 
(For[j = 1, j < Max && Divisible[N, Prime[j]] == False, j++]; 
 Return[j == Max])

For example for testing a 300K digits number with first 100,000 primes, it took 4 seconds which is very slow since a real testing is applied for billions of primes.

In[179]:= Timing[TrialFactorFreeQ[3^1000000 + 2, 100000]]

Out[179]= {4.031, True}

Can it be optimized?

---

I wrote the parallel version and its timing seems good but when I wrap it in a function, the timing goes to sky:

p = 3^1000000 + 2;
max = Prime[100000];

In[85]:= ParallelTrialFactorFreeQ[Num_, Maxi_] := (IsFree = True; 
  ParallelDo[
    If[Divisible[Num, i], IsFree = False; Break[]], {i, 
     Prime[Range[1, PrimePi[Maxi]]]}] If[! IsFree, AbortKernels[]]; 
  Return[IsFree])

In[82]:= AbsoluteTiming[ParallelTrialFactorFreeQ[p, max]]

Out[82]= {18.2821829, True}

In[84]:= Timing[Num = p; Maxi = max;
     IsFree = True; 
     ParallelDo[
       If[Divisible[Num, i], IsFree = False; Break[]], 
       {i,Prime[Range[1, PrimePi[Maxi]]]}];
       If[! IsFree, AbortKernels[]]; IsFree]

Out[84]= {0.609, True}

---

Updated

Based on the Daniel Lichtblau answer to my Fast Sieve Implementation question, I wrote a very fast Trial Division function which can be used for numbers that are larger than products of primes in the given range:

prod = Product[i, {i, Prime[Range[10^5]]}];

TrialFactorFreeQ2[n_] := GCD[n, prod] == 1

And a 10x speedup

n = 3^1000000 + 2;

Timing[TrialFactorFreeQ2[n]]

{0.375, True}

Answer 1

For me

 SetAttributes[ParallelTrialFactorFreeQ, HoldAll]

does the trick and reduces the measured time your function needs to evaluate down to what you would expect.

Additionally I noticed that specifying the method used by ParallelDo brings down the timing slightly. By try and error i figured that Method -> "CoarsestGrained" works best on my machine. For other methods have a look in the documentation of Parallelize in the MORE INFORMATION- section. Code looks as follows:

p = 3^1000000 + 2;
maxIndex = 100000;

.

ParallelTrialFactorFreeQ[Num_, MaxIndex_] := (
  IsFree = True;
  ParallelDo[
    If[Divisible[Num, i], IsFree = False; Break[]], 
    {i, Prime[Range[1, MaxIndex]]}, 
    Method -> "CoarsestGrained"] 
  If[! IsFree, AbortKernels[]];
  Return[IsFree])

.

SetAttributes[ParallelTrialFactorFreeQ, HoldAll]

.

AbsoluteTiming[ParallelTrialFactorFreeQ[p, maxIndex]]

{2.611613, True} (* First Run *)

{2.0153128, True} (* Second Run *)

Answer 2

Preface

I want to emphasis the comment of @DanielLichtblau

I think the Miller-Rabin type of testing makes more sense than endless divisibility tests

Furthermore, I won't go into discussing your programming style with setting global variables in Modulearized functions, using parameter-names like N which are clearly built-in symbols, etc.

Solution

I will use two things in my solution

1. When you have $KernelCount parallel subkernels available you can parallelize it as follows: The first kernel starts with the first prime and makes steps of $KernelCount. The second kernel starts with the second prime and has the same step-size. Therefore, on a subkernel a loop is running which large steps (usually of size 4, when you have 4 subkernels).

2. ParallelTry[f,{arg1, arg2, ...}] has the nice property, that it stops when the first kernel comes up with a valid result. Invalid is, when f returns $Failed. Therefore, we can come up with an f which takes a start prime and the step-size and tests all numbers in a loop. If it finds a valid divisor, it returns False (because the possible prime is none) and otherwise $Failed. With this we ensure that ParallelTry runs as long as no subkernel found a divisor.

First we define f

f[{possPrime_, start_, end_, step_}] := Module[
  {result = $Failed},
  Do[If[Divisible[possPrime, Prime[i]],
    result = False;
    Break[]
    ],
   {i, start, PrimePi[end], step}];
  result
  ]

DistributeDefinitions[f]

In TrialFactorPrimeQ we first create as many parameter sets for f as we have kernels. Note, that I use max as input for NextPrime and therefore, max is the close to the maximum number we test. In your two examples you handle it not consistently. Since ParallelTry throws a message when no kernel comes up with a result, we have to check this:

TrialFactorPrimeQ[possPrime_, max_] := With[
  {chunks = Table[{possPrime, i, NextPrime[max], $KernelCount}, {i, $KernelCount}]},
  Quiet[
   Check[ParallelTry[f, chunks], True],
   ParallelTry::toofew
   ]
  ]

Testing

LaunchKernels[]
p = 3^1000000 + 2;
max = Prime[100000];

AbsoluteTiming@TrialFactorPrimeQ[p, max]

This runs about 2.5 seconds here, but you should test it on your machine. Compared to the serial version which needs 11 seconds

AbsoluteTiming[f[{p, 1, max, 1}]]