PARI's commands

[CRGreathouse]

Originally Posted by kar_bon

Evaluate the Sum of p^p, p prime, up to m-th prime:

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
};

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

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)
};

[kar_bon]

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:

addhelp(WritePPowerPFile, "WritePPowerPFile(m): Write PPowerP(m) into a file.");
WritePPowerPFile(m)={
  write("PPowerP/PPowerP"m".txt",SumPPowerP(m));
};

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

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);
};

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.

[axn]

Originally Posted by CRGreathouse

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

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

[CRGreathouse]

Originally Posted by axn

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

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

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

[science_man_88]

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)));for(b=1,#x,if(b=<#x && x[b]==0,x=vector((b-1),n,x[n]);break()))

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.

[CRGreathouse]

I get a different error message:

 ***   syntax error, unexpected '<': for(b=1,#x,if(b=<#x&&x[b]==0,x=vecto

The black section of your code includes a =< which should be <=. Once I change that it works for me.

[science_man_88]

Originally Posted by CRGreathouse

I get a different error message:

 ***   syntax error, unexpected '<': for(b=1,#x,if(b=<#x&&x[b]==0,x=vecto

The black section of your code includes a =< which should be <=. Once I change that it works for me.

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.

[science_man_88]

I named my code to show you:

(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

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

[science_man_88]

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

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

[axn]

Originally Posted by science_man_88

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

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

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.

[kar_bon]

Like this:

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);
};