6
$\begingroup$

Let's say that we have some $300\,\text{K}$ digits (arbitrary function) and want to trial factor with $100{,}000{,}000$ first prime numbers.

PRP = {}; 

Func[Q_] := 77^Q + 2;

Do[If[TimeConstrained[PrimeQ[Func[i]], 1, True], AppendTo[PRP, i]], {i, 160000, 161000}]

PRP: ($64$ numbers out of $1000$ test samples)

{160001,160025,160037,160041,160045,160053,160065,160097,160101,160133,160137,160141,160161,160177,160181,160201,160241,160253,160257,160261,160285,160297,160305,160317,160321,160361,160393,160397,160405,160413,160421,160425,160441,160485,160541,160557,160577,160581,160601,160621,160633,160637,160645,160657,160685,160693,160713,160721,160725,160753,160761,160801,160821,160833,160845,160853,160857,160861,160881,160897,160921,160933,160945,160965}

Then I use the following to compute product of primes in $5$ part ($20{,}000{,}000$ range)

Pro = Product[i, {i, Prime[Range[1, 20000000]]}]; // about 170 million digits

And for each product I run:

Output ={};

Do[If[GCD[Func[i], Pr] == 1, AppendTo[Output, i]], {i, PRP}]

Currently each product (of $20{,}000{,}000$ primes) take $1\,\text{GB}$ of RAM and nearly $30$ seconds for each GCD

In[86]:= Timing[GCD[Func[160001], Pr]]

Out[86]= {28.562500, 117414067}

Is my configuration and setup the best possible one in memory management and speed?

Update

I modified the @tchronis answer a little to remove composite as early as possible. Not the best code (can be much better using Reap and Sow), but it works:

ClearAll[gcd2];
gcd2[list_, noprimes_: 20000000, noclusters_: 1] := 
 Module[{clusters, n, pro}, 
  clusters = 
   Table[{(i - 1)*noprimes/noclusters + 1, i*noprimes/noclusters}, {i,
      noclusters}];
  innerList = list;
  Do[
   pro = Product[i, {i, Prime[Range @@ clusters[[j]]]}];
   newList = {};
   Do[
    If[GCD[func[m], pro] == 1, AppendTo[newList, m]],
    {m, innerList}
    ]; innerList = newList,
   {j, 1, noclusters}
   ]; newList
  ]

Result

Timing[gcd[list, 5000000, 4]]

{327.9706, { .... 54 out of 64 ... } }

Timing[gcd2[list, 5000000, 4]]

{290.328125, { .... 54 out of 64 ... } }

The higher the number of primes and the larger the list, the speedup by stopping computation upon first non-one GCD gains more.

Improve this question
$\endgroup$
8
  • $\begingroup$ Why are you taking the product since you are dealing with prime numbers? GCD[Func[i],Pr] will be equal to 1 if and only if GCD[Func[i],everyprime] result will be equal to 1 for every prime of those 20000000. In this way you will save a lot of Memory but maybe will be slower... $\endgroup$
    – tchronis
    Dec 15, 2013 at 12:02
  • $\begingroup$ @tchronis, if p be a factor of n, then GCD[p,n] != 1 $\endgroup$
    – Mohsen Afshin
    Dec 15, 2013 at 12:14
  • $\begingroup$ Yes that is true. I am talking about the theorem: Gcd(a,p1*p2*p3*...pn)==1 if and only if gcd(a,pi)==1 for all i from 1 to n where pi are distinct primes. $\endgroup$
    – tchronis
    Dec 15, 2013 at 12:20
  • $\begingroup$ @tchronis, I know that. The same operation using Divisible or GCD by all primes is 10 times slower (230 seconds) $\endgroup$
    – Mohsen Afshin
    Dec 15, 2013 at 12:29
  • 2
    $\begingroup$ Maybe you can split the product in let's say 4 sets and with a loop check them all. In that case you will use 1/4 of the memory and become less slower... $\endgroup$
    – tchronis
    Dec 15, 2013 at 12:32

2 Answers 2

Reset to default
8
$\begingroup$

Try to use lowercase names at least for the initial letter of your functions. My timings are different than yours. Probably your computer is much faster than mine.

In your code pro is calculated in 144 seconds and GCD[func[160001],pro] needs another 34 seconds so a total of 178 seconds. Of course pro is calculated only once but the same happens in the clustered products below.

prp = {};
func[q_] := 77^q + 2;
Do[If[TimeConstrained[PrimeQ[func[i]], 1, True], 
  AppendTo[prp, i]], {i, 160000, 161000}]

The following function divides in clusters the product to gain memory but also speed.

