Sow within While Loop (Mathematica)

[ramsey2879]

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 < ((Caa-1)/2),Temp = Mod[6 S1 - S0 - 6,Caa];If[Temp <= 1,Break[]];
S0=S1;S1=Temp;Co++];
n = ((Caa-1)/2) - Co; If[n<2,If[n>0, Sow[Caa]]];
Caa++;Caa++]]

 

[Bill Simpson]

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.

 

[ramsey2879]

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.

Thanks for the help. I used your suggestion and also cut the time in half by rewriting the code as follows to avoid variable n:
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++]]

 

[Bill Simpson]

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++]]

While[Co > 0 && Temp = Mod[6 S1 - S0 - 6, Caa] >= 1,
S0 = S1;S1 = Temp;Co--];
If[Co == 1, Sow[Caa]];
Caa+=2]]

 

[ramsey2879]

While[Co > 0 && Temp = Mod[6 S1 - S0 - 6, Caa] >= 1,
S0 = S1;S1 = Temp;Co--];
If[Co == 1, Sow[Caa]];
Caa+=2]]

I get a set write error, i.e. "And" in "Co > 0 && Temp" is protected. What am I doing wrong?

 

[Bill Simpson]

I get a set write error, i.e. "And" in "Co > 0 && Temp" is protected. What am I doing wrong?

I apologize. I wrote that without testing it. I try to never do that.

That line needed a pair of () for precedence and it should have been:

