mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   PARI/GP (https://www.mersenneforum.org/forumdisplay.php?f=155)
-   -   PARI's commands (https://www.mersenneforum.org/showthread.php?t=13636)

science_man_88 2010-12-14 23:01

[QUOTE=axn;241856]This is decidedly _not_ Sieve of Eratosthenes. This is called "trial division". And a slow version at that. You'll get good reference pseudocode for SoE on the web. Try to implement one of those.[/QUOTE]

[CODE]sieve_of_eratosthenes(x) = d=vector((x-1),n,n+1);index=1;p=2;for(y=1,#d,if(d[y]!=0,index=y;p=d[y]);if(p^2>x,break(2));forstep(i=index+p,#d,p,if(d[i]!=0,d[i]=0)))[/CODE] this better ?

The main difference to what I did here and the other is that the numbers were represented poorly in the indexes of a vector.

science_man_88 2010-12-14 23:28

[QUOTE=science_man_88;241863][CODE]sieve_of_eratosthenes(x) = d=vector((x-1),n,n+1);index=1;p=2;for(y=1,#d,if(d[y]!=0,index=y;p=d[y]);if(p^2>x,break(2));forstep(i=index+p,#d,p,if(d[i]!=0,d[i]=0)))[/CODE] this better ?

The main difference to what I did here and the other is that the numbers were represented poorly in the indexes of a vector.[/QUOTE]

okay this kinda takes advantage of Euler's Sieve, never mind it just has the same idea lol.

science_man_88 2010-12-14 23:57

not sure how to turn it into Atkin or Sundaram

CRGreathouse 2010-12-15 00:36

[QUOTE=science_man_88;241844]I named my code to show you:[/QUOTE]

How about
[code]primesUPTO(x)=primes(primepi(x));[/code]
?

Also: I'm proud of you for making the function return a vector instead of just printing values. You're learning!

science_man_88 2010-12-15 01:09

[QUOTE=CRGreathouse;241880]How about
[code]primesUPTO(x)=primes(primepi(x));[/code]
?

Also: I'm proud of you for making the function return a vector instead of just printing values. You're learning![/QUOTE]

You're too good!!!

does your semiprime detector have anything to do with:

[CODE]fa = factor(N);
isprime( fa[,1] )[/CODE] which is listed in the PARI FAQ lol.

CRGreathouse 2010-12-15 01:42

[QUOTE=science_man_88;241887]You're too good!!!

does your semiprime detector have anything to do with:

[CODE]fa = factor(N);
isprime( fa[,1] )[/CODE] which is listed in the PARI FAQ lol.[/QUOTE]

No.

Right now it looks like this:
[code]long
issemiprime(GEN n)
{
if (typ(n) != t_INT)
pari_err(arither1, "issemiprime");
if (signe(n) <= 0)
return 0;

ulong nn = itou_or_0(n);
if (nn)
return uissemiprime(nn);

pari_sp ltop = avma;
if (!mpodd(n)) {
long ret = mod4(n) && isprime(shifti(n, -1));
avma = ltop;
return ret;
}

long p = 0;
byteptr primepointer = diffptr;
NEXT_PRIME_VIADIFF(p, primepointer); // Skip 2
for (;;)
{
NEXT_PRIME_VIADIFF(p, primepointer);
if (p > 997) // What is a good breakpoint here?
break;
if (dvdis(n, p)) // "1 if p divides n and 0 otherwise"
{
long ret = isprime(diviuexact(n, p));
avma = ltop;
return ret;
}
}

if (isprime(n))
return 0;

GEN fac = Z_factor_until(n, shifti(n, -1)); // Find a nontrivial factor -- returns just the factored part
GEN expo = gel(fac, 2);
GEN pr = gel(fac, 1);
long len = glength(expo);
if (len > 2) {
avma = ltop;
return 0;
}
if (len == 2) {
if (cmpis(gel(expo, 1), 1) > 0 || cmpis(gel(expo, 2), 1) > 0) {
avma = ltop;
return 0;
}
GEN p = gel(pr, 1);
GEN q = gel(pr, 2);
long ret = !cmpii(mulii(p, q), n) && isprime(p) && isprime(q);
avma = ltop;
return ret;
}
if (len == 1) {
long e = itos(gel(expo, 1));
if (e != 2) {
avma = ltop;
return 0;
}
GEN p = gel(pr, 1);
long ret = !cmpii(sqri(p), n) && isprime(p);
avma = ltop;
return ret;
}

pari_err(talker, "Z_factor_until returned an unexpected value %Ps at n = %Ps, exiting...", fac, n);
avma = ltop;
return NEVER_USED;
}[/code]

but I'm not done with it. Note that this is Pari, not GP.

science_man_88 2010-12-15 15:20

vector basics ?
 
[CODE](11:17)>v=vector(100,n,20)
%359 = [20, 20, 20, 20, 20, 20,
(11:18)>v=vector(100,n,n)
%360 = [1, 2, 3, 4, 5, 6, 7, 8,
(11:18)>v=vector(#v-1,n,v[n])
%361 = [1, 2, 3, 4, 5, 6, 7, 8,
(11:18)>v=vector(#v-1,n,v[n+1])
%362 = [2, 3, 4, 5, 6, 7, 8, 9,[/CODE]
[CODE](11:22)>vecsort(v)
%377 = [16, 17, 18, 19[/CODE]

science_man_88 2010-12-15 15:41

nice I found a faster way to reverse a vector lol, might have to update reverse with this one lol:

[CODE]v=vector(#v,n,v[#v-(n-1)])[/CODE]

science_man_88 2010-12-15 15:51

[QUOTE=science_man_88;241973]nice I found a faster way to reverse a vector lol, might have to update reverse with this one lol:

[CODE]v=vector(#v,n,v[#v-(n-1)])[/CODE][/QUOTE]

reversal of a Ascii string made easier:

[CODE]c=Vecsmall("hello");c=vector(#c,n,c[#c-(n-1)]);print(Strchr(c))[/CODE]

axn 2010-12-15 15:54

[QUOTE=science_man_88;241968][CODE]
(11:18)>v=vector(#v-1,n,v[n])
%361 = [1, 2, 3, 4, 5, 6, 7, 8,
(11:18)>v=vector(#v-1,n,v[n+1])
%362 = [2, 3, 4, 5, 6, 7, 8, 9,[/CODE]
[/QUOTE]

vecextract might possibly be better for the these two.

science_man_88 2010-12-15 16:01

[QUOTE=axn;241975]vecextract might possibly be better for the these two.[/QUOTE]

only thing I've been able to accomplish is throwing out every index I wanted to use and left with the one I wanted off.

never mind I got one half of those 2 to work but how do i do the other side ?

[code]c=vecextract(c,vector(#c-1,n,n))[/code]
[code]c=vecextract(c,vector(#c-1,n,n+1))[/code]

which is pretty much adding a command to what I have lol.


All times are UTC. The time now is 23:05.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.