ClearAll[gcd];
gcd[list_List, noprimes_: 20000000, noclusters_: 1] :=
 Module[{clusters, n, pro},
  clusters = 
   Table[{(i - 1)*noprimes/noclusters + 1, i*noprimes/noclusters}, {i,
      noclusters}];

  Times @@ ParallelTable[
    pro = Product[i, {i, Prime[Range @@ clusters[[j]]]}];
    GCD[func[#], pro] & /@ list
    , {j, 1, noclusters}]
  ]

Notice : noclusters must divide exactly noprimes

So using 40 clusters I managed to gain speed and lot's of memory.

AbsoluteTiming[
 gcd[{160001}, 200000000, 40]
 ]

run in 140 seconds (in total for a single kernel).

You can use it also as : gcd[prp,20000000,40] to get all results.

Remark: You can also gain a speedup x3 (it runs in 45 seconds in my machine) in a quad core if you use ParallelTable inside the gcd function. Times @@ ParallelTable[pro = Product[i, {i, .....

Results: (parallel table on in 644 sec only) So now you can grow your search space.

{643.072782, {117414067, 1, 1002109, 1, 1, 86192723, 1, 1445419, 1, 1,
   1, 1, 1, 123791117, 1, 14474231, 100253533, 1, 1, 1, 113404303, 1, 
  1, 1227181, 1, 1, 1, 1, 1, 1, 27367889, 1, 1, 1, 1, 41178451, 
  10287491, 1, 1, 1, 367885753, 1604461, 1, 1, 1, 1, 196558247, 1, 1, 
  582153380112268780892543053, 152758451, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
  207751097, 117692293, 42702377, 1}}

Results: for 100M primes

{3332.700619, {117414067, 1, 1002109, 1, 1, 86192723, 1, 1445419, 1, 
  1, 1, 1, 1, 123791117, 1, 14474231, 167856571665498613, 1, 1, 1, 
  113404303, 1, 1, 1227181, 1, 1, 1097490881, 1, 1, 1, 27367889, 1, 
  1271310361, 1, 712991261, 41178451, 10287491, 1, 1, 1, 367885753, 
  1604461, 1758063169, 552414241, 1, 1, 196558247, 1, 1, 
  582153380112268780892543053, 152758451, 1, 1, 1635548081, 1, 1, 1, 
  1, 1, 1, 207751097, 117692293, 42702377, 1}}

Results: for 1 billion primes AbsoluteTiming[gcd[prp, 1000000000, 2000]]

{55026.855356, {117414067, 1, 1002109, 1, 1, 86192723, 1, 1445419, 1, 
  1, 1, 1, 1, 554101616002949381, 1, 14474231, 167856571665498613, 1, 
  15121883417, 1, 113404303, 1, 2071728343, 1227181, 1, 1, 1097490881,
   1, 1, 1, 27367889, 1, 7718516532572391937, 1, 712991261, 41178451, 
  10287491, 1, 1, 1, 367885753, 1604461, 1758063169, 552414241, 1, 1, 
  196558247, 1, 1, 582153380112268780892543053, 152758451, 1, 
  3743040013, 1635548081, 11228560151, 1, 1, 1, 1, 1, 207751097, 
  743891268843592201, 42702377, 1}}

It seems that for more primes (10 billion) there will be even less 1's

Improvements

The two versions below return the list of exponents that lead to GCD=1 (for the selected number of primes). In the second version gcd2 exponents that lead to GCD greater than 1 are dropped quickly. Also there is an additional parameter namely startcluster in case someone wants to start from greater primes than 2,3,...

ClearAll[gcd];
gcd[list_List, noprimes_: 20000000, noclusters_: 1, 
  startcluster_: 1] := 
 Module[{clusters, pro}, 
  clusters = 
   Table[{(i - 1)*noprimes/noclusters + 1, i*noprimes/noclusters}, {i,
      startcluster, noclusters}];

  Pick[list,
   Times @@ 
    ParallelTable[pro = Product[i, {i, Prime[Range @@ clusters[[j]]]}];
     GCD[func[#], pro] & /@ list, {j, 1, noclusters - startcluster +1}]
   , 1]]

I noticed that the code below does not exploit the CPU kernels efficiently so it seems some times inferior than "simple" gcd above. This happens because some work may left only for one kernel before moving to the next product.

ClearAll[gcd2];
gcd2[list_List, noprimes_: 20000000, noclusters_: 1, 
  startcluster_: 1] := 
 Module[{clusters, pro, res}, 
  clusters = 
   Table[{(i - 1)*noprimes/noclusters + 1, i*noprimes/noclusters}, {i,
      startcluster, noclusters}];

  res = Table[1, {Length[list]}];

  Pick[list,
   Times @@ 
    ParallelTable[pro = Product[i, {i, Prime[Range @@ clusters[[j]]]}];
     If[res[[#]] != 1, 1, res[[#]] = GCD[func[list[[#]]], pro]] & /@ 
      Range[Length[list]],
     {j, 1, noclusters - startcluster + 1}]
   , 1]]

Also tuning the number of clusters is a hard nut and needs some exploration time.

After all the above my approach would be using the new gcd (not the gcd2) and construct my clusters progressively (maybe using a Fold).

Improve this answer
$\endgroup$
4
  • $\begingroup$ @thchronis, Thank you for great response, but does it break the operation by the first GCD result greater than 1? Since if a number has a factor in of the chunks, then computing the rest is just a matter of time wasting as I don't require the exact GCD result but rather equality or inequality to 1. If I understood correctly, your code multiplies all (40) GCD results for each number. $\endgroup$
    – Mohsen Afshin
    Dec 15, 2013 at 17:02
  • $\begingroup$ @MohsenAfshin thanks, you are right it doesn't break the operation to save more time in case of a GCD greater than 1. It needs some modifications to do that although I am not so sure about how much this will benefit you (you see the most time consuming part is the product calculation - GCD is extremely fast). I can modify it if you like this but not today :-) I am running for 1 billion primes just from curiosity and will post an update in a few hours. $\endgroup$
    – tchronis
    Dec 15, 2013 at 17:08
  • $\begingroup$ @tchrnois, It's a tradeoff between trial factor and deterministic tests. 100M primes left 38 numbers while 1B primes left 34 numbers, only 4 numbers removed by the additional 14 hours of test while performing long deterministic tests on those 4 numbers would take 4 hours, so there's a waste of 9 hours. There should be a balance. $\endgroup$
    – Mohsen Afshin
    Dec 16, 2013 at 7:31
  • $\begingroup$ @MohsenAfshin you are right for long computations you need to drop quickly the numbers with certain GCD>1. Also my choice of 2000 clusters in this case slowed down the calculations. I am working on a new function now to improve it further. $\endgroup$
    – tchronis
    Dec 16, 2013 at 7:42
4
+50
$\begingroup$

I post here my latest and final code. I am adding an additional post to expose a clear code.

At first we define func and we put in prp list the candidate exponents.

prp = {};
func[q_] := 77^q + 2;
Do[If[TimeConstrained[PrimeQ[func[i]], 1, True], 
  AppendTo[prp, i]], {i, 160000, 161000}]

Then the following function returns the exponents from the list argument that are relative prime with all primes from starcluster first prime to noprimes. In this way you can omit already calculated ranges if you want to distribute the calculation to several computers.

ClearAll[gcd];
gcd[list_List, noprimes_: 10000000, noclusters_: 16, 
  startcluster_: 1] := 
 Module[{clusters, pro}, 
  clusters = 
   Table[{(i - 1)*noprimes/noclusters + 1, i*noprimes/noclusters}, {i,
      startcluster, noclusters}];
  Pick[list, 
   Times @@ 
    ParallelTable[pro = Product[i, {i, Prime[Range @@ clusters[[j]]]}];
     GCD[func[#], pro] & /@ list, {j, 1, 
      noclusters - startcluster + 1}], 1]]

Now for not carrying along exponents that have been characterized with GCD>1 the best approach I came through was a repetitive calculation in ranges (The function fold arranges them).

(*res stores all numbers that correspond to GCD\[Equal]1*)
res = prp;

(*This number must be a multiple of the number of available kernels*)
noclusters = 16;
(*noprimesinacluster/noclusters is bounded for the calculations to \
fit in memory*)
noprimesinacluster = 10000000;

ClearAll[fold];
fold[i_] := (res = 
    gcd[res, i*noprimesinacluster, 
     i*noclusters, (i - 1)*noclusters + 1]);

current = 0;
Manipulate[
 Dynamic[
  stop = u;
  Column[{current, res, Length[res], Complement[prp, res]}]
  ]
 , {{u, 10000, "Stop at ? x " <> ToString[noprimesinacluster]}, 1, 
  10000}]

Now with the following code you can calculate res using the first 30*10000000 primes. (Remember we have set noprimesinacluster = 10000000)

AbsoluteTiming[
 Catch[Scan[If[# > stop, Throw[#], current = #; fold[#]] &, 
   Range @@ {1, 30}]]
 ]

You can stop the evaluation by moving the slider above to 0 (as left as you can). You will wait for at most 1 minute for it to stop because it will be calculating the latest 100000000 group. If you want to continue you have to move the slider again to the right.

Finally I tested for the first 2 billion prime numbers and I got 33 exponents (from 64):

res={160025, 160041, 160045, 160065, 160101, 160133, 160137, 160141, \
160161, 160181, 160253, 160297, 160321, 160361, 160397, \
160405, 160413, 160425, 160485, 160581, 160601, 160621, 160685, \
160693, 160721, 160725, 160801, 160853, 160857, 160861, 160881, \
160897, 160965}

Maybe I should take the time to put more comments but if you run the code in the same order all gray areas will be lightened.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 318 seconds on my test sample. It seems that it reached the limits of speed and performance $\endgroup$
    – Mohsen Afshin
    Dec 17, 2013 at 19:54

Your Answer

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

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