Register to reply 
Sow within While Loop (Mathematica) 
Share this thread: 
#1
Feb1012, 05:28 PM

P: 894

I posted this in number theory "Ramsey Primes" but I don't think this qualifies as an improper use of this forum. I need help with my code. My code for Mathematical is short so I will include it below for reference. There are some flaws. For instance the Reap/Sow function generates a list of the form {Null,{{prime3,prime5,prime6,prime8,...}}}. I want to generate a simple list {prime3,prime5,prime6,prime8...}. Also, the code took about an hour to pick out 31 primes between 1049000 and 1050000. Overall, about 2/3 of the prime numbers are not picked out, but that is what I am looking for code that will generate only primes, not code that will generate all primes but also include some composites. So far some 30,000 prime numbers between 5 and 1050000 were generated with no false hits. Any help improving my code will be appreciated. The code follows:
Caa = 5 Reap[While[Caa < 1050001, Co = 1; S0 = 2; S1 = 3; While[Co < ((Caa1)/2),Temp = Mod[6 S1  S0  6,Caa];If[Temp <= 1,Break[]]; S0=S1;S1=Temp;Co++]; n = ((Caa1)/2)  Co; If[n<2,If[n>0, Sow[Caa]]]; Caa++;Caa++]] 


#2
Feb1012, 06:40 PM

P: 1,029

The documentation for Reap tells that it returns a list of two items, the first being the expression contained within the Reap and the second the list of items Sown, with sublists for each tag you used, including using no tag.
Your enclosed While returns the Null you are seeing and the nested list are your primes Sown with no tag. What I always write is Reap[expressionContainingSow][[2,1]] because I never care about the expression and almost never use tags. 


#3
Feb1012, 11:23 PM

P: 894

Caa = 1049001; Reap[While[Caa < 1050001, Co = (Caa  3)/2; S0 = 2; S1 = 3; While[Co > 0, Temp = Mod[6 S1  S0  6, Caa]; If[Temp < 1, Break[]]; S0 = S1; S1 = Temp; Co]; If[Co < 2, If[Co > 0, Sow[Caa]]]; Caa++; Caa++]] 


#4
Feb1112, 03:36 PM

P: 1,029

Sow within While Loop (Mathematica)
S0 = S1;S1 = Temp;Co]; If[Co == 1, Sow[Caa]]; Caa+=2]] 


#5
Feb1312, 04:33 PM

P: 894




#6
Feb1312, 09:30 PM

P: 1,029

That line needed a pair of () for precedence and it should have been: While[Co>0 && (Temp=Mod[6 S1S06, Caa])>=1,S0=S1; S1=Temp; Co]; I just tested this and believe that is correct now. Your code is very slow and spends almost all its time inside that inner while loop. You might try to find a more efficient way of computing that. Or you could try using an external C program through MathLink. Or you could try using Mathematica's Compile[] function, which has nothing to do with compiling external programs and is one of the least well documented features I ever use. Every time I try to write a Compile[] I want to tear my hair out. But after three attempts I offer you this. Note that I reduced the range it was searching through because I couldn't wait as long as it would have taken your original code to finish. In[1]:= Timing[ Caa=1049001; Reap[ While[Caa<1049101, Co=(Caa3)/2;S0=2;S1=3; While[Co>0&&(Temp=Mod[6 S1S06,Caa])>=1,S0=S1;S1=Temp;Co]; If[Co==1,Sow[Caa]]; Caa+=2 ] ] ] Out[1]= {336.904 Second,{Null,{{1049011,1049051,1049077}}}} In[2]:= Timing[ f=Compile[{{Co,_Integer},{S0,_Integer},{S1,_Integer},{Caa,_Integer}}, Module[{xCo=Co,xS0=S0,xS1=S1,Temp}, While[Temp=Mod[6 xS1xS06, Caa]; xCo>0 && Temp>=1;xS0=xS1;xS1=Temp;xCo]; xCo ] ]; Caa=1049001; Reap[ While[Caa<1049101, Co=(Caa3)/2;S0=2;S1=3; If[f[Co,S0,S1,Caa]==1,Sow[Caa]]; Caa+=2 ] ] ] Out[2]= {22.983 Second, {Null, {{1049011, 1049051, 1049077}}}} Using Compile[] is fraught with difficulties and limitations, but if you can live with these you see that you can have about 10x the speed and the same output. Probably the one limitation you should watch out for is that all the integer values within Compile should be within +/ 2^31. If you just strip the Timing[] off this then you should have the same result as your original code. Tiny changes in this can make big differences, I moved one statement and gained another 50% in speed. Type declarations in Compile[] for Temp and Mod[_] might help, or hurt, but I've never figured out how to do that correctly. But if you fiddle with this and break it then I probably won't be the one to fix it. If anyone would write the world's greatest set of clear guidelines for writing a Compile[] that will correctly do what I want on the first try then I might buy that person a cake. I hope this is correct. Please point out any errors. If you scrape that off the screen and paste it into Mathematica then watch out for the space inserted by the forum software between the I and n in {Caa,_Integer} above. I don't know why it does that but I've seen that before. 


