nope vecsmall only stores it as integers by the looks of it but I can't get all of them to work.

[QUOTE=CRGreathouse;240430]You can't store anything but ASCII characters in Pari strings (because, underneath, they're just C strings which have the same limitation).[/QUOTE]

[url]http://evanjones.ca/unicode-in-c.html[/url] has how they do it in C wonder if we can implement the same in Pari lol.

[QUOTE=science_man_88;240432][url]http://evanjones.ca/unicode-in-c.html[/url] has how they do it in C wonder if we can implement the same in Pari lol.[/QUOTE]

I alluded to this in post #2035. It's quite possible, though it would require perhaps a thousand changes in the Pari source code (mostly small).

[QUOTE=science_man_88;240431]nope vecsmall only stores it as integers by the looks of it but I can't get all of them to work.[/QUOTE]

Reread my post.

[QUOTE=CRGreathouse;240426]
[quote]
- The expression "sup[n]=sup[n]+32" won't work with "sup[n]+=32" !
[/quote]
It should...

[CODE]> v=vector(3,n,n);v[2]+=10;v
%1 = [1, 12, 3]
> v=vectorsmall(3,n,n);v[2]+=10;v
%2 = Vecsmall([1, 12, 3])[/CODE]

[/QUOTE]

I used that one directly in the Strchr-function:

[code]
s=concat(s,Strchr(sup[n]-=32))
[/code]

and that will produce the error:

[code]
*** _-=_: bug in PARI/GP (Segmentation Fault), please report
[/code]

Perhaps something really to report!?

But the code
[code]
s=concat(s,Strchr(sup[n]=sup[n]-32))
[/code]
works here fine.

[QUOTE=kar_bon;240520]I used that one directly in the Strchr-function:

[code]
s=concat(s,Strchr(sup[n]-=32))
[/code]

and that will produce the error:

[code]
*** _-=_: bug in PARI/GP (Segmentation Fault), please report
[/code]

Perhaps something really to report!?

But the code
[code]
s=concat(s,Strchr(sup[n]=sup[n]-32))
[/code]
works here fine.[/QUOTE]

Interesting. I may report this.

Regardless, you don't need to do that here; you can just do
[CODE]StrToUp(str)={
my(sup=Vecsmall(str),s="");
for(n=1,#sup,
if(sup[n]>=97 && sup[n]<=122,
s=concat(s,Strchr(sup[n]-32))
,
s=concat(s,Strchr(sup[n]))
)
);
s
};[/CODE]
since you don't use sup later.

It's been somewhat difficult to find a factor for sigma(4477204945765526980349352291678647292170001931988428673↑2), or 20045364126387295389488107373344241624560499868937574184788150529738565510740963223091652652278528299596969603.

So far, I've found no factors ≤ 30 digits.

anyone want to see the encoding I made to show for 40 questions I don't think it will hurt it too much lol.

