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
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 *)
Comments
Post a Comment