#7
Feb1412, 10:11 AM

P: 894

S(0) = 2, S(1) = 6, S(N) = 6S(N1) S(N2); the test is that S((N1)/2) = 6 Mod N and no S(i) prior to S((N1)/2) = 6 mod N. The above recursive series has a direct formulation for the N th term and thus it is easy to first test if S((N1) /2) = 6 mod N prior to doing the second part of the test. By the way, I noted that this new test series is directly related to the Pell series in that if we denote P(0) = 1, P(1) = 2 P(n) = 2*P(n1) + P(n2) then S(n)*P(2n)1 = P(4n). 


#8
Feb1512, 11:36 AM

P: 894

In[2]:= Timing[ f=Compile[{{Co,_Integer},{S0,_Integer},{S1,_Integer},{Caa,_Integer}}, Module[{xCo=Co,xS0=S0,xS1=S1,Temp}, While[Temp=Mod[6 xS1xS06, Caa]; xCo>0 && Temp>=1;xS0=xS1;xS1=Temp;xCo]; xCo ] ]; Caa=1049001; Reap[ While[(Mod[Round[N[(3 + 2 Sqrt[2])^((Caa1)/2) + (3  2 Sqrt[2])^((Caa1)/2)]],Caa]) == Caa6 && Caa<1049101, Co=(Caa3)/2;S0=2;S1=3; If[f[Co,S0,S1,Caa]==1,Sow[Caa]]; Caa+=2 ] ] ] I am currently running old code in Mathematica and don't want to interrupt it as I have too much time invested in it. If you verify the timing and results of the modified code, I would greatly appreciate it. I used twice the direct formula from the Online Encyclopedia of Integer Sequences for A001541 i.e.{1,3,17,99,577,3366,19601} which is my old sequence with each term S_New(n) = 2 S_Old(n) 3. That is why the test is M(S(n),Caa)== Caa6 instead of ==Caa3. PS I am uncertain whether my modification will work with the compile function, If a further modification is needed please advise me. 


#9
Feb1512, 06:03 PM

P: 1,029

Unfortunately at this point I am completely confused.
If I run your latest modification I get Out[1]= {0. Second, {Null, {}}} This appears to be because the outer While condition fails on the first evaluation because (Mod[Round[N[(3 + 2 Sqrt[2])^((Caa  1)/2) + (3  2 Sqrt[2])^((Caa  1)/2)]], Caa]) == 544707 and that is not equal to Caa6 == 1048995. So you are never even getting to f. Perhaps we can simplify all this by going back to the beginning. Can you just list the first dozen terms of your sequence, not starting at the 1049001'th term? Can you show the sequence number of the existing OEIS sequence that this is a slight variation of? Perhaps with that we can find a very simple modification of the OEIS code to give you what you want. Perhaps not. 


#10
Feb1512, 11:07 PM

P: 894