[CODE]coding(string) = a=Vecsmall(string);for(y=1,#a,vecsort(a);for(x=1,#a,a[x]=a[x]+x));for(z=1,#a,if(z%2==0,print1(Strchr(a[z]%256)),print1(a[z])));print1("\n"string)[/CODE]

I got this recently lol :

[CODE](21:17)>coding(string) = a=Vecsmall(string);vecsort(a);for(x=1,#a,a[x]+=x)
%32 = (string)->a=Vecsmall(string);vecsort(a);for(x=1,#a,a[x]+=x)
(21:17)>coding("string")
[COLOR="Red"]*** _+=_: bug in PARI/GP (Segmentation Fault), please report[/COLOR][/CODE]

[QUOTE=3.14159;240629]It's been somewhat difficult to find a factor for sigma(4477204945765526980349352291678647292170001931988428673↑2), or 20045364126387295389488107373344241624560499868937574184788150529738565510740963223091652652278528299596969603.

So far, I've found no factors ≤ 30 digits.[/QUOTE]

Nothing to 45 digits (though I'm not quite finished with the ECM). It's almost surely a semiprime, and a hard one at that. You might have to use NFS if you want to crack it.

[CODE](14:14)>forprime(x=1,6000,c=vector(1000,n,0);c[1]=x;for(y=2,#c,c[y]=2*c[y-1]+1;for(z=1,#mersenne,if(c[y]!=mersenne[z],,print(y","z","x);break(2)))))
2,3,2
2,4,3
3,8,7
7,15,19
3,12,31
5,14,37
2,11,53
5,15,79
3,14,151
(14:16)>
[/CODE]

okay so there are exceptions to my thoughts of for each exception all the rest never hit mersenne exponents in this way however 4 in the primes under 6000 isn't that bad is it ? I might test under 10000 next.

funny how the start of the 2^x-1 chain that seems to work for this starts with 2047 and I posted post 2047.

[QUOTE=science_man_88;241130][CODE](14:14)>forprime(x=1,6000,c=vector(1000,n,0);c[1]=x;for(y=2,#c,c[y]=2*c[y-1]+1;for(z=1,#mersenne,if(c[y]!=mersenne[z],,print(y","z","x);break(2)))))
2,3,2
2,4,3
3,8,7
7,15,19
3,12,31
5,14,37
2,11,53
5,15,79
3,14,151
(14:16)>
[/CODE]

okay so there are exceptions to my thoughts of for each exception all the rest never hit mersenne exponents in this way however 4 in the primes under 6000 isn't that bad is it ? I might test under 10000 next.

funny how the start of the 2^x-1 chain that seems to work for this starts with 2047 and I posted post 2047.[/QUOTE]


made my code faster I think lol and tested it with all primes under 500000.

[CODE]forprime(x=1,2000,c=vector(1000,n,0);c[1]=x;for(y=2,#c,c[y]=2*c[y-1]+1;if(c[y]<mersenne[#mersenne],for(z=1,#mersenne,if(c[y]==mersenne[z],print(y","z","x);break(2))),break())))[/CODE] I have tested it under 500000 and got about 9 seconds I think. mind you a lot go into the untested range.

[QUOTE=CRGreathouse;240704]Nothing to 45 digits (though I'm not quite finished with the ECM). It's almost surely a semiprime, and a hard one at that. You might have to use NFS if you want to crack it.[/QUOTE]

Or three to five days of SIQS..

[QUOTE=3.14159;241167]Or three to five days of SIQS.[/QUOTE]

SIQS would work but would be quite a bit slower than NFS.

alt(x)=if(x<128,Strchr(x),a=Vec(read("E:\\alt.txt"));print(a[x-127]))

I tried this for doing proper alt characters after 127 as my PARI didn't on it's own.

the problem is my PARI won't read anything not in a gp binary anymore and even when i renamed it .gp it told me it wasn't a gp binary but my codes file worked.

[QUOTE=science_man_88;241280][CODE]alt(x)=if(x<128,Strchr(x),a=Vec(read("E:\\alt.txt"));print(a[x-127]))
[/CODE]
I tried this for doing proper alt characters after 127 as my PARI didn't on it's own.

the problem is my PARI won't read anything not in a gp binary anymore and even when i renamed it .gp it told me it wasn't a gp binary but my codes file worked.[/QUOTE]

forgot the code tags.

From the other [url=http://www.mersenneforum.org/showthread.php?t=13181]thread[/url]:

Evaluate the Sum of p^p, p prime, up to m-th prime:
[code]
addhelp(SumPPowerP, "SumPPowerP(m): Returns the sum of p^p, p prime up to m-th prime p.");
SumPPowerP(m)={
my(SumPPowerP=0);
for(n=1,m,SumPPowerP+=prime(n)^prime(n));
SumPPowerP
};
[/code]

and

find the value n for which SumPPowerP(n) mod 10^m == 0
[code]
addhelp(FindSumPPPMod10, "FindSumPPPMod10(m): Return n of (SumPPowerP(n) mod 10^m == 0).");
FindSumPPPMod10(m)={
n=1;
until(SumPPowerP(n)%10^m==0,
n++
);
n
};
[/code]

so
prime(FindSumPPPMod10(1)) finds 11,
prime(FindSumPPPMod10(2)) finds 751,
prime(FindSumPPPMod10(3)) finds 1129

With an own procedure (with loop and printing) you got the values with one call.

Note: The value for n=4 will take some time!

[QUOTE=kar_bon;241663]From the other [url=http://www.mersenneforum.org/showthread.php?t=13181]thread[/url]:

<snip>

Note: The value for n=4 will take some time![/QUOTE]

Oh, don't do that. Seriously. That's very inefficient. Sum(n+1) = Sum(n) + prime(n+1)^(prime(n+1). So computing the sum every time from scratch is really really bad (N^2 vs N).

[CODE](18:15)>v=vector(1000,n,(prime(n)^prime(n))%10)
%189 = [4, 7, 5, 3, 1, 3, 7, 9, 7, 9, 1, 7, 1, 7, 3, 3, 9, 1, 3, 1, 3, 9, 7, 9, 7, 1, 7, 3, 9
(18:16)>a=0;for(m=1,100,for(x=1,#v,a=a+v[x];b=10^m;if(a%b==0,print(prime(x));break()));a=0)
11
661
4397[/CODE]

future [URL="http://www.youtube.com/watch?v=17DPJHNVx2Q&feature=player_embedded"]thec[/URL]

[QUOTE=axn;241667]Oh, don't do that. Seriously. That's very inefficient. Sum(n+1) = Sum(n) + prime(n+1)^(prime(n+1). So computing the sum every time from scratch is really really bad (N^2 vs N).[/QUOTE]

Not really... the values increase (heuristically) super-exponentially (by which I mean not in E), so the expected penalty is ~11% rather than (N-1) * 100%.

[QUOTE=kar_bon;241663]Evaluate the Sum of p^p, p prime, up to m-th prime:
[code]
addhelp(SumPPowerP, "SumPPowerP(m): Returns the sum of p^p, p prime up to m-th prime p.");
SumPPowerP(m)={
my(SumPPowerP=0);
for(n=1,m,SumPPowerP+=prime(n)^prime(n));
SumPPowerP
};
[/code][/QUOTE]

I'll give you the same advice I gave sm on the other thread: whenever you find yourself writing

for(n=1, _, ... prime(n) ...)

you should probably be writing

forprime(p=2, prime(_), ... p ...)

instead. It's significantly faster. (Actually, in this case it isn't because the powerings take a long time, but in general this can make a big difference.)

Also, since the numbers get so large so quickly, I recommend working mod the appropriate power of ten rather than storing the whole number and reducing. So
[code]addhelp(SumPPowerPModM, "SumPPowerP(lim,m): Returns the sum of p^p, p prime and <= lim, mod m.");
SumPPowerPModM(lim, m)={
my(SumPPowerP=Mod(0,m));
forprime(p=2,lim,SumPPowerP+=Mod(p,m)^p);
lift(SumPPowerP)
};[/code]

Ok, thanks for the advice.

The function SumPPowerP(m) was meant to be as it is.

I would create a new function like you did then for the MOD-values.

And further more it's possible, to use file reading and writing for higher m-values to calculate the next power from a saved file.

I've done this before your post, so a combination of file-read/write and MOD should be best.

Here's what I've done for my original functions:

- Create a folder named "PPowerP" (in PARI-installation folder) to hold those save-files

To create a new file for PPowerP(m) use:
[code]
addhelp(WritePPowerPFile, "WritePPowerPFile(m): Write PPowerP(m) into a file.");
WritePPowerPFile(m)={
write("PPowerP/PPowerP"m".txt",SumPPowerP(m));
};
[/code]

For example "WritePPowerPFile(1000)" will create a file called "PPowerP1000.txt" with the sum (whole number = 30874 digits).

To create from this file the next say for m=2000 use
[code]
addhelp(WritePPowerPFile, "WritePPowerPFile(m): Write PPowerP(n) to file with template from PPowerP(m)-file.");
WritePPowerPFromFile(m,n)={
my(PP=read("PPowerP/PPowerP"m".txt"));
for(x=m+1,n,PP+=prime(x)^prime(x));
write("PPowerP/PPowerP"n".txt",PP);
};
[/code]

by calling "WritePPowerPFromFile(1000,2000)" and it will create the file "PPowerP2000.txt".

I've not yet inserted the 'forprime' usage and the MOD but this should not such big problem here.
Perhaps a trap should be used to exclude errors when a file is not found.

[QUOTE=CRGreathouse;241692]Not really... the values increase (heuristically) super-exponentially (by which I mean not in E), so the expected penalty is ~11% rather than (N-1) * 100%.[/QUOTE]

Huh?! His posted code had a nested loop, while the "proper" way to do is a straight loop. N^2 vs N.

[QUOTE=axn;241719]Huh?! His posted code had a nested loop, while the "proper" way to do is a straight loop. N^2 vs N.[/QUOTE]

Oh, I see what you're talking about. Yes, that's inefficient, though still not N vs. N[SUP]2[/SUP] since the work grows greatly with the increasing numbers; I get more like N vs. 0.47N[SUP]1.5[/SUP] taking that into account.

In short: you were (essentially) right, I was wrong. :blush:

[CODE][COLOR="Lime"]x=vector(100,n,if(isprime(n),n));x=vecsort(x);for(y=1,#x-1,if(y<=#x && x[y]==0,x=vector(#x-1,n,x[n+1]);x=concat(x,0)))[/COLOR];for(b=1,#x,if(b=<#x && x[b]==0,x=vector((b-1),n,x[n]);break()))[/CODE]

I don't see what's wrong with this code it works doing the green section alone, but won't work for the last section tells me I have unused characters. I don't get why it doesn't work. The code finds the primes then vecsorts all the 0's to the front then gets it to push the primes forward . But now I can't cleave off the 0's like i used to.

I get a different error message:
[code] *** syntax error, unexpected '<': for(b=1,#x,if(b=<#x&&x[b]==0,x=vecto[/code]
The black section of your code includes a =< which should be <=. Once I change that it works for me.

[QUOTE=CRGreathouse;241832]I get a different error message:
[code] *** syntax error, unexpected '<': for(b=1,#x,if(b=<#x&&x[b]==0,x=vecto[/code]
The black section of your code includes a =< which should be <=. Once I change that it works for me.[/QUOTE]

same here now lol.

but now it seems to only works for random lengths that I have no patern to lol 1068 works as a limit 1000 doesn't 100 worked but 108 or 109 wouldn't.

I named my code to show you:

[CODE](17:28)>for(x=1,100,print(primesUPTO(x)","x))
[],1
[2],2
[2, 3],3
[],4
[],5
[2, 3, 5],6
[2, 3, 5, 7],7
[],8
[2, 3, 5, 7],9
[],10
[],11
[2, 3, 5, 7, 11],12
[2, 3, 5, 7, 11, 13],13
[],14
[2, 3, 5, 7, 11, 13],15
[],16
[],17
[2, 3, 5, 7, 11, 13, 17],18
[2, 3, 5, 7, 11, 13, 17, 19],19
[],20
[2, 3, 5, 7, 11, 13, 17, 19],21
[],22
[],23
[2, 3, 5, 7, 11, 13, 17, 19, 23],24
[],25
[2, 3, 5, 7, 11, 13, 17, 19, 23],26
[],27
[2, 3, 5, 7, 11, 13, 17, 19, 23],28
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29],29
[],30
[],31
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31],32
[],33
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31],34
[],35
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31],36
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37],37
[],38
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37],39
[],40
[],41
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41],42
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43],43
[],44
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43],45
[],46
[],47
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47],48
[],49
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47],50
[],51
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47],52
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53],53
[],54
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53],55
[],56
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53],57
[],58
[],59
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59],60
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61],61
[],62
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61],63
[],64
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61],65
[],66
[],67
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67],68
[],69
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67],70
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71],71
[],72
[],73
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73],74
[],75
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73],76
[],77
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73],78
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79],79
[],80
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79],81
[],82
[],83
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83],84
[],85
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83],86
[],87
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83],88
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89],89
[],90
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89],91
[],92
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89],93
[],94
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89],95
[],96
[],97
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97],98
[],99
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97],100[/CODE]

