![]() |
[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] |
| All times are UTC. The time now is 23:05. |
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.