As for my original test sequence, it is S = 2,3,10,51,... where S(n) = 6*S(n1)  S(n2)  6. As I said each term is multiplied by 2 and then 3 is subtracted to get S_New = 1,3,17,99,... S(n) = 6(S(n1)  S(n2) and the new test should be S((Caa1)/2) = 3 Mod Caa, or 2*S((Caa1)/2) = 6 Mod Caa. The sequence of Ramsey Primes begins 3,5,11,13,19,... but 11 or above has to be the starting number of my program since the first Temp term checked is 2*S(2) =2* S((51)/2) = 4 Mod 5. Although 4 == 56 Mod 5, I am not sure they are equal as far as the test goes. I am beginning to doubt whether the Mathematica code for A(n) in the Online Encyclopedia is correct. I am now going to try to see what is wrong. 


#11
Feb1512, 11:52 PM

P: 1,029

While[False,stuff] isn't going to do stuff. Did you perhaps get that from OEIS A001541 here? "MATHEMATICA Table[Round[N[(1/2) (3 + 2 Sqrt[2])^n + (1/2) (3  2 Sqrt[2])^n]], {n, 0, 20}] [From Artur Jasinski..."? But I don't understand what you are doing changing the OEIS code to generate something different. Can you just have a moment of patience for me, pretend to forget all the things you know about this and just tell me the first two terms of the sequence you really want and how to generate the nth term from the(n1)th and (n2)th terms? I'm trying to cut through all the things you are saying to get to exactly what is essential and no more. I can't tell from all the things you are saying what is absolutely important and what makes no difference. Is all you want just A001541 and if so then why are you making changes? And if A001541 isn't what you want then precisely what are the first 20 terms of what you want and how do you generate the next from the previous two? Perhaps imagine you are dealing with someone who doesn't have much concentration left and knows nothing about what you have done. You need to give them a recipe so they can leave the room, work, come back and give you exactly what you want. Every detail that is not essential for them to do this reduces your chance of getting what you want by 10%. Look at the first line or two of several OEIS sequences. That is a recipe. Or maybe you are looking for something more complicated than that. But I just can't tell. Thank you 


#12
Feb1612, 03:31 PM

P: 894

Bill, I am sorry for the confusion and will try to respond to each of your questions without the web site timing out on me as it just did and I lost my reply.
First there is no change in the test sequence or in the code within the inner loop. I am only trying to get a direct formulation for the (Caa1)/2 term of my origional series and test to see if that equals 0 mod Caa so we can avoid needless time spent within the inner loop. Since there is no changes made to the variables by the test there should be no changes to either the module or the inner loop even if we add the test. OEIS A001541 is related to my test sequence in that each term plus 3 is twice the corresponding term of my sequence. So I need to add 3 to the code for the (Caa1)/2 term and divide by 2 to get S((Caa1)/2). However the step of dividing by 2 is not neccessary since if Caa divides S((Caa1)/2) it will also divide twice that. Likewise, the steps of division in the code for the (Caa1)/2 th term of OEIS A001541 is not necessary but then 6 must be added to get 4*S((Caa1)/2) instead of adding 3 to get twice S((Caa1)/2. Others have noted that only primes of the form 8n +/ 3 are picked out by my test, but not all primes of that form. I am accepting that and so intend to run my code twice starting at 8n+3 and 8n+5 respectively over each range checked so Caa+=2 should be changed to Caa+=8. That change should double the speed. Yours truely, Ken 


#13
Feb1912, 08:12 PM

P: 894

Here is the modified code I talked about. It does increase the speed since it no longer searches numbers of the form 8n +/ 1.
Timing[ f = Compile[{{Co, _Integer}, {S0, _Integer}, {S1, _Integer}, {Caa, _Integer}}, Module[{xCo = Co, xS0 = S0, xS1 = S1, Temp}, While[Temp = Mod[6 xS1  xS0  6, Caa]; xCo > 0 && Temp >= 1, xS0 = xS1; xS1 = Temp; xCo]; xCo ] ]; P = 3000000; Caa = P+3; List3 = Reap[ While[ Caa < P+400000, Co = (Caa  3)/2; S0 = 2; S1 = 3; If[f[Co, S0, S1, Caa] == 1, Sow[Caa]]; Caa +=2 Co = (Caa  3)/2; S0 = 2; S1 = 3; If[f[Co, S0, S1, Caa] == 1, Sow[Caa]]; Caa +=6 ] ][[2,1]]; Export["aaazzz.csv",List3]; PrimeQ[List3] ] A search over the above range took took around 40,000 seconds to run. Still no composites returned after searching all numbers 8n +/ 3 below 3400000. 


Register to reply 
Related Discussions  
Mathematica: Summation in the DoLoop  Math & Science Software  25  
Mathematica Loop with Dot Product  Math & Science Software  3  
For Loop Problem in Mathematica  Math & Science Software  3  
Loop in Mathematica  Introductory Physics Homework  5  
Mathematica loop problem.  Math & Science Software  8 