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)

CRGreathouse 2010-12-14 00:05

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

kar_bon 2010-12-14 01:53

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.

axn 2010-12-14 03:18

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

CRGreathouse 2010-12-14 04:19

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

science_man_88 2010-12-14 20:07

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

CRGreathouse 2010-12-14 20:15

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.

science_man_88 2010-12-14 21:12

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

science_man_88 2010-12-14 21:28

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)

science_man_88 2010-12-14 22:13

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

axn 2010-12-14 22:21

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

kar_bon 2010-12-14 22:57

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]


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

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