While[Co>0 && (Temp=Mod[6 S1-S0-6, 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=(Caa-3)/2;S0=2;S1=3;
While[Co>0&&(Temp=Mod[6 S1-S0-6,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 xS1-xS0-6, Caa]; xCo>0 && Temp>=1;xS0=xS1;xS1=Temp;xCo--];
xCo
]
];
Caa=1049001;
Reap[
While[Caa<1049101,
Co=(Caa-3)/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.

 

[ramsey2879]

I apologize. I wrote that without testing it. I try to never do that.

That line needed a pair of () for precedence and it should have been:

While[Co>0 && (Temp=Mod[6 S1-S0-6, 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=(Caa-3)/2;S0=2;S1=3;
While[Co>0&&(Temp=Mod[6 S1-S0-6,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 xS1-xS0-6, Caa]; xCo>0 && Temp>=1;xS0=xS1;xS1=Temp;xCo--];
xCo
]
];
Caa=1049001;
Reap[
While[Caa<1049101,
Co=(Caa-3)/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.

Thanks For your work, it is very much appreciated. I was also looking for a way to speed up my code and came up with the following equivalent recursive series and test.
S(0) = 2, S(1) = 6, S(N) = 6S(N-1) -S(N-2); the test is that S((N-1)/2) = -6 Mod N and no S(i) prior to S((N-1)/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((N-1) /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(n-1) + P(n-2) then S(n)*P(2n)-1 = P(4n).

 

[ramsey2879]

I apologize. I wrote that without testing it. I try to never do that.

That line needed a pair of () for precedence and it should have been:

While[Co>0 && (Temp=Mod[6 S1-S0-6, 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=(Caa-3)/2;S0=2;S1=3;
While[Co>0&&(Temp=Mod[6 S1-S0-6,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 xS1-xS0-6, Caa]; xCo>0 && Temp>=1;xS0=xS1;xS1=Temp;xCo--];
xCo
]
];
Caa=1049001;
Reap[
While[Caa<1049101,
Co=(Caa-3)/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.

Here is my modification based upon a new Initial Test:


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 xS1-xS0-6, Caa]; xCo>0 && Temp>=1;xS0=xS1;xS1=Temp;xCo--];
xCo
]
];
Caa=1049001;
Reap[
While[(Mod[Round[N[(3 + 2 Sqrt[2])^((Caa-1)/2) + (3 - 2 Sqrt[2])^((Caa-1)/2)]],Caa]) == Caa-6 && Caa<1049101,
Co=(Caa-3)/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)== Caa-6 instead of ==Caa-3.
PS I am uncertain whether my modification will work with the compile function, If a further modification is needed please advise me.

 

[Bill Simpson]

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 Caa-6 == 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.

 

[ramsey2879]

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 Caa-6 == 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.

Bill, I think my initial test fails but I don't know why. Also, I think I found a typo in the module as well. Isn't the semicolon before xS0 = xS1 supposed to be a coma? I had to do that change to get your code to work. Stil I could not get my initial test to work and had to delete it. I haven't given up yet as your code runned in the range 349049001-349049499 took 11001.1 sec to pick out the numbers 340049059, '069, '179, '291, '293, '317, '333, '347, '347, '411 and '443 as prospective primes. All checked out as prime.
As for my original test sequence, it is S = 2,3,10,51,... where S(n) = 6*S(n-1) - S(n-2) - 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(n-1) - S(n-2) and the new test should be S((Caa-1)/2) = -3 Mod Caa, or 2*S((Caa-1)/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((5-1)/2) = 4 Mod 5. Although 4 == 5-6 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.

 

[Bill Simpson]

Bill, I think my initial test fails but I don't know why.

I can't tell where you came up with this

(Mod[Round[N[(3 + 2 Sqrt[2])^((Caa - 1)/2) + (3 - 2 Sqrt[2])^((Caa - 1)/2)]], Caa]) == Caa - 6

from but that is where I would start looking first. That is simply false before anything else happens and
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.

Also, I think I found a typo in the module as well. Isn't the semicolon before xS0 = xS1 supposed to be a coma? I had to do that change to get your code to work.

If I think I am looking at what you are then that is a good find, you are correct. That was my error introduced as I was fumbling around trying to get Compile[] to work, moving statements back and forth, dealing with it speeding up, slowing down, failing with a stream of errors, etc, etc. So that was my mistake and thank you for catching that.

Still I could not get my initial test to work and had to delete it.

I tried to warn you in a previous post about fiddling around with Compile. I don't know what to say and I'm still not sure I understand where to begin after reading your last post.

I haven't given up yet as your code runned in the range 349049001-349049499 took 11001.1 sec to pick out the numbers 340049059, '069, '179, '291, '293, '317, '333, '347, '347, '411 and '443 as prospective primes. All checked out as prime.

Ok.

As for my original test sequence, it is S = 2,3,10,51,... where S(n) = 6*S(n-1) - S(n-2) - 6.

So does your original sequence begin 2, 3, 10, 51, 290, 1683, 9802, 57123, 332930, 1940451? And is that not in the OEIS catalog? And if this is not what you really want then why is this being described? I'm more confused with each new post.

As I said each term is multiplied by 2 and then 3 is subtracted to get S_New = 1,3,17,99,...

So does your new sequence begin 1, 3, 17, 99, 577, 3363, 19601, 114243, 665857, 3880899? And is that OEIS sequence A001541? And this is what you want to generate? Or do you want something else?

S(n) = 6(S(n-1) - S(n-2) and the new test should be S((Caa-1)/2) = -3 Mod Caa, or 2*S((Caa-1)/2) = -6 Mod Caa.

Is this a brand new sequence that has roughly nothing to do with A001541? And if so then why was A001541 brought up? Can you just write down S(n)=and give the function of S(n-1) and S(n-2) to generate this?

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((5-1)/2) = 4 Mod 5. Although 4 == 5-6 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.

It is not impossible that code in the OEIS might be incorrect, but I would be astonished to see that for a sequence like A001541 with that many references and contributors. There are cranks and nut jobs who have found OEIS and submit dozens or hundreds of junk sequences just to get their name printed somewhere. Those often have mistakes in them. But I would think that for this sequence it is some mistake or misunderstanding that we have made.

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(n-1)th and (n-2)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

 

[ramsey2879]

Bill, I am sorry for the confusion and will try to respond to each of your questions without the website 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 (Caa-1)/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 (Caa-1)/2 term and divide by 2 to get S((Caa-1)/2). However the step of dividing by 2 is not necessary since if Caa divides S((Caa-1)/2) it will also divide twice that. Likewise, the steps of division in the code for the (Caa-1)/2 th term of OEIS A001541 is not necessary but then 6 must be added to get 4*S((Caa-1)/2) instead of adding 3 to get twice S((Caa-1)/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

 

[ramsey2879]

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.