Implementing the Atkin-Bernstein sieve efficiently is nontrivial. I wouldn't recommend it for you -- not until you've at least implemented
1. modular exponentiation by squaring
2. Miller-Rabin
3. the sieve of Eratosthenes
4. modular exponentiation with a sliding window
5. Karatsuba multiplication
as warm-ups. You'll want to work in some suitably lowish-level language, certainly no higher than C++ (not gp, not Java, not Perl).

[QUOTE=CRGreathouse;265046]3. the sieve of Eratosthenes[/QUOTE]


[CODE]Eratosthenes(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),for(y=1,#b,if(b[y]!=0&&b[y]%a[#a]==0,b[y]=0));for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));for(d=a[#a],#b,if(b[d]!=0,a=concat(a,b[d])));a[/CODE]
is this close enough or do I have to find 2 prime as well ?

[QUOTE=science_man_88;265058][CODE]Eratosthenes(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),for(y=1,#b,if(b[y]!=0&&b[y]%a[#a]==0,b[y]=0));for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));for(d=a[#a],#b,if(b[d]!=0,a=concat(a,b[d])));a[/CODE]
is this close enough or do I have to find 2 prime as well ?[/QUOTE]

yeah I know I didn't optimize it lol.

[QUOTE=science_man_88;265058]is this close enough[/QUOTE]

It's not the sieve of Eratosthenes. You do
[code]b[y]!=0 && b[y]%a[#a]==0[/code]
about [TEX]x\pi(\sqrt x)\approx x^{3/2}/\log x[/TEX] times, which is more operations right there than the entire SoE uses. Critically, though, you shouldn't be dividing in the inner loop at all -- you should just be adding and marking 0s.

But again, GP isn't the right language for this -- as good as it is on high-level operations it's quite poor at low-level ones like b[y]=0. Consider what it actually does when faced with that statement:
1. Parse the bytecode into machine instructions
2. Go to the location of y and read its codewords.
3. Do some quick bitwise ops to check the type of y.
4. Do some quick bitwise ops to find the location of the length of y.
5. Check that the length of y is at most one machine word.
6. Check that the sign of y is positive.
7. Do some quick bitwise ops to find the value of the first word of the value of y.
8. Read this value into a register, treating it as an unsigned long.
9. Go to the location of b and read its codewords.
10. Do some quick bitwise ops to check the type of b.
11. Do some quick bitwise ops to find the location of the length of b.
12. Check that a is at most the length of b.
13. Do some quick arithmetic to calculate the address pointed to by b[y]
14. Check if the value pointed to by b[y] is at the end of the PARI stack. (?)
15. If so, move the stack pointer avma. (?)
16. Overwrite b[y] with a precomputed value that refers to gen_0.

If the value was anything other than one of the special values, a new location would need to be created, the value written, and avma moved. But fortunately I don't think that has to happen in this case.

Compare this to what C would have done:
1. Do some quick arithmetic to calculate the address pointed to by b[y]
2. Set this equal to 0.

[QUOTE=CRGreathouse;265070]It's not the sieve of Eratosthenes. You do
[code]b[y]!=0 && b[y]%a[#a]==0[/code]
about [TEX]x\pi(\sqrt x)\approxx^{3/2}/\log x[/TEX] times, which is more operations right there than the entire SoE uses. Critically, though, you shouldn't be dividing in the inner loop at all -- you should just be adding and marking 0s.[/QUOTE]

it's not optimal but it does the exact thing the sieve does and gets the correct answer.

[QUOTE=science_man_88;265071]it's not optimal but it does the exact thing the sieve does and gets the correct answer.[/QUOTE]

but I'll do it over .

[QUOTE=science_man_88;265071]it's not optimal but it does the exact thing the sieve does and gets the correct answer.[/QUOTE]

No, actually it doesn't. The key operation in your function is the modular division % I quoted, but that operation is barely used at all in the SoE. In fact, it doesn't even need to use it at all!

What you need to do is, for each prime p <= sqrt(x) is to go from p^2, p^2 + p, p^2 + 2p, ..., all the way up to x and mark each of these values as 0.


By the way, I edited some information into #2248 explaining why GP isn't good at this. But still, I think you'll find that if you make the change I suggested above and replace the vector with a vectorsmall you can get the program to be at least 1000 times faster. :)

[QUOTE=science_man_88;265073]but I'll do it over .[/QUOTE]

[CODE]erastosthenes(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),forstep(y=a[#a],#b,a[#a],b[y]=0);for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));a[/CODE]

is this good enough or did I forget something else ? I know an optimization on this if it's correct though.

[QUOTE=science_man_88;265076]is this good enough or did I forget something else ?[/QUOTE]

The core of it works properly. Unfortunately your copying between arrays isn't working quite right and it's giving the primes as 2, 3, 3, 3, 3, ....

If your optimization involves not copying between arrays (very, very slow in GP) then why don't you do that and see if it fixes that problem. Two birds with one stone and all that.

[QUOTE=CRGreathouse;265075]No, actually it doesn't. The key operation in your function is the modular division % I quoted, but that operation is barely used at all in the SoE. In fact, it doesn't even need to use it at all!

What you need to do is, for each prime p <= sqrt(x) is to go from p^2, p^2 + p, p^2 + 2p, ..., all the way up to x and mark each of these values as 0.


By the way, I edited some information into #2248 explaining why GP isn't good at this. But still, I think you'll find that if you make the change I suggested above and replace the vector with a vectorsmall you can get the program to be at least 1000 times faster. :)[/QUOTE]

is that so:

[CODE]
l(a[#a]>sqrt(x),forstep(y=a[#a],#b,a[#a],b[y]=0);for(z=1,#b,if(b[z]!=
= (x)->a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),forstep(y=a[#a],
47)>Eratosthenes(10)
* user interrupt after 33,172 ms.
48)>Eratosthenes3(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),f
= (x)->a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),for(y=1,#b,if(b[
48)>Eratosthenes3(10)
= [2, 3, 5, 7]
48)>##
* last result computed in 0 ms.
49)>[/CODE]

actually my old code gave me a decrease of at least 1/33000 by the looks of these.

[QUOTE=science_man_88;265081]is that so:[/QUOTE]

Yes, I had to interrupt the code to see what it was doing as it had fallen into an infinite loop. If your third version works, good for you.