by the looks of it it fails for numbers of 6x+4 6x+5 ,6x+8(+2)

[CODE]sieve_of_eratosthenes(x) = a=0;c=[2,3];for(y=4,x,for(z=1,#c,if(y%c[z]==0,break(),a=a+1));if(a==#c,c=concat(c,y));a=0);c[/CODE]

This is a lot slower but I can't find a failure of it.

[QUOTE=science_man_88;241855][CODE]sieve_of_eratosthenes(x) = a=0;c=[2,3];for(y=4,x,for(z=1,#c,if(y%c[z]==0,break(),a=a+1));if(a==#c,c=concat(c,y));a=0);c[/CODE]

This is a lot slower but I can't find a failure of it.[/QUOTE]

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.

Like this:

[code]
SoE(n)={
my(v=vector(n,m,1));
v[1]=0;
forstep(i=4,n,2,v[i]=0);
for(i=3,n\2,
if(v[i]==1,
forstep(j=2*i,n,i,v[j]=0)));
v
};

PrintSoE(n)={
my(v=SoE(n),j=0);
for(i=1,n,if(v[i],print1(i" ");j++));
print;
print(j" primes from 1 to "n);
};
[/code]

[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.

[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.

not sure how to turn it into Atkin or Sundaram

[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!

[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.

[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.

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]

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

[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.

[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.

[CODE]c=vector(10, n, n)
vecextract(c, concat("1..", #c-1))
vecextract(c, concat("2..", #c))[/CODE]

[CODE]c=concat(c,c[1]);c=vector(#c-1,n,c[n+1])[/CODE]

I think a couple of this add-on to the basics would be nice as 2 of them could possibly be used to crack a basic shift cipher.

[QUOTE=axn;241984][CODE]c=vector(10, n, n)
vecextract(c, concat("1..", #c-1))
vecextract(c, concat("2..", #c))[/CODE][/QUOTE]

didn't know I could do that syntax.

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

Meh. I don't quite know the code to make a custom vector that returns something other than 1 and 0, other than insert code into the vector command.

[QUOTE=3.14159;242066]Meh. I don't quite know the code to make a custom vector that returns something other than 1 and 0, other than insert code into the vector command.[/QUOTE]

what exactly do you want to return ?

[CODE]x=vector(100,n,n) ; makes x have 100 indexes with each indexes value equals to the index.[/CODE]
[CODE]x=vector(100,n,if(isprime(n,n))) ; makes a vector of 100 indexes with the index value equal to the index if the index value is prime. otherwise the index value is 0.[/CODE]

[QUOTE=3.14159;242066]Meh. I don't quite know the code to make a custom vector that returns something other than 1 and 0, other than insert code into the vector command.[/QUOTE]

Well, you need to make a vector, either from scratch with the vector command or by converting something to a vector with Vec. Once you have that you can put in values at appropriate indices as needed, then return it.

Certainly vectors aren't appropriate for all tasks, but it's good to use them when they are.

----

Here's a trick* I just learned. When you want to build a vector element-by-element, a better way than
v = concat(v, new_element)
is to make a list, store elements with listput, then convert (if needed) to a vector with Vec. So
[code]natural(N)={
my(l=List());
for(n=1,N,listput(l,n));
Vec(l)
};[/code]
is faster than
[code]natural1(N)={
my(v=[]);
for(n=1,N,v=concat(v,n));
v
};[/code]

Actually, with the changes made to concat() a few versions ago, the difference isn't great unless you aren't doing too much calculation, and it's better to make the whole vector at once when you can (vector(N,i,i) will be faster than either in this case). But I thought I'd share in case anyone was interested.

* Not "something clever" but "the way I always should have been doing it".

[QUOTE=CRGreathouse;242108]So
[code]natural(N)={
my(l=List());
for(n=1,N,listput(l,n));
Vec(l)
};[/code]
is faster than
[code]natural1(N)={
my(l=List());
for(n=1,N,listput(l,n));
Vec(l)
};[/code]
[/QUOTE]

I seriously doubt if the first one is faster than the second :wink:

Hah! Fixed it.

[QUOTE=axn;242113]I seriously doubt if the first one is faster than the second :wink:[/QUOTE]

[CODE](09:42)>natural(100000)
%3 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
1851, 1852, 1853, 1854, 1855, 1856, 1857, 185
7, 3518, 3519, 3520, 3521, 3522, 3523, 3524,
5184, 5185, 5186, 5187, 5188, 5189, 5190, 519
0, 6851, 6852, 6853, 6854, 6855, 6856, 6857,
8517, 8518, 8519, 8520, 8521, 8522, 8523, 852
157, 10158, 10159, 10160, 10161, 10162, 10163
, 11586, 11587, 11588, 11589, 11590, 11591, 1
3014, 13015, 13016, 13017, 13018, 13019, 1302
2, 14443, 14444, 14445, 14446, 14447, 14448,
15871, 15872, 15873, 15874, 15875, 15876, 158
99, 17300, 17301, 17302, 17303, 17304, 17305,
18728, 18729, 18730, 18731, 18732, 18733, 18
156, 20157, 20158, 20159, 20160, 20161, 20162
, 21585, 21586, 21587, 21588, 21589, 21590, 2
3013, 23014, 23015, 23016, 23017, 23018, 2301
1, 24442, 24443, 24444, 24445, 24446, 24447,
25870, 25871, 25872, 25873, 25874, 25875, 258
98, 27299, 27300, 27301, 27302, 27303, 27304,
28727, 28728, 28729, 28730, 28731, 28732, 28
155, 30156, 30157, 30158, 30159, 30160, 30161
, 31584, 31585, 31586, 31587, 31588, 31589, 3
3012, 33013, 33014, 33015, 33016, 33017, 3301
0, 34441, 34442, 34443, 34444, 34445, 34446,
35869, 35870, 35871, 35872, 35873, 35874, 358
97, 37298, 37299, 37300, 37301, 37302, 37303,
38726, 38727, 38728, 38729, 38730, 38731, 38
154, 40155, 40156, 40157, 40158, 40159, 40160
, 41583, 41584, 41585, 41586, 41587, 41588, 4
3011, 43012, 43013, 43014, 43015, 43016, 4301
9, 44440, 44441, 44442, 44443, 44444, 44445,
45868, 45869, 45870, 45871, 45872, 45873, 458
96, 47297, 47298, 47299, 47300, 47301, 47302,
48725, 48726, 48727, 48728, 48729, 48730, 48
153, 50154, 50155, 50156, 50157, 50158, 50159
, 51582, 51583, 51584, 51585, 51586, 51587, 5
3010, 53011, 53012, 53013, 53014, 53015, 5301
8, 54439, 54440, 54441, 54442, 54443, 54444,
55867, 55868, 55869, 55870, 55871, 55872, 558
95, 57296, 57297, 57298, 57299, 57300, 57301,
58724, 58725, 58726, 58727, 58728, 58729, 58
152, 60153, 60154, 60155, 60156, 60157, 60158
, 61581, 61582, 61583, 61584, 61585, 61586, 6
3009, 63010, 63011, 63012, 63013, 63014, 6301
7, 64438, 64439, 64440, 64441, 64442, 64443,
65866, 65867, 65868, 65869, 65870, 65871, 658
94, 67295, 67296, 67297, 67298, 67299, 67300,
68723, 68724, 68725, 68726, 68727, 68728, 68
151, 70152, 70153, 70154, 70155, 70156, 70157
, 71580, 71581, 71582, 71583, 71584, 71585, 7
3008, 73009, 73010, 73011, 73012, 73013, 7301
6, 74437, 74438, 74439, 74440, 74441, 74442,
75865, 75866, 75867, 75868, 75869, 75870, 758
93, 77294, 77295, 77296, 77297, 77298, 77299,
78722, 78723, 78724, 78725, 78726, 78727, 78
150, 80151, 80152, 80153, 80154, 80155, 80156
, 81579, 81580, 81581, 81582, 81583, 81584, 8
3007, 83008, 83009, 83010, 83011, 83012, 8301
5, 84436, 84437, 84438, 84439, 84440, 84441,
85864, 85865, 85866, 85867, 85868, 85869, 858
92, 87293, 87294, 87295, 87296, 87297, 87298,
88721, 88722, 88723, 88724, 88725, 88726, 88
149, 90150, 90151, 90152, 90153, 90154, 90155
, 91578, 91579, 91580, 91581, 91582, 91583, 9
3006, 93007, 93008, 93009, 93010, 93011, 9301
4, 94435, 94436, 94437, 94438, 94439, 94440,
95863, 95864, 95865, 95866, 95867, 95868, 958
91, 97292, 97293, 97294, 97295, 97296, 97297,
98720, 98721, 98722, 98723, 98724, 98725, 98
(09:42)>##
*** last result computed in 125 ms.
(09:42)>natural1(N)={
my(v=[]);
for(n=1,N,v=concat(v,n));
v
};
(09:42)>natural1(100000)
*** user interrupt after 3mn, 39,875 ms.[/CODE]

if 65536 was the limit for Vec ( which I thought it would be ) then the faster script also gets around that lol.

those test are with minimal memory with my maximum I was able to do natural(10000000) in under 17 seconds according to the timing it gave. I wonder how much isprime would slow something like this down lol. looks like me adding a isprime check and trying up to 10000000 actually sped it up, listed all primes under 10 million in under 10 seconds according to the timer.That's without turning it into a forstep loop to only check 6n+1 of 6n-1.

[QUOTE=science_man_88;242171]those test are with minimal memory with my maximum I was able to do natural(10000000) in under 17 seconds according to the timing it gave.[/QUOTE]

Yes -- natural1 would take a long, long time to do that. It took 4:59.8 to do 10^5; based on curve-fitting with smaller numbers I estimate it would take about a month for it to do 10^7.

[QUOTE=CRGreathouse;242184]Yes -- natural1 would take a long, long time to do that. It took 4:59.8 to do 10^5; based on curve-fitting with smaller numbers I estimate it would take about a month for it to do 10^7.[/QUOTE]

yeah with my isprime check and forstep loop put/editted in this seems even faster as it reduces the work. I was just trying:

[CODE]natural2(10000000)[/CODE] which is natural with those replacements etc. it comes back with:

[CODE]*** last result computed in 7,344 ms.[/CODE] and even when i tried 100 million I was able to get every prime in under 1 min 20 I think.

Let me see the code.

[QUOTE=CRGreathouse;242193]Let me see the code.[/QUOTE]

[CODE]natural2(N)=my(l=List());forstep(n=5,N,[2,4],if(isprime(n),listput(l,n)));Vec(l);[/CODE]

did I mention I'm running cityville on Facebook as well lol.

[CODE](12:43)>##
*** last result computed in 1mn, 18,328 ms.
(12:43)>natural2(100000000)[/CODE]

pretty quick.

[QUOTE=science_man_88;242194][CODE]natural2(N)=my(l=List());forstep(n=5,N,[2,4],if(isprime(n),listput(l,n)));Vec(l);[/CODE]

did I mention I'm running cityville on Facebook as well lol.[/QUOTE]

So it's based on the fast natural, not the slow natural1. No surprise it's faster than natural1.

[QUOTE=CRGreathouse;242235]So it's based on the fast natural, not the slow natural1. No surprise it's faster than natural1.[/QUOTE]

last I checked it was faster than the fast natural

[QUOTE=science_man_88;242237]last I checked it was faster than the fast natural[/QUOTE]

Sure -- it doesn't add as many numbers.

[QUOTE=CRGreathouse;242239]Sure -- it doesn't add as many numbers.[/QUOTE]

outside of the fact that the default list seems faster if I could get it fast enough could it ever replace the current one ?

[QUOTE=science_man_88;242244]outside of the fact that the default list seems faster if I could get it fast enough could it ever replace the current one ?[/QUOTE]

I'm not sure what you're asking.

If the question is, "could some version of natural2 replace natural?", I'm not sure how to answer -- natural2 gives different output from natural, and natural itself is just a less-efficient version of vector(N,i,i) which doesn't have many uses anyway.

If the question is, "could some version of List replace List?", it's possible -- but that would have to be in C/PARI not gp.

[QUOTE=CRGreathouse;242248]I'm not sure what you're asking.

If the question is, "could some version of natural2 replace natural?", I'm not sure how to answer -- natural2 gives different output from natural, and natural itself is just a less-efficient version of vector(N,i,i) which doesn't have many uses anyway.

If the question is, "could some version of List replace List?", it's possible -- but that would have to be in C/PARI not gp.[/QUOTE]

I'm talking [CODE]default(primelimit,x)[/CODE] seems faster than this script .

[QUOTE=science_man_88;242291]I'm talking [CODE]default(primelimit,x)[/CODE] seems faster than this script .[/QUOTE]

Wow, that was random.

Sure. So is the program
[code]1[/code]

[CODE]a=0;for(i=1,#mersenne,for(x=a+1,100,if(prime(x)==mersenne[i],a=x;break(),print1(x",");a=x;break())))[/CODE]

Basically checking for [url]http://oeis.org/A135980[/url] , but it's not working as intended. I'm using a so x doesn't repeat a value, this seems to work however it reports 6 as the second term not 9 so I've made a major error. doh I should have reversed the order of the loops I think. I see why it prints six because it checks 11==13 and that's false so it prints x then puts 5 in a so it then increments without finding the next match and hence it returns 6 on finding 13!=17.

[QUOTE=science_man_88;243003][CODE]a=0;for(i=1,#mersenne,for(x=a+1,100,if(prime(x)==mersenne[i],a=x;break(),print1(x",");a=x;break())))[/CODE]

Basically checking for [url]http://oeis.org/A135980[/url] , but it's not working as intended. I'm using a so x doesn't repeat a value, this seems to work however it reports 6 as the second term not 9 so I've made a major error. doh I should have reversed the order of the loops I think. I see why it prints six because it checks 11==13 and that's false so it prints x then puts 5 in a so it then increments without finding the next match and hence it returns 6 on finding 13!=17.[/QUOTE]

note putting a i=i-1 in the if false part worked.

[CODE]aliquot(x,y)=b=0;v=vector(y);v[1]=x;for(i=2,#v,for(a=1,i-1,if(sigma(v[i-1])-v[i-1]==0 || sigma(v[i-1])-v[i-1]==v[a],break(2),b=b+1;if(b!=a,b=0;break(2),v[i]=sigma(v[i-1])-v[i-1];b=0;break()))))[/CODE]

this seem to work for aliquot checking repeat lol.Now I just have to eliminate the y variable and cleave it to the final length if not less than infinite I think. too bad i can't get it to start at any index already calculated.

[QUOTE=science_man_88;244430][CODE]aliquot(x,y)=b=0;v=vector(y);v[1]=x;for(i=2,#v,for(a=1,i-1,if(sigma(v[i-1])-v[i-1]==0 || sigma(v[i-1])-v[i-1]==v[a],break(2),b=b+1;if(b!=a,b=0;break(2),v[i]=sigma(v[i-1])-v[i-1];b=0;break()))))[/CODE]

this seem to work for aliquot checking repeat lol.Now I just have to eliminate the y variable and cleave it to the final length if not less than infinite I think. too bad i can't get it to start at any index already calculated.[/QUOTE]

crap a failure.

[CODE]for(x=1,100, print(aliquot(x,5)))[/CODE]

gave 95,25,6,6,6, but it should have terminated after finding the second 6.

[CODE]aliquot(x,y)=b=0;v=vector(y);v[1]=x;for(i=2,#v,for(a=1,i-1,if(sigma(v[i-1])-v[i-1]==0 || sigma(v[i-1])-v[i-1]==v[a],break(2),b=b+1;if(b==i-1,v[i]=sigma(v[i-1])-v[i-1];b=0;break()))))[/CODE] this works better for now.

[CODE]aliquot(x,y)=b=0;v=vector(y);v[1]=x;for(i=2,#v,for(a=1,i-1,if(sigma(v[i-1])-v[i-1]==0 || sigma(v[i-1])-v[i-1]==v[a],v=vector(i-1,n,v[n]);break(2),b=b+1;if(b==i-1,v[i]=sigma(v[i-1])-v[i-1];b=0;break()))))[/CODE] this not only works but trims off excess.

[QUOTE=CRGreathouse;229472]kar_bon, I understand what an aliquot sequence is, but that doesn't tell me what sm88 wants.

I take it from his last post that he wants a vector with last term 1 and first term the input, where v[i+1] = sigma(v[i]) - v[i]. But of course it's not even proven that such a vector exists for all n...![/QUOTE]

I see one possible exception but I don't know how to calculate when it would happen, a sequence with all abundant numbers as it would forever increase.

[QUOTE=science_man_88;244558]I see one possible exception but I don't know how to calculate when it would happen, a sequence with all abundant numbers as it would forever increase.[/QUOTE]

good news I've been testing x upto 3 million until the 4th spot in the sequence each and no x so far has made it past 3 terms that are abundant. the problem is if it decreases to another abundant number spiking it back up. Doh calculation error. redone and still accurate.

done up to 7 million now without a trace of a triply ? abundant number.

I couldn't get the addhelp to add help I've had it working before but I just tried the autosave that lavalamp got me doing and I can't get what I want to work.

[CODE](15:45)>for(x=1,400,if(((2^x)/(log(2^x)-1))-((2^x-1)/(log(2^x-1)-1))>=1,print(x)))
108
116[/CODE] this is my attempt at something new if only I could get it to align with A000043.

[QUOTE=science_man_88;250325]I couldn't get the addhelp to add help I've had it working before but I just tried the autosave that lavalamp got me doing and I can't get what I want to work.[/QUOTE]

Can you give an example? For example, what happens if you type
[code]foo()="bar";
addhelp(foo, "foo(): Bar.")
?foo[/code]
?

[QUOTE=CRGreathouse;250363]Can you give an example? For example, what happens if you type
[code]foo()="bar";
addhelp(foo, "foo(): Bar.")
?foo[/code]
?[/QUOTE]

if works but I can't get it to work if i save it in file like that I bet. ? how'd it work I couldn't get mine to work grrr. I bet I messed up. also can you save vectors in file and call them back ? I can't get it or I'd save the MPE's it might save me a while lol. got it thanks again you are very useful.

[QUOTE=science_man_88;250367]also can you save vectors in file and call them back ?[/QUOTE]

Yep. I have a few common ones in my auto-loaded script.

To sm88: are you working on mersennes, perfects, abundants,
deficients, aliquots, riesels, sierpinskis, and brilliants too?

They are intimately connected.

Why not give us your versions of their definitions, just for fun.

pari mutual question
 

How does one run pari?

[QUOTE=science_man_88;243003][CODE]a=0;for(i=1,#mersenne,for(x=a+1,100,if(prime(x)==mersenne[i],a=x;break(),print1(x",");a=x;break())))[/CODE]Basically checking for [URL]http://oeis.org/A135980[/URL] , but it's not working as intended. I'm using a so x doesn't repeat a value, this seems to work however it reports 6 as the second term not 9 so I've made a major error. doh I should have reversed the order of the loops I think. I see why it prints six because it checks 11==13 and that's false so it prints x then puts 5 in a so it then increments without finding the next match and hence it returns 6 on finding 13!=17.[/QUOTE]

Not working as intended because why?

Re-think your whole programming methodology.

[QUOTE=davar55;250409]How does one run pari?[/QUOTE]

Download the binary at
[url]http://pari.math.u-bordeaux.fr/[/url]
(or use your package manager, if on *nix) and run.

CRG do you think I'm able enough in PARI to attempt the BWT algorithm of encoding ? it uses RLE and I've coded that. the main hard part is remembering all the vector alteration codes. I have to rotate through all combinations ( putting them in a vector may help) then I'd have to vecsort, alphabetically then take the last character of each element in the vector and then perform RLE. I also have to have a state so I can tell if compression or decompression is wanted I could use a variable there. I think the decompression si the hard part and I don't know if I ever tried decompression for RLE.

hi is there a way in PARI to tell a function to only accept vectors and rename them to use them in the script I'm working on something of a aredisjoint(A,B) script. but I want the arguments to be treated as vectors:

[CODE]aredisjoint(A,B) = a=0;for(i=1,#A,for(j=1,#B,if(A[i]==B[j],return(0),a=a+1;if(a==#A,return(1)))))[/CODE]
I gave it a test and it failed it it always seems to return 1 even if the values are the same. if I get this to work I may start on a ispd (is pairwise disjoint) function. never mind I found my error it should be if(a==#A*#B

final script with one other change and I get:

[CODE]aredisjoint(A=[],B=[]) = if(#A!=0 && #B!=0,a=0;for(i=1,#A,for(j=1,#B,if(A[i]==B[j],return(0),a=a+1;if(a==#A*#B,return(1))))),return(0))[/CODE]

[CODE]ispd(C) = b=0;for(i=1,#C,for(j=1,#C,if(i!=j && aredisjoint(C[i],C[j]),b=b+1)));if(b==#C*#C-#C,return(1),return(0))[/CODE] I checked this on a vector with 2 other vectors inside it.

I was also thinking of making a Cartesian product making script but my main problem is I don't know how to allow for infinite arguments in the function. I guess I could make one that worked for 2 at a time, and puts the parentheses around it with a complete argument. what am i saying I can make it so the person has to give a collection of sets doh.

[CODE](12:17)>p=[1,2];c=[1,2];d=[];for(i=1,#p,for(j=1,#c,x=[p[i],c[j]];d=concat(d,[x])))
(12:17)>d
%201 = [[1, 1], [1, 2], [2, 1], [2, 2]]
(12:17)>p=[1,2];e=[];for(i=1,#p,for(j=1,#d,x=[d[j],p[i]];e=concat(e,[x])))
(12:17)>e
%202 = [[[1, 1], 1], [[1, 2], 1], [[2, 1], 1], [[2, 2], 1], [[1, 1], 2], [[1, 2], 2], [[2, 1], 2], [[2, 2], 2]][/CODE]

found the basics but the problem I have is how to make this more like tuples. I think I have a way but it would take forever with the method I'm thinking of.

[CODE](12:34)>Cp(p,c)= d=[];for(i=1,#p,for(j=1,#c,x=[p[i],c[j]];d=concat(d,[x])))
%205 = (p,c)->d=[];for(i=1,#p,for(j=1,#c,x=[p[i],c[j]];d=concat(d,[x])))
(12:39)>Cp([1,3],[1,2])
(12:39)>d
%206 = [[1, 1], [1, 2], [3, 1], [3, 2]]
(12:39)>CPC(C)= A=Cp(C[1],C[2]);B=[];for(i=3,#C,if(i%2==1,B=Cp(A,C[i]),A=Cp(B,C[i])));if(#C%2==1,return(B),return(A))
%207 = (C)->A=Cp(C[1],C[2]);B=[];for(i=3,#C,if(i%2==1,B=Cp(A,C[i]),A=Cp(B,C[i])));if(#C%2==1,return(B),return(A))
(12:43)>C=[v,c]
%208 = [[1, 2, 3, 4, 5, 6, 7, 8], [1, 2]]
(12:44)>CPC(C)
%209 = 0[/CODE]

not sure why the last result and I know I forgot the empty set.

[QUOTE=science_man_88;251276][CODE](12:34)>Cp(p,c)= d=[];for(i=1,#p,for(j=1,#c,x=[p[i],c[j]];d=concat(d,[x])))
%205 = (p,c)->d=[];for(i=1,#p,for(j=1,#c,x=[p[i],c[j]];d=concat(d,[x])))
(12:39)>Cp([1,3],[1,2])
(12:39)>d
%206 = [[1, 1], [1, 2], [3, 1], [3, 2]]
(12:39)>CPC(C)= A=Cp(C[1],C[2]);B=[];for(i=3,#C,if(i%2==1,B=Cp(A,C[i]),A=Cp(B,C[i])));if(#C%2==1,return(B),return(A))
%207 = (C)->A=Cp(C[1],C[2]);B=[];for(i=3,#C,if(i%2==1,B=Cp(A,C[i]),A=Cp(B,C[i])));if(#C%2==1,return(B),return(A))
(12:43)>C=[v,c]
%208 = [[1, 2, 3, 4, 5, 6, 7, 8], [1, 2]]
(12:44)>CPC(C)
%209 = 0[/CODE]

not sure why the last result and I know I forgot the empty set.[/QUOTE]

got it working:

[CODE]CPC(C)= Cp(C[1],C[2]);A=d;for(i=3,#C,Cp(A,C[i]);A=d);A;[/CODE] Now I'd still like it to be more proper output (using notation).

[CODE](14:31)>intersection(x,y) = C=[];for(i=1,#x,for(j=1,#y,if(x[i]==y[j],C=concat(C,x[i]))));C;
(14:33)>intersection([1,2,3,1],[2,3,1])
%259 = [1, 2, 3, 1][/CODE]

I think this is good.

[QUOTE=science_man_88;251292]got it working[/QUOTE]

Yes. You should be returning the result,though, not storing it in a variable. What if you needed to call it twice, like

dostuff(Cp(v1, v2), Cp(v3, v4))

?

[QUOTE=science_man_88;251292]Now I'd still like it to be more proper output (using notation).[/QUOTE]

Once you have a working function (as you do!), the right way -- what you seem to be doing -- is to have a second function that formats it "nicely". (An alternate, also good, method is to have an argument that, if set to a particular value, "nicens" the result -- though in this case that wouldn't be as good, I think.)

Maybe this:
[code]printVectorAsSet(v)={
if(type(v) != "t_VEC",
print1(v)
,
print1("{");
if (#v, printVectorAsSet(v[1]));
for(i=2,#v,
print1(", ");
printVectorAsSet(v[i])
);
print1("}")
)
}[/code]

So after fixing your function Cp to return the value, Cp([10,20,30],[1,2]) returns
[code][[10, 1], [10, 2], [20, 1], [20, 2], [30, 1], [30, 2]][/code]
and printVectorAsSet(Cp([10,20,30],[1,2])) writes
[code]{{10, 1}, {10, 2}, {20, 1}, {20, 2}, {30, 1}, {30, 2}}[/code]

[QUOTE=science_man_88;251298][CODE](14:31)>intersection(x,y) = C=[];for(i=1,#x,for(j=1,#y,if(x[i]==y[j],C=concat(C,x[i]))));C;
(14:33)>intersection([1,2,3,1],[2,3,1])
%259 = [1, 2, 3, 1][/CODE]

I think this is good.[/QUOTE]

Yep. It could be made faster, but that works fine.

I'll admit, I don't like the repetition. I suppose for both functions you could add vecsort(__,,8) in the appropriate places to fix that.

[QUOTE=CRGreathouse;251322]Yes. You should be returning the result,though, not storing it in a variable. What if you needed to call it twice, like

dostuff(Cp(v1, v2), Cp(v3, v4))

?



Once you have a working function (as you do!), the right way -- what you seem to be doing -- is to have a second function that formats it "nicely". (An alternate, also good, method is to have an argument that, if set to a particular value, "nicens" the result -- though in this case that wouldn't be as good, I think.)

Maybe this:
[code]printVectorAsSet(v)={
if(type(v) != "t_VEC",
print1(v)
,
print1("{");
if (#v, printVectorAsSet(v[1]));
for(i=2,#v,
print1(", ");
printVectorAsSet(v[i])
);
print1("}")
)
}[/code]

So after fixing your function Cp to return the value, Cp([10,20,30],[1,2]) returns
[code][[10, 1], [10, 2], [20, 1], [20, 2], [30, 1], [30, 2]][/code]
and printVectorAsSet(Cp([10,20,30],[1,2])) writes
[code]{{10, 1}, {10, 2}, {20, 1}, {20, 2}, {30, 1}, {30, 2}}[/code]



Yep. It could be made faster, but that works fine.

I'll admit, I don't like the repetition. I suppose for both functions you could add vecsort(__,,8) in the appropriate places to fix that.[/QUOTE]

I couldn't remember if I knew that or not lol.[CODE]intersection(x,y) = C=[];for(i=1,#x,for(j=1,#y,if(x[i]==y[j],C=concat(C,x[i]))));printVectorAsSet(vecsort(C,,8));[/CODE]

applying your make to a set function to the CPC gives back:

[CODE]{{{1, 1}, 1}, {{1, 1}, 2}, {{1, 2}, 1}, {{1, 2}, 2}, {{2, 1}, 1}, {{2, 1}, 2}, {{2, 2}, 1}, {{2, 2}, 2}, {{3, 1}, 1}, {{3, 1}, 2}, {{3, 2}, 1}, {{3, 2}, 2}, {{4, 1}, 1}, {{4, 1}, 2}, {{4, 2}, 1}, {{4, 2}, 2}, {{5, 1}, 1}, {{5, 1}, 2}, {{5, 2}, 1}, {{5, 2}, 2}, {{6, 1}, 1}, {{6, 1}, 2}, {{6, 2}, 1}, {{6, 2}, 2}, {{7, 1}, 1}, {{7, 1}, 2}, {{7, 2}, 1}, {{7, 2}, 2}, {{8, 1}, 1}, {{8, 1}, 2}, {{8, 2}, 1}, {{8, 2}, 2}}[/CODE]

but this comes from the fact that I never got each index to be one vector. on restart of PARI I got:

[QUOTE]*** syntax error, unexpected $end, expecting KPARROW or ',' or ')': ...=[],B=[])=if(#A!=0&&#B!=0,[/QUOTE] I don't understand why this is in a function called aredisjoint .

[CODE]aredisjoint(A=[],B=[]) =
if(#A!=0 && #B!=0,
a=0;
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0),
a=a+1;
if(a==#A*#B,
return(1)
)
)
)
),
return(0)
)
addhelp(aredisjoint,"a function to tell when 2 sets (vectors(which are actually sequences)) are considered disjoint.")[/CODE]

not sure why it's happened.

Looks like you missed a parenthesis.

[QUOTE=CRGreathouse;251346]Looks like you missed a parenthesis.[/QUOTE]

5 opening 5 closing for the for and if statements and I can't see anything else not closed.

[QUOTE=science_man_88;251329]

[CODE]aredisjoint[COLOR="Lime"]([/COLOR]A=[],B=[][COLOR="lime"])[/COLOR] =
if[COLOR="Red"]([/COLOR]#A!=0 && #B!=0,
a=0;
for[COLOR="DarkOrange"]([/COLOR]i=1,#A,
for[COLOR="Purple"]([/COLOR]j=1,#B,
if[COLOR="Blue"]([/COLOR]A[i]==B[j],
return[COLOR="Silver"]([/COLOR]0[COLOR="Silver"])[/COLOR],
a=a+1;
if[COLOR="DarkRed"]([/COLOR]a==#A*#B,
return[COLOR="Gray"]([/COLOR]1[COLOR="Gray"])[/COLOR]
[COLOR="DarkRed"])[/COLOR]
[COLOR="blue"])[/COLOR]
[COLOR="Purple"])[/COLOR]
[COLOR="DarkOrange"])[/COLOR],
return[COLOR="Magenta"]([/COLOR]0[COLOR="Magenta"])[/COLOR]
[COLOR="red"])[/COLOR]
addhelp(aredisjoint,"a function to tell when 2 sets (vectors(which are actually sequences)) are considered disjoint.")[/CODE][/QUOTE]

color coded