PARI's commands
 

As I posted elsewhere, I only know the rather basic commands of PARI:

Here are some examples of what I am able to do:

a(n) = {
n = n^n
}

Which, in PARI, generates n raised to the nth power.
Or, alternatively:

a(n) = {
n = n^(random(10^3)+2*10^3)
}

Which generates n raised to a random power between 2000 and 3000.

Or: a(n) = {
n = nextprime(n^n)
}
Which gives the smallest prime number greater than n[sup]n[/sup]
Any tips on other beginner/intermediate-level commands?

Admittedly, those 3-lined commands allow for many, many functions.

Loops are important.
[code]for(n=1,10,print(n))[/code]
[code]sumUpTo(n)={
\\ Faster version: n*(n-1)/2
sum(k=1,n,k)
};[/code]
[code]prod(i=1,20,i) \\ Faster version: 20![/code]
[code]forprime(p=2,97,print(p, "^2 = ", p^2))[/code]
[code]fordiv(120,d, \\ Notice different syntax!
print1(d" ")
)[/code]

if(condition, doIfTrue, doIfFalse) is useful. If the condition evaluates to zero or [0] or something like that, it runs doIfFalse; otherwise it runs doIfTrue. If doIfFalse is omitted, do nothing in that case.
[code]forprime(p=2,10,if(14%p == 0, print(p" divides 14")))[/code]
[code]if(7^2-49, print("this never happens"), print("7^2 minus 49 is 0"))[/code]

CRG: For prime numbers: The prod function can be useful:

prod(i=1, 20, prime(i))?

[QUOTE=3.14159;222279]CRG: For prime numbers: The prod function can be useful:

prod(i=1, 20, prime(i))?[/QUOTE]

See my comment on the other thread. This works but is generally inefficient. In this particular case it doesn't matter, since the numbers are large enough that the multiplications take far longer than the overhead -- but in general, especially when going to large primes, this is to be avoided.

Guess for printing n^2 + 1 primes:

if(isprime(n^2+1) = 1, print(n))?

[QUOTE=3.14159;222285]Guess for printing n^2 + 1 primes:

if(isprime(n^2+1) = 1, print(n))?[/QUOTE]

Close. "Equals" is ==, not =; a single = sign is assignment. "a = 2" puts the value 2 into the variable a; "a == 2" checks to see if a is equal to 2. So
[code]if(isprime(n^2+1) == 1, print(n))[/code]
will work. But there's no need for the " == 1" part: isprime already returns the appropriate thing (zero if not prime, nonzero if prime). So
[code]if(isprime(n^2+1), print(n))[/code]
is better.

Alright: Loops:

In general: for(n=x[sub]1[/sub],x[sub]2[/sub] print(f(n))?

Products: prod(i=x[sub]1[/sub],x[sub]2[/sub], i)?

[QUOTE=3.14159;222287]Alright: Loops:

In general: for(n=x[sub]1[/sub],x[sub]2[/sub] print(f(n))?[/QUOTE]

You can do much more than print inside loops, that's just a convenient way to demonstrate the inner workings.
[code]for(a=1,10,for(b=a,10,print(a" "b)))[/code]
[code]for(a=1,10,for(b=a,10,if(a%b==0,print(a" "b))))[/code]

I can also print for(n=1,10,for(b=1,10,print(<insert operator between n and b here>))

Here's a more complicated example of loops from an old .gp file of mine:
[code]Crocker(lim)={
my(part=[3*5*17*257*65537],npart,d);
if(lim<3,return([]));
for(n=5,lim-1,
d=divisors(2^(2^n)+1);
npart=[];
for(i=2,#d,npart=concat(npart,vector(#part,j,part[j]*d[i])));
part=npart
);
npart=[];
lim=2^(2^lim);
for(i=1,#part,
npart=concat(npart,vector((lim\part[i] - 1)>>4,j,part[i]*bitor(1,(j-1)<<4)))
);
npart
};
addhelp(Crocker, "Crocker(lim): Returns a vector of numbers shown by R. Crocker (1971) to be not of the form 2^a + 2^b + p for a < b and p prime.");[/code]

You can see that a lot can go on inside a loop!

Also note my use of addhelp. If you paste this program into gp and type ?Crocker it will give you my custom help message. It's useful when you can't remember the exact syntax of a command, or if you wrote the program so long ago you had no idea what it did (as in this case).

Within PARI, can you use trial division via looping?

[QUOTE=3.14159;222294]Within PARI, can you use trial division via looping?[/QUOTE]

Sure.

[code]trialDivide(n)={
forprime(p=2,sqrtint(n),
if(n%p == 0, return(p))
);
n \\ n seems to be prime
};[/code]

question which of these is fastest and why:


for(n=1,200,if(factor(n)!=[n 1],print(n" is prime")))
for(n=1,200,if(divisors(n)!=[1,n],print(n" is prime")))
for(n=1,200,if(isprime(n),print(n" is prime")))

[QUOTE=science_man_88;222302]question which of these is fastest and why:


for(n=1,200,if(factor(n)!=[n 1],print(n" is prime")))
for(n=1,200,if(divisors(n)!=[1,n],print(n" is prime")))
for(n=1,200,if(isprime(n),print(n" is prime")))[/QUOTE]

For possibly-large numbers, the third is vastly faster since the others involve factoring n which is difficult.

next question
 

according to Wikipedia there isn't a way to do arrays in Pari GP I think I know a way I just need to have a safe area in memory to point getstack() at as really arrays are just a line of memory addresses with the data inside them I learned this when i took a look at ASM. so is there a way you could help implement arrays ?

[QUOTE=science_man_88;224067]according to Wikipedia there isn't a way to do arrays in Pari GP[/QUOTE]

The usual way to do array stuff in gp itself is with vectors. To work with it like an array, just do

[code]v=vector(100);for(i=1,#v,v[i]=i^2)[/code]
which creates a vector then fills it (in this case, with the squares). You can get the length of a vector with the # operator, as in the above code.

There's also a shortcut:
[code]vector(100,n,n^2)[/code]

[QUOTE=science_man_88;224067]I think I know a way I just need to have a safe area in memory to point getstack() at as really arrays are just a line of memory addresses with the data inside them I learned this when i took a look at ASM. so is there a way you could help implement arrays ?[/QUOTE]

If you need a literal array, and vectors won't do, then you need to program in C [using the Pari library as needed]. In that case you can create them directly.

I was thinking something along:

array(n) = for(x=1,n,getstack())) but I don't know how do put it in the position of the variable it uses and I'd need a variable to keep the original place where I added it into also in ASM going down in a loop is more efficient (probably from numbers getting smaller), and i ASM I'd need to put the last one in first so I can read it first as it's last in = first out.

actually would array(n) = for(x=1,n,write(getstack(),data))) work ? the hard part is how to get it in all at once

What are you actually trying to do?

[QUOTE=CRGreathouse;224077]What are you actually trying to do?[/QUOTE]

my loop was basically to get write writing data into memory locations instead of what we think of as files. like I said the hard part I find is passing the array data without having an array already made.

and the fact write only writes to file lol.

though if we can read and write to file we can do file arrays or a array into a readable file to do something with.

[CODE]array(n) = for(x=1,n,write1("E:\\"getstack()".txt",hi))[/CODE]

if we can give the array a name with something we can write to E:\\array name.txt.

[QUOTE=science_man_88;224078]my loop was basically to get write writing data into memory locations instead of what we think of as files. like I said the hard part I find is passing the array data without having an array already made.[/QUOTE]

(This is frustrating.)

What are you actually trying to do? What data are you trying to put in memory?

getstack() is for debugging; it's not directly useful for memory access in gp.

[QUOTE=CRGreathouse;224084](This is frustrating.)

What are you actually trying to do? What data are you trying to put in memory?

getstack() is for debugging; it's not directly useful for memory access in gp.[/QUOTE]

I was just thinking of a way to get an array to work getstack only give a variable standing for a pointer to memory address I wanted the pointer to point to address to write to lol but alas this can't work.

[QUOTE=science_man_88;224094]I was just thinking of a way to get an array to work[/QUOTE]

So, use vectors.

can we make vectors and then search them for a number couldn't we then sieve through them and check if it's in them if not it's prime ?

[QUOTE=science_man_88]can we make vectors and then search them for a number couldn't we then sieve through them and check if it's in them if not it's prime ?
[/QUOTE]

factor(%) function.

so factor(x)==vector() is the only way ?

[QUOTE=science_man_88]so factor(x)==vector() is the only way ?
[/QUOTE]

Or you could individually factor each number that showed up in the vector.

[QUOTE=science_man_88;224102]can we make vectors and then search them for a number couldn't we then sieve through them and check if it's in them if not it's prime ?[/QUOTE]

Are you asking if you could write a sieve in Pari?

If so, yes. For example:

[code]sieveTo(n)={
v=vector(n,i,1); \\ Sets up a vector n long, all filled with 1s.
v[1]=0; \\ 1 isn't prime
forprime(p=2,sqrt(n),
forstep(i=2*p,n,p, \\ Mark 2p, 3p, ..., n as nonprime
v[i] = 0
)
);
for(i=1,#v,
if(v[i], print1(i" ")) \\ This is bad form, printing instead of just returning the vector... but I'll do it here so you can see what's going on.
);
};[/code]

[CODE]forprime(p=2,100,print1(vector(100,n,6*n*p+p)))[/CODE]

this is the start of a pari code for my idea using vector(s) then.

okay maybe not but really it compares:

vector(100,n,6*n*p+p)
vector(100,n,14*m+7)
and then the indexes that match up with should be either compared with 002450 or a less thick sequence of
A121290 which seems to be connected to other sequences but shows only 002450 that would give prime exponents to begin with.

if m and n weren't different I'd just use if (vector(100,n,14*m+7)-vector(100,n,6*n*p+p)==0,print(m)) type thing.

[QUOTE=science_man_88;224136][CODE]forprime(p=2,100,print1(vector(100,n,6*n*p+p)))[/CODE]

this is the start of a pari code for my idea using vector(s) then.[/QUOTE]

I don't know what you're trying to do, but this almost certainly isn't it.

[QUOTE=science_man_88;224136]okay maybe not but really it compares:

vector(100,n,6*n*p+p)
vector(100,n,14*m+7)
and then the indexes that match up with should be either compared with 002450 or a less thick sequence of
A121290 which seems to be connected to other sequences but shows only 002450 that would give prime exponents to begin with.[/QUOTE]

There's rarely a good reason to use the vector() command inside a loop. Make the vector and then manipulate the vector just like you would an array:
[code]v=vector(100);
for(i=1,#v,
v[i] = 1
);[/code]

[QUOTE=CRGreathouse;224142]I don't know what you're trying to do, but this almost certainly isn't it.[/QUOTE]

loop makes the vectors on the fly.

[QUOTE=CRGreathouse;224142]There's rarely a good reason to use the vector() command inside a loop. Make the vector and then manipulate the vector just like you would an array:
[code]v=vector(100);
for(i=1,#v,
v[i] = 1
);[/code][/QUOTE]

okay so how do I get one to exactly line up with the other so that my vector()-vector() method can find 0 at the right m and get m such that if can compare it to a formula linking it with a index in one of the 2 sequences I gave ?

[CODE]v=vector(1000);forprime(p=2,2000,for(i=1,#v,v[i] = 6*i*p+p;if(v[i]%24==7,f=v[i]\24;print(p"n+"f%p);break())))[/CODE]

this is the best i can do and it doesn't find my original 23n+16 example. I have no idea where I messed up.

[QUOTE=science_man_88;224145]okay so how do I get one to exactly line up with the other so that my vector()-vector() method can find 0 at the right m and get m such that if can compare it to a formula linking it with a index in one of the 2 sequences I gave ?[/QUOTE]

Perhaps you can explain yourself by giving examples of the arrays you'd like to construct (in a small case) and what you'd put in them. Use whatever language you like, or just write out arrays as lists of numbers. Then explain what you do with them in your small example.

[QUOTE=CRGreathouse;224185]Perhaps you can explain yourself by giving examples of the arrays you'd like to construct (in a small case) and what you'd put in them. Use whatever language you like, or just write out arrays as lists of numbers. Then explain what you do with them in your small example.[/QUOTE]


6*i*p+p is my vector for 6np+p and I was checking when they hit 24m+7 though i forgot to include 6np-p as well which is probably why I didn't get my result. take those m that they intersect 24m+7 at and check if it's in one of the 2 sequences I listed earlier about this in this thread.

[CODE] v=vector(1000);forprime(p=2,2000,for(i=1,#v,v[i] = 6*i*p+p;if(v[i]%24==7,f=v[i]\24;print(p"n+"f%p);break())));forprime(p=2,2000,for(i=1,#v,v[i] = 6*i*p-p;if(v[i]%24==7,f=v[i]\24;print(;break())))[/CODE]

this found 23n+16 now instead of printing them I'd like to check them against one of the sequences I talked of earlier.

Your code is doing something it shouldn't
[CODE]v=vector(1000);
forprime(p=2,2000,
for(i=1,#v,
v[i] = 6*i*p+p;
if(v[i]%24==7,
f=v[i]\24;
print(p"n+"f%p);
break()
)
)
);
forprime(p=2,2000,
for(i=1,#v,
v[i] = 6*i*p-p;
if(v[i]%24==7,
f=v[i]\24;
print([COLOR="Red"];break()[/COLOR])
)
)[/CODE]

In the future, please format your code so that these sorts of bugs (1) don't happen, or at least (2) are easier to notice.

[QUOTE=science_man_88;224188]this found 23n+16 now instead of printing them I'd like to check them against one of the sequences I talked of earlier.[/QUOTE]

This doesn't help me because I don't have any idea of what the code is supposed to do. As written it doesn't even really use the vector/array. Look:
[CODE]v=vector(1000);
forprime(p=2,2000,
for(i=1,#v,
[COLOR="Navy"]not_a_vector[/COLOR] = 6*i*p+p;
if([COLOR="Navy"]not_a_vector[/COLOR]%24==7,
f=[COLOR="Navy"]not_a_vector[/COLOR]\24;
print(p"n+"f%p);
break()
)
)
);
forprime(p=2,2000,
for(i=1,#v,
[COLOR="Navy"]not_a_vector[/COLOR] = 6*i*p-p;
if([COLOR="Navy"]not_a_vector[/COLOR]%24==7,
f=[COLOR="Navy"]not_a_vector[/COLOR]\24;
print([COLOR="Red"];break()[/COLOR])
)
)[/CODE]
It can be re-written do do just the same thing without actually ever using the vector (well, except to check its length, which could just be a constant instead). So there's really no hope that I could figure out what you want to do with vectors based on this code that isn't... you know... using vectors.

I see no reason it isn't a vector according to your words on vectors.

[QUOTE=science_man_88;224190]I see no reason it isn't a vector according to your words on vectors.[/QUOTE]

You say you want to use arrays (or vectors) but you post a program which doesn't actually appear to use them. (It looks like you just copied my example program and put your program inside it.) You don't explain what you're trying to do in any way I can understand.

Apparently, you can't understand the point I made in my last post, and I can't understand you at all. If you figure out what I was saying, find a different way to explain yourself, or get someone to explain to me what you're trying to do, great. Until then I don't think we can move forward.

I'm trying to use arrays to enact my idea about A002450 posted in the wheel factorization thread.

[QUOTE=science_man_88;224193]I'm trying to use arrays to enact my idea about A002450 posted in the wheel factorization thread.[/QUOTE]

Yes. But I didn't understand your idea there any more than I do here. And in particular I don't know what you actually want to do with the arrays.


Let me explain how my example works. Maybe then you can similarly explain what you have in mind.

I started with a list of, say, 12 numbers, all set to 1:
v=vector(12, i, 1)
%1 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

I then change the first to a 0, since 1 isn't prime:
v[1] = 0; v
%2 = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

I go through the list, marking every second element 0, since even numbers (other than 2) aren't prime:
forstep(i=4,12,2,v[i]=0);v
%3 = [0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]

Then I do the same for every third element, since numbers divisible by 3 (other than 3) aren't prime:
forstep(i=6,12,3,v[i]=0);v
%4 = [0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0]

The prime numbers are marked 1 in the vector, while the nonprimes are 0. I can print them if I want to:
for(i=1,12,if(v[i],print1(i" ")))
2 3 5 7 11


This uses an array. That is, I store values in different places in the vector and refer back to those locations. Your code used one cell at a time as scratch space, never again looking back at the change it made to the vector.

okay:

1) I stated 6np+/-p intersecting 24m+7 means for that m 24m+7 can't be prime,

2) I also stated that the m that given odd indexed Mersenne numbers is A002450 by the looks of it (A121290 seems to be the list for prime exponents)

3) I found patterns to suggest the m values take on m=p*x+c form.

so I want to find out when A002450 or A121290 fit this form for a calculated c(in more steps than needed probably, according to you I can calculated c based on p last I checked), and i wanted to try to use array form to figure it out like holding the sequence I want to check against etc.

[QUOTE=science_man_88;224198]1) I stated 6np+/-p intersecting 24m+7 means for that m 24m+7 can't be prime,[/QUOTE]

I *suspect* that this means something along the lines of: "If 24m + 7 is of the form 6np + p or 6np - p, then 24m + 7 is composite", where m and n are positive integers and p is a prime.

In that case, as discussed before, it's true, with 6n + 1 or 6n - 1 being one factor and p being another.

Here's a snippet of Pari code that checks if 24m + 7 is of the form [TEX]p(6n\pm1)[/TEX]:
[code]twentyFourMPlusSevenIsPTimesSixNPlusOrMinusOne(m)=!isprime(abs(24*m+7))[/code]

[QUOTE=science_man_88;224198]2) I also stated that the m that given odd indexed Mersenne numbers is A002450 by the looks of it (A121290 seems to be the list for prime exponents)[/QUOTE]

I don't understand anything you wrote here.

[QUOTE=science_man_88;224198]3) I found patterns to suggest the m values take on m=p*x+c form.[/QUOTE]

Which m values? The m values in 24m + 7, which can be any positive integer?

And, dare I ask, what patterns?

[QUOTE=science_man_88;224198]so I want to find out when A002450 or A121290 fit this form for a calculated c(in more steps than needed probably, according to you I can calculated c based on p last I checked), and i wanted to try to use array form to figure it out like holding the sequence I want to check against etc.[/QUOTE]

You say you want to use an array to figure "it" out, but I don't know what "it" is, nor how to find it, nor whether it actually has anything to do with arrays.

@CRG: A lack of definition for many terms.. I don't see the appeal to the ambiguity, but I think, to some extent, that it's just me. What is it that sm88 is attempting to do concerning PARI's vectors? By the way, searching for a prime greater than 100k digits is a trouble not worth going through. Would have taken me 17 years to find something in the range I was looking for.

6np+/-p==24m+7

m=p*x+ c(start value of m that work for the p value)

take p*x+c and figure out what it hits in either A002450 or A121290 eliminate the index it hits. continue sieve.

@sm88: I wish you luck.

could we use vectors to use do vector(x,n,n)-vector() stuff to eliminate multiples of primes greater than their squares ?

[QUOTE=science_man_88]could we use vectors to use do vector(x,n,n)-vector() stuff to eliminate multiples of primes greater than their squares ?
[/QUOTE]

A sieve already eliminates every multiple of p, regardless of whether or not n> p[sup]2[/sup]

[QUOTE=science_man_88;224547]could we use vectors to use do vector(x,n,n)-vector() stuff to eliminate multiples of primes greater than their squares ?[/QUOTE]

I don't know what you want to do, or whether vectors would be a useful way to do it. But there are certainly many different ways to sieve. If you wanted to find all squarefree numbers up to x, for example, you could simply remove multiples of p^2 where p runs from 2 to sqrt(N).

[QUOTE=CRGreathouse;224552]I don't know what you want to do, or whether vectors would be a useful way to do it. But there are certainly many different ways to sieve. If you wanted to find all squarefree numbers up to x, for example, you could simply remove multiples of p^2 where p runs from 2 to sqrt(N).[/QUOTE]

i basically want to try something like the sieve of Eratosthenes with vectors to eliminate multiples greater than p^2

[QUOTE=science_man_88;224554]i basically want to try something like the sieve of Eratosthenes with vectors to eliminate multiples greater than p^2[/QUOTE]

So you have a list of numbers from 1 to N and you want to remove multiples of p^2 for what range of p? (I assume you mean for p to be prime.)

[QUOTE=CRGreathouse;224555]So you have a list of numbers from 1 to N and you want to remove multiples of p^2 for what range of p? (I assume you mean for p to be prime.)[/QUOTE]

not multiples of p^2 multiples of p > p^2

[QUOTE=science_man_88;224556]not multiples of p^2 multiples of p > p^2[/QUOTE]

Fine. Which p?

[QUOTE=CRGreathouse;224557]Fine. Which p?[/QUOTE]

well i was hoping for a script that could be applied regardless of prime p but up to you I guess.

once we delete multiples of 2 is there a way to go through and say go to 3 and delete all those multiples left. etc.

[QUOTE=science_man_88;224558]well i was hoping for a script that could be applied regardless of prime p but up to you I guess.[/QUOTE]

You only want it to remove *one* value of p?

[code]doit(n,p)={
v=vector(n,i,1);
forstep(k=p^2+p,n,p,
v[k]=0
);
v
};[/code]

[CODE]forprime(x=2,10,print(vector(100,n,n)%x))[/CODE]

is the best i can do but it doesn't keep the last ones eliminated.

[QUOTE=science_man_88]once we delete multiples of 2 is there a way to go through and say go to 3 and delete all those multiples left. etc.
[/QUOTE]

Please, see [URL="http://en.wikipedia.org/wiki/Sieve_of_eratosthenes"]here[/URL].

sm88, once again I tell you: don't put vector() commands inside a loop. It's a sign that you're doing something wrong.

Rule of thumb: create a vector outside a loop, then access its members inside a loop with the []s.

[QUOTE=CRGreathouse;224562]sm88, once again I tell you: don't put vector() commands inside a loop. It's a sign that you're doing something wrong.

Rule of thumb: create a vector outside a loop, then access its members inside a loop with the []s.[/QUOTE]
then i have no idea how to use crap all in pari which means I'm useless also it does knock them out just not forever and I have no idea how to print the result of eliminating all sequences with yours. and when I do it on it's own I can't understand it's result.

[QUOTE=science_man_88;224563]then i have no idea how to use crap all in pari which means I'm useless also it does knock them out just not forever and I have no idea how to print the result of eliminating all sequences with yours. and when I do it on it's own I can't understand it's result.[/QUOTE]

I suspect the underlying problem is that you don't understand how to use arrays, rather than having trouble with Pari specifically.

[QUOTE=CRGreathouse;224565]I suspect the underlying problem is that you don't understand how to use arrays, rather than having trouble with Pari specifically.[/QUOTE]

I just don't see why using a vector inside a loop when to check all results one after another it seems necessary is a sin.

I found my struggle with your code and I know why it's my fault but I still can't figure out how to print the results(indexes) that don't contain 0 afterwards.

[QUOTE=science_man_88;224567]I just don't see why using a vector inside a loop when to check all results one after another it seems necessary is a sin.[/QUOTE]

You wrote here that your code (with vectors initialized in a loop) "doesn't keep" the results:

[QUOTE=science_man_88;224560]is the best i can do but it doesn't keep the last ones eliminated.[/QUOTE]

This is because, rather than keep the results, you [i]overwrite them at every step[/i]. I'm telling you: don't overwrite them!

If it bothers you that I tell you to not do this, fine -- feel free to write code however you like. But don't be surprised that the results are eliminated when [i]you[/i] are eliminating them.

[QUOTE=science_man_88;224568]I found my struggle with your code and I know why it's my fault but I still can't figure out how to print the results(indexes) that don't contain 0 afterwards.[/QUOTE]

If you have a vector v and you want to print the indexes where the vector is nonzero, just loop through the entries of the vector and print when, uh, the value is nonzero:
[code]for(i=1,#v,if(v[i] != 0, print(i)))[/code]

[QUOTE=CRGreathouse;224571]You wrote here that your code (with vectors initialized in a loop) "doesn't keep" the results:



This is because, rather than keep the results, you [i]overwrite them at every step[/i]. I'm telling you: don't overwrite them!

If it bothers you that I tell you to not do this, fine -- feel free to write code however you like. But don't be surprised that the results are eliminated when [i]you[/i] are eliminating them.[/QUOTE]

then how can i eliminate all multiples>p^2 of multiple p values without
a loop ?

[QUOTE=science_man_88;224573]then how can i eliminate all multiples>p^2 of multiple p values without
a loop ?[/QUOTE]

I didn't say not to use loops, I said not to create vectors in loops. Create a vector then loop through the primes, removing multiples as needed from the vector (that you create only once).

Example:
[code]doit(n)={
v=vector(n,i,1);
forprime(p=2,sqrt(n),
forstep(k=p^2+p,n,p,
v[k]=0
)
);
v
};[/code]

[QUOTE=CRGreathouse;224574]I didn't say not to use loops, I said not to create vectors in loops. Create a vector then loop through the primes, removing multiples as needed from the vector (that you create only once).

Example:
[code]doit(n)={
v=vector(n,i,1);
forprime(p=2,sqrt(n),
forstep(k=p^2+p,n,p,
v[k]=0
)
);
v
};[/code][/QUOTE]


okay and this code when I do print(v) after it's in give me back every number.
it's seems only to eliminate them when i use a forprime loop but then it's back to eliminating results only to overwrite.

[QUOTE=science_man_88;224575]okay and this code when I do print(v) after it's in give me back every number.
it's seems only to eliminate them when i use a forprime loop but then it's back to eliminating results only to overwrite.[/QUOTE]

The type of loop I use (forprime, for, forstep, fordiv...) doesn't matter in terms of overwriting.

This code doesn't overwrite the vector once it's created. The only time that I write "v=" is at the beginning.

[QUOTE=CRGreathouse;224577]The type of loop I use (forprime, for, forstep, fordiv...) doesn't matter in terms of overwriting.

This code doesn't overwrite the vector once it's created. The only time that I write "v=" is at the beginning.[/QUOTE]

got it working for what i wanted except now i want to get rid of all the 0's

[CODE]doit(n){=v=vector(n,i,i);forprime(p=1,sqrt(n),forstep(k=p^2,n,p,v[k]=0);v});[/CODE]

oh and the number one.

the forprime loop limits how high I can go though, wish I could use v as primelist then I could square each time and go higher.

[QUOTE=science_man_88;224578]got it working for what i wanted except now i want to get rid of all the 0's[/QUOTE]

It won't surprise you that I don't know what you mean. The code you posted is what I posted, except that
1. Instead of storing 1 for every number, you store the value of the number at each number; if you work with big numbers this will waste space.
2. Rather than return the vector when the function finishes, you 'return' it as the output of each forprime() step, which discards the result.
3. You remove multiples of p >= p^2 rather than multiples of p > p^2.
4. You removed the formatting that made my code easy to read.

[QUOTE=CRGreathouse;224580]It won't surprise you that I don't know what you mean. The code you posted is what I posted, except that
1. Instead of storing 1 for every number, you store the value of the number at each number; if you work with big numbers this will waste space.
2. Rather than return the vector when the function finishes, you 'return' it as the output of each forprime() step, which discards the result.
3. You remove multiples of p >= p^2 rather than multiples of p > p^2.
4. You removed the formatting that made my code easy to read.[/QUOTE]

actually I think you'll find trying to alter it after first declaring it left it in this line, okay but how do i know which number it is then ?
if the result is known why keep it in memory why not print it ? (though I don't see how either is relevant until I invoke print(v) after entering this again)

yeah I realized p^2 wasn't prime and that's what caused me to think yours didn't work.

If he wants a sieve: He should program one! Programming a sieve is kiddie work! Why so much confusion?

Also: I made a π(x) approximation function, although it only works when n > 10[sup]30[/sup]

[QUOTE=science_man_88;224582]actually I think you'll find trying to alter it after first declaring it left it in this line, okay but how do i know which number it is then ?[/QUOTE]

I'm altering individual members of the vector (by setting them equal to 0, in this case), but not replacing the entire array. You had been creating entirely new arrays and either overwriting the old one or never storing it in the first place.

I don't know what you mean by "how do i know which number it is then".

[QUOTE=science_man_88;224582]if the result is known why keep it in memory why not print it ?[/QUOTE]

That's a good general question. You want to keep it in memory so that you can do things with it later. For example, say you make a program that prints twin primes. You print the twin primes up to 100... great. Now you want the twin primes up to a million, but why put them on the screen? You won't even be able to read them all! Instead you store them in memory and do something else with them.

Even when the output is small enough to conveniently read, you want to keep it in memory when possible so that you can do things with it.

[QUOTE=science_man_88;224582]yeah I realized p^2 wasn't prime and that's what caused me to think yours didn't work.[/QUOTE]

I had the program do what you told me you wanted it to do. You didn't tell me the purpose -- maybe you were trying to compute [url=http://oeis.org/classic/A000430]A000430[/url].

Of course since you were just generating the primes, I could have told you to not bother with generating them and to use the built-in commands instead. :smile:

[QUOTE=3.14159;224583]If he wants a sieve: He should program one! Programming a sieve is kiddie work! Why so much confusion?[/QUOTE]

Your question, "Why so much confusion?", is an interesting and difficult question. I'm not an educator, though, and don't have an answer for you.

[QUOTE=3.14159;224583]Also: I made a π(x) approximation function, although it only works when n > 10[sup]30[/sup][/QUOTE]

Can I see? I have one as well.

[QUOTE=CRGreathouse]Your question, "Why so much confusion?", is an interesting and difficult question. I'm not an educator, though, and don't have an answer for you.
[/QUOTE]

I was curious how there would be confusion over making a simple sieve.

[QUOTE=CRGreathouse]Can I see? I have one as well.
[/QUOTE]

px(n) = if(n>10^30, return(floor(n/log(n))).

(I don't see why they eliminated ln.)

[QUOTE=3.14159;224587]I was curious how there would be confusion over making a simple sieve.[/QUOTE]

I take it that you're not an educator, either. :smile: It's hard to teach most anything, and harder still to see what parts of things will be confusing.

[QUOTE=3.14159;224587]px(n) = if(n>10^30, return(floor(n/log(n))).[/QUOTE]

Mine's better. :razz: I'll snip out the part where it gives sources, since I haven't updated that part to account for the newest paper it uses.

[code]piBounds(x)={
my(lower,upper,lx,t,lowerRH,upperRH,roundFlag=default(realprecision)>sizedigit(x));

if(roundFlag, x = floor(x));
if (x < default(primelimit),
lower=primepi(x);
print("primepi("x") = "lower);
return([lower,lower]);
);

lx = log(x);
lower = x/lx * (1 + 1/lx + 1.8/lx^2); \\ Dusart 1998, x >= 32299
upper = x/(lx - 1.1); \\ Dusart 1998, x >= 60184

if (x >= 88783,
lower = x/lx * (1 + 1/lx + 2/lx^2); \\ Dusart 2010, x >= 88783
);

if (x >= 355991,
t = x/lx * (1 + 1/lx + 2.51/lx^2); \\ Dusart 1998, x >= 355991
if (upper > t, upper = t)
);

if (x >= 2953652302,
t = x/lx * (1 + 1/lx + 2.334/lx^2); \\ Dusart 2010, x >= 2 953 652 302
if (upper > t, upper = t)
);

if (x >= 1.332e10,
t = x/lx * (1 + 1.0992/lx); \\ Dusart 1998, x >= 1.332e10
if (upper > t, upper = t)
);

if (roundFlag,
lower = ceil(lower);
upper = floor(upper);
);

t = sqrt(x)/8/Pi*log(x);
lowerRH = li(x)-t; \\ Schoenfield, x >= 2657
upperRH = li(x)+t; \\ Schoenfield, x >= 2657
if (roundFlag,
lowerRH = ceil(lowerRH);
upperRH = floor(upperRH);
);

print("For primepi("x"):");
print(lower" (lower bound)");
if (lowerRH > lower,
print(lowerRH" (lower bound under the RH)");
);
if(roundFlag,
print(round(R(x))" (Riemann R, approximate)");
print(round(Li(x))" (logarithmic integral, apx)");
,
print(R(x)" (Riemann R, approximate)");
print(Li(x)" (logarithmic integral, apx)");
);
if (upperRH < upper,
print(upperRH" (upper bound under the RH)");
);
print(upper" (upper bound)");

[lower, upper]
};
addhelp(piBounds, "piBounds(x, verbose=0): Bounds on primepi(x).");


li(x)=-eint1(-log(x));
addhelp(li, "Logarithmic integral.");


Li(x)=eint1(-2) - eint1(-log(x));
addhelp(Li, "Offset logarithmic integral, an estimate for primepi. Crandall and Pomerance call this li_0.");


\\Euler + log(x) + suminf(n=1,x^n/(n*n!))
\\Euler + x + suminf(n=1,x^n/(n*n!))
Ei(x)={
-eint1(-x)
};
addhelp(Ei, "Exponential integral.");[/code]

[QUOTE=CRGreathouse]I take it that you're not an educator, either. :smile: It's hard to teach most anything, and harder still to see what parts of things will be confusing.
[/QUOTE]

Far from it, to be exact. I's best to learn on one's own.

[QUOTE=CRGreathouse]Mine's better. :razz: I'll snip out the part where it gives sources, since I haven't updated that part to account for the newest paper it uses.
[/QUOTE]

:orly owl: Is this a conspiracy to [URL="http://www.mersenneforum.org/showpost.php?p=223602&postcount=44"]steal my work[/URL]?

[code]>px(1e31)
%1 = 140094994162339299242299651263[/code]

[code]>piBounds(1e31)
For primepi(10000000000000000000000000000000):
142112646690226626210386971421 (lower bound)
142115097348071950748902528525 (lower bound under the RH)
142115097348080886394439772958 (Riemann R, approximate)
142115097348080932014151407887 (logarithmic integral, apx)
142115097348089913279400287260 (upper bound under the RH)
142121830318981733579096747865 (upper bound)
%2 = [142112646690226626210386971421, 142121830318981733579096747865][/code]

So yours is a little light. :smile:

[QUOTE=CRGreathouse]So yours is a little light. :smile:
[/QUOTE]

[B]Coughcoughbutminedoesnotdependonunprovenspeculations/conjectures,youseecoughcough.
Coughcoughseeingasyouhavezeroproofsoftheconjecturesthatyoudependon,Iwouldn'tbesoquicktomakesmugstatementscoughcough.[/B]

[QUOTE=3.14159;224599]Coughcoughbutminedoesnotdependonunprovenspeculations/conjectures,youseecoughcough.[/QUOTE]

Wait... kettle, pot?

My results rely on no unproven conjectures. The results returned are exact; further, additional information is printed (like R(x)) that provide guesses that are probably better than the bounds themselves. In particular, my program does [b]not[/b] depend on the Riemann Hypothesis. (This shouldn't be surprising, considering that the key paper on which the program is based is called "Estimates of Some Functions Over Primes without R.H.".)

Yours provides (as far as I can tell) only a guess, not stating whether it's an upper bound, a lower bound, or whether it can be both low and high at different points.

new code problem
 

[QUOTE]Originally Posted by science_man_88
Code:
[CODE]v=vector(100);for(i=1,#v,v[i]=1;if(i%5==4,v[i]==0))[/CODE]
I'm trying to eliminate obvious patterns to match isprime(6n+1) and maybe i'll try 6*n-1 after that.

or is this all trivial. probably
That code actually works. Post it on the forums for comment,[SPOILER] though, not to PM -- I have too many PMs to keep track of, they're hard to review, and I'm running out of space in my inbox.[/SPOILER][/QUOTE]

like I said I'm trying to eliminate obvious(even to me) patterns from 6n+1 without doing math on the numbers themselves.

the one's that seem to work so far are 5x+4,7x+8, and 10x+9, not sure if it's the law of small numbers but these all have a difference 1 between the 2 coefficients.

oh wait 10x+9 is a off branch of 5x+4 doh.

Your code marks the indexes that are 4 mod 5 as 0 and the others as 1. To do this it requires one division per vector element.

[QUOTE=science_man_88;224612]like I said I'm trying to eliminate obvious(even to me) patterns from 6n+1 without doing math on the numbers themselves.[/QUOTE]

I'm not sure what patterns you're looking for, or what you mean by 'without doing math on the numbers themselves'. But your code runs and doesn't do anything obviously counterproductive like overwriting arrays.

[QUOTE=CRGreathouse;224619]Your code marks the indexes that are 4 mod 5 as 0 and the others as 1. To do this it requires one division per vector element.



I'm not sure what patterns you're looking for, or what you mean by 'without doing math on the numbers themselves'. But your code runs and doesn't do anything obviously counterproductive like overwriting arrays.[/QUOTE]

I could replace division with a whole lot of subtraction but I doubt that's better.

the reason I'm using patterns is I seem to have found some in:

[CODE](14:33) gp > vector(100,n,isprime(6*n+1))
%75 = [1, 1, 1, [COLOR="Red"]0[/COLOR], 1, 1, 1, [COLOR="Lime"]0[/COLOR], [COLOR="red"]0[/COLOR], 1, 1, 1, 1,[COLOR="red"] 0[/COLOR], [COLOR="lime"]0[/COLOR], 1, 1, 1, [COLOR="red"]0[/COLOR], 0, 1, [COLOR="lime"]0[/COLOR], 1, [COLOR="red"]0[/COLOR], 1, 1, 1, 0, [COLOR="red"]0[/COLOR], 1, 0, 1, 1, [COLOR="red"]0[/COLOR], 1, 0, 1, 1, [COLOR="red"]0[/COLOR], 1, 0, 0, 0, [COLOR="red"]0[/COLOR], 1, 1, 1, 0, [COLOR="red"]0[/COLOR], 0, 1, 1, 0, [COLOR="red"]0[/COLOR], 1, 1, 0, 1, [COLOR="red"]0[/COLOR], 0, 1, 1, 1, [COLOR="red"]0[/COLOR], 0, 1, 0, 1, [COLOR="red"]0[/COLOR], 1, 0, 1, 1, [COLOR="red"]0[/COLOR], 0, 1, 1, 0, [COLOR="red"]0[/COLOR], 0, 1, 0, 1, [COLOR="red"]0[/COLOR], 0, 0, 1, 0, [COLOR="red"]0[/COLOR], 1, 1, 0, 0, [COLOR="red"]0[/COLOR], 1, 1, 0, 0, [COLOR="red"]0[/COLOR], 1][/CODE]

legend: [COLOR="Red"]#[/COLOR] = 5x+4 indexes [COLOR="Lime"]#[/COLOR] = 8x+7 indexes (not all filled in)

I think I find 8x+20 as well which bust my pattern of patterns lol.

the ones are supposed to represent numbers in vector(100,n,6*n+1) in this case.

Here are your indices with composites of the form [url=http://factordb.com/search.php?query=6*n%2B1&v=n&n=1&EC=1&E=1&C=1&FF=1&CF=1&of=H&pp=50&sw=Update]6n+1[/url]

[QUOTE=kar_bon;224622]Here are your indices with composites of the form [url=http://factordb.com/search.php?query=6*n%2B1&v=n&n=1&EC=1&E=1&C=1&FF=1&CF=1&of=H&pp=50&sw=Update]6n+1[/url][/QUOTE]

thanks maybe I can find a few patterns to help eliminate them then apply to a higher range.

is this trying to tell me "we don't care" seems like it for some reason lol (maybe it's just because I find you helpful lol)

[QUOTE=science_man_88;224621]I could replace division with a whole lot of subtraction but I doubt that's better.[/QUOTE]

Better would be to start at the first number that is 4 mod 4 and increase the index by 5 at a time, so instead of 5 divisions for 5 numbers you do 1 addition for 5 numbers.

[QUOTE=CRGreathouse;224625]Better would be to start at the first number that is 4 mod 4 and increase the index by 5 at a time, so instead of 5 divisions for 5 numbers you do 1 addition for 5 numbers.[/QUOTE]
wouldn't 4 mod 4 give a first of 0 if it has to be above that then it's 4 I don't have to divide then i can just turn them 0. the hard part is applying it to all patterns found care to elaborate more ?

since i know the patterns will likely intersect multiple times i can put in a check to see if it's already knocked out by checking if v[i] == 1 if it is continue if not break and try the next one.

anymore optimizations ?

[QUOTE=science_man_88;224626]wouldn't 4 mod 4 give a first of 0[/QUOTE]

It's a typo for 4 mod 5. 4 is right next to 5 on my keyboard.

[QUOTE=science_man_88;224626]since i know the patterns will likely intersect multiple times i can put in a check to see if it's already knocked out by checking if v[i] == 1 if it is continue if not break and try the next one.[/QUOTE]

If you do it my way you spend less time by just marking it 0 than you would spend testing.

[code]foo(n)={
v=vector(n,unused,1);
forstep(i=4,n,5,v[i]=0); \\ remove 4 mod 5
forstep(i=2,n,7,v[i]=0); \\ remove 2 mod 7
v
};

foo1(n)={
v=vector(n,unused,1);
for(i=1,n,if(i%5==4,v[i]=0)); \\ remove 4 mod 5
for(i=1,n,if(i%7==2,v[i]=0)); \\ remove 2 mod 7
v
};

foo2(n)={
v=vector(n,unused,1);
for(i=1,n,if(i%5==4,v[i]=0)); \\ remove 4 mod 5
for(i=1,n,if(v[i] && i%7==2,v[i]=0)); \\ remove 2 mod 7
v
};[/code]

foo is the way I'd write it (well... not really, but close enough). foo1 is how you might have written it, and foo2 is your suggestion on how to improve on foo1 (by not checking numbers that are already eliminated).

foo1(10^7) takes 13.5 seconds. foo2, your improvement, takes 12.8 seconds. foo takes 6.3 seconds.

[QUOTE=CRGreathouse;224627]It's a typo for 4 mod 5. 4 is right next to 5 on my keyboard.



If you do it my way you spend less time by just marking it 0 than you would spend testing.

[code]foo(n)={
v=vector(n,unused,1);
forstep(i=4,n,5,v[i]=0); \\ remove 4 mod 5
forstep(i=2,n,7,v[i]=0); \\ remove 2 mod 7
v
};

foo1(n)={
v=vector(n,unused,1);
for(i=1,n,if(i%5==4,v[i]=0)); \\ remove 4 mod 5
for(i=1,n,if(i%7==2,v[i]=0)); \\ remove 2 mod 7
v
};

foo2(n)={
v=vector(n,unused,1);
for(i=1,n,if(i%5==4,v[i]=0)); \\ remove 4 mod 5
for(i=1,n,if(v[i] && i%7==2,v[i]=0)); \\ remove 2 mod 7
v
};[/code]

foo is the way I'd write it (well... not really, but close enough). foo1 is how you might have written it, and foo2 is your suggestion on how to improve on foo1 (by not checking numbers that are already eliminated).

foo1(10^7) takes 13.5 seconds. foo2, your improvement, takes 12.8 seconds. foo takes 6.3 seconds.[/QUOTE]

would checking for ones knocked out slow foo(n) down ?.

when did 2 mod 7 come in ? check [url]http://www.research.att.com/~njas/sequences/A046954[/url]

[CODE]foo(n)={
v=vector(n,unused,1);
forstep(i=4,n,5,v[i]=0);
forstep(i=8,n,7,v[i]=0);
v
};[/CODE]

I'm guessing you meant

[QUOTE=science_man_88;224628]would checking for ones knocked out slow foo(n) down ?.[/QUOTE]

Yes, probably by at least 20% on modern deeply pipelined CPUs.

[QUOTE=science_man_88;224628]when did 2 mod 7 come in ?[/QUOTE]

I'm showing you how to sieve. I don't know what significance numbers that are 4 mod 5 have to you; I just picked a different residue class mod a small prime so I could show you the timing difference. (Actually, it's very surprising that they're so close; usually you'd expect the sieve to be 20 times faster, not just 2 times faster.)

Obviously you replace my choices of residue class with those relevant to the problem you're working on. I don't know what primes you want to work with, nor what congruence classes you want to remove mod those primes.

[QUOTE=CRGreathouse;224633]Yes, probably by at least 20% on modern deeply pipelined CPUs.



I'm showing you how to sieve. I don't know what significance numbers that are 4 mod 5 have to you; I just picked a different residue class mod a small prime so I could show you the timing difference. (Actually, it's very surprising that they're so close; usually you'd expect the sieve to be 20 times faster, not just 2 times faster.)

Obviously you replace my choices of residue class with those relevant to the problem you're working on. I don't know what primes you want to work with, nor what congruence classes you want to remove mod those primes.[/QUOTE]

for the timing yeah especially with me involved lol :) I'm sieving 6*n+1 for this one then i can apply the answer to 6n-1 with alteration.

with the [CODE]if(v[i]==0,break())[/CODE] alteration I tested to [CODE]foo(100000) [/CODE]for 5x+4 and 7x+8 and got 47 ms can't test any higher until I use [CODE]allocatemem()[/CODE]. just need more patterns to eliminate from it now lol

[QUOTE=science_man_88;224635]with the [CODE]if(v[i]==0,break())[/CODE][/QUOTE]

That's not the right thing to do. That means that when you hit the first 0 you stop doing anything with that prime, which will leave bad numbers in your vector.

[QUOTE=CRGreathouse;224637]That's not the right thing to do. That means that when you hit the first 0 you stop doing anything with that prime, which will leave bad numbers in your vector.[/QUOTE]

[CODE]foo(n)=v=vector(n,unused,1);forstep(i=4,n,5,if(v[i]==0,break());v[i]=0);forstep(i=8,n,7,if(v[i]==0,break());v[i]=0);v; [/CODE]


I haven't found any I'll make a code to check but I haven't found any.

[CODE]for(i=4,100000,if(i%5==4 || i%7==8,if(v[i]!=0,print(i))))[/CODE]

my checking code, any critics ?

and it didn't find anything.

[QUOTE=science_man_88;224638][CODE]foo(n)=v=vector(n,unused,1);forstep(i=4,n,5,if(v[i]==0,break());v[i]=0);forstep(i=8,n,7,if(v[i]==0,break());v[i]=0);v; [/CODE]


I haven't found any I'll make a code to check but I haven't found any.[/QUOTE]

It looks like your code is trying to remove numbers that are 4 mod 5 or 1 mod 7. The break in the second forstep means that the second loop stops when it reaches 29, so 36, 43, ... are never marked off.

[QUOTE=science_man_88;224638][CODE]for(i=4,100000,if(i%5==4 || i%7==8,if(v[i]!=0,print(i))))[/CODE]

my checking code, any critics ?

and it didn't find anything.[/QUOTE]

Considering that n%k gives a number in the range [0, k-1], i%7 == 8 is always false.

[QUOTE=CRGreathouse;224639]It looks like your code is trying to remove numbers that are 4 mod 5 or 1 mod 7. The break in the second forstep means that the second loop stops when it reaches 29, so 36, 43, ... are never marked off.[/QUOTE]

so how do i handle it though my check code found nothing.

[QUOTE=science_man_88;224641]so how do i handle it though my check code found nothing.[/QUOTE]

You could write Mod(i,7) == 8. I would just write i%7 == 1, since [TEX]1\equiv8\pmod7[/TEX].

You'll find that the exceptions are what I told you they would be in post #93.

[QUOTE=CRGreathouse;224642]You could write Mod(i,7) == 8. I would just write i%7 == 1, since [TEX]1\equiv8\pmod7[/TEX].

You'll find that the exceptions are what I told you they would be in post #93.[/QUOTE]

first exception = 1 then

how to fix it ?