# Fast Trial Division

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}

-
It's possible to have it do the divisibility checks in parallel at least. – ssch Nov 29 '12 at 15:23
Don't use variable names like Max, since it is a built-in function ! – Artes Nov 29 '12 at 16:29
I think the Miller-Rabin type of testing makes more sense than endless divisibility tests, in terms of efficiency. – Daniel Lichtblau Nov 29 '12 at 18:49

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 *)

-
Did my suggestion to use HoldAll settle your issue completely? I wondered about the output {2.5453935, True} you added to my answer. Have a look at your answer where Out[84]={0.609, True}. Or is it just another machine? – Sascha Nov 30 '12 at 15:20
Timing is 2.5026742. My guess is that if the body of function is executed directly with no function call, timing would be equal to function usage with HoldAll. Maybe the expanded version is HoldAll by default. – Mohsen Afshin Nov 30 '12 at 15:44

## 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}]]

-
On my machine this takes approximately 2.6 seconds while Mohsens solution with the additions I made takes around 2.15 seconds. – Sascha Nov 30 '12 at 18:04
Curious note: First time evaluation of Mohsens solution takes approximately 2.7 seconds while multiple evaluation afterwards take only the 2.15 seconds mentioned above. – Sascha Nov 30 '12 at 18:10
@Sascha, I get different timings here. Please see my comments here. – halirutan Nov 30 '12 at 19:40