|
[QUOTE=science_man_88;259782][CODE]pieces=[N,B,K,R,Q,a,b,c,d,e,f,g,h];files=[a,b,c,d,e,f,g,h];rows=[1,2,3,4,5,6,7,8];captures=[x,""];for(piece=1,#pieces,for(capture=1,#captures,for(file=1,#files,for(row=1,#rows,print1(pieces[piece]captures[capture]files[file]rows[row]",")))))[/CODE]
once again I had it working ( i thought) but can't find my error this code is for: [url]http://www.mersenneforum.org/showthread.php?p=259371[/url] accidently defined c as a vector that's part of it.[/QUOTE] okay I've fixed one. |
[CODE]pick(n,k)=if(n>=k,return(n!/(k!*(n-k)!)),return((n+k-1)!/((k!)*(n-1)!)))[/CODE]
just thought I'd post one again. |
Faster version:
[code]pick(n,k)=binomial(if(n<k,n+k-1,n),k)[/code] Basically, your function is binomial(n+k-1, k) if n<k and binomial(n,k) otherwise. I could have written [code]pick(n,k)=if(n<k,binomial(n+k-1,k), binomial(n,k))[/code] but I chose to combine the similar parts and use just [code]pick(n,k)=binomial([COLOR="Red"]if(n<k,n+k-1,n)[/COLOR],k)[/code] When PARI can compute the binomial directly instead of indirectly through the factorial it can save a lot of steps. When you calculate [TEX]{9\choose2}=\frac{9!}{7!2!}=\frac{9\cdot8\cdot7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1}{(7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1)(2\cdot1)}[/TEX] you probably just cross out the 7 * ... * 1 on top and bottom. When you give PARI binomial(9,2) it can do the same; when you give it 9!/(7!*2!) it's forced to compute the large numbers and then divide. |
[code]binom(n,k)={ my(t=i=1);
if(2*k>n, k=n-k); while(i<=k, t=(t*n)/i; i++; n-- ); t }[/code]----------------------------- |
[QUOTE=CRGreathouse;262186]Faster version:
[code]pick(n,k)=binomial(if(n<k,n+k-1,n),k)[/code] Basically, your function is binomial(n+k-1, k) if n<k and binomial(n,k) otherwise. I could have written [code]pick(n,k)=if(n<k,binomial(n+k-1,k), binomial(n,k))[/code] but I chose to combine the similar parts and use just [code]pick(n,k)=binomial([COLOR="Red"]if(n<k,n+k-1,n)[/COLOR],k)[/code] When PARI can compute the binomial directly instead of indirectly through the factorial it can save a lot of steps. When you calculate [TEX]{9\choose2}=\frac{9!}{7!2!}=\frac{9\cdot8\cdot7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1}{(7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1)(2\cdot1)}[/TEX] you probably just cross out the 7 * ... * 1 on top and bottom. When you give PARI binomial(9,2) it can do the same; when you give it 9!/(7!*2!) it's forced to compute the large numbers and then divide.[/QUOTE] thanks:[CODE](19:14)>allocatemem(1) *** allocatemem: Warning: new stack size = 1024 (0.001 Mbytes). (19:15)>pick(11,1000000000) %269 = 275573207396384843474431809413623724923088088349707046408843005955895496034675000001 (19:15)>allocatemem(931987423) *** allocatemem: Warning: new stack size = 931987420 (888.812 Mbytes). (19:15)>pick1(11,1000000000) *** _!: the PARI stack overflows ! current stack size: 931987420 (888.812 Mbytes) [hint] you can increase GP stack with allocatemem()[/CODE] I changed my original to pick1 for this test. |
found something I find interesting now CRG:
[CODE](20:20)>for(a=1,38,print1(sum(x=1,a,pick(Me(x)%6,Me(x+1)%6)%4)%4",")) [COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,[/COLOR][/CODE] okay it works out this way for all primes it looks like. |
[QUOTE=science_man_88;262207]found something I find interesting now CRG:
[CODE](20:20)>for(a=1,38,print1(sum(x=1,a,pick(Me(x)%6,Me(x+1)%6)%4)%4",")) [COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,[/COLOR][/CODE] okay it works out this way for all primes it looks like.[/QUOTE] I made an extended version of MeVec to test this and so far the exponents for x=40 to 47( or guessing) fit this pattern. |
What happens if you reverse the order of the last 8?
|
[QUOTE=CRGreathouse;262228]What happens if you reverse the order of the last 8?[/QUOTE]
nothing I bet because late last night I realized why the pattern is set up it's because pick(a,b) with a = 1 or 5, b=1,5 will be 1 mod 4. so the sum is mostly a sum of ones ( which will only cycle mod anything) the 2 that start it 2 and 3 and 3 and 5 are 0 and 1 mod 4 so 0+1 is the start and +1 for each step after. so really all this states is that the first 3 exponents are used to set it up. |
[QUOTE=science_man_88;262207]found something I find interesting now CRG:
[CODE](20:20)>for(a=1,38,print1(sum(x=1,a,pick(Me(x)%6,Me(x+1)%6)%4)%4",")) [COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,[/COLOR][/CODE] okay it works out this way for all primes it looks like.[/QUOTE] I've got another pattern to figure out before revealing. |
[QUOTE=science_man_88;262323]I've got another pattern to figure out before revealing.[/QUOTE]
okay I gave in: [CODE](19:34)>for(x=1,47,print1(pick(prime(x),x^2-10)%x",")) 0,0,0,3,0,0,0,0,0,4,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,13,0,20,0,0,0,0,0,0,30,0,0,0,0,0,0,0,0,0,0,0,0, (19:34)>for(x=1,47,print1(pick(prime(x),x^2-5)%x",")) 0,0,2,0,0,2,0,0,0,4,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,0,0,25,0,0,0,13,0,0,28,0,0,0,0,0, (19:35)>for(x=1,47,print1(pick(MPE[x],x^2-10)%x",")) 0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15,0,0,0,0,0,0,0,0,0,0, (19:35)>for(x=1,47,print1(pick(MPE[x],x^2-5)%x",")) 0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, (19:35)>for(x=1,47,print1(MPE[x]<x^2-5",")) 0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, (19:47)>for(x=1,47,print1(MPE[x]<x^2-10",")) 0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,[/CODE] [QUOTE]MPE[x]!/(k!*(MPE[x]-k)!) (MPE[x]+k-1)!/((k!)*(MPE[x]-1)!)[/QUOTE] is what I consider important |
[QUOTE=science_man_88;262324]okay I gave in:
[CODE][COLOR="Red"](19:34)>for(x=1,47,print1(pick(prime(x),x^2-10)%x",")) 0,0,0,3,0,0,0,0,0,4,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,13,0,20,0,0,0,0,0,0,30,0,0,0,0,0,0,0,0,0,0,0,0, (19:34)>for(x=1,47,print1(pick(prime(x),x^2-5)%x",")) 0,0,2,0,0,2,0,0,0,4,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,0,0,25,0,0,0,13,0,0,28,0,0,0,0,0,[/COLOR] [COLOR="Lime"](19:35)>for(x=1,47,print1(pick(MPE[x],x^2-10)%x",")) 0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15,0,0,0,0,0,0,0,0,0,0, (19:35)>for(x=1,47,print1(pick(MPE[x],x^2-5)%x",")) 0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,[/COLOR] [COLOR="DarkOrange"](19:35)>for(x=1,47,print1(MPE[x]<x^2-5",")) 0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, (19:47)>for(x=1,47,print1(MPE[x]<x^2-10",")) 0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,[/COLOR][/CODE] is what I consider important[/QUOTE] the red is for proving not all primes follow the patterns 1 after the other. the green shows the patterns I think I have ( pick(Me(x),x^2-5) is 0 mod 2 and pick(Me(x),x^2-10) is 0 mod 3). the orange is to show which path of pick each Me(x) brings us through. |
sieve of atkins
not sure how to create this , I know I ignore 2,3,and 5, and I check mod 60 but I have no idea ( the pseudo-code in the Wikipedia article didn't help me as far as I can tell) otherwise.
|
Implementing the Atkin-Bernstein sieve efficiently is nontrivial. I wouldn't recommend it for you -- not until you've at least implemented
1. modular exponentiation by squaring 2. Miller-Rabin 3. the sieve of Eratosthenes 4. modular exponentiation with a sliding window 5. Karatsuba multiplication as warm-ups. You'll want to work in some suitably lowish-level language, certainly no higher than C++ (not gp, not Java, not Perl). |
[QUOTE=CRGreathouse;265046]3. the sieve of Eratosthenes[/QUOTE]
[CODE]Eratosthenes(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),for(y=1,#b,if(b[y]!=0&&b[y]%a[#a]==0,b[y]=0));for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));for(d=a[#a],#b,if(b[d]!=0,a=concat(a,b[d])));a[/CODE] is this close enough or do I have to find 2 prime as well ? |
[QUOTE=science_man_88;265058][CODE]Eratosthenes(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),for(y=1,#b,if(b[y]!=0&&b[y]%a[#a]==0,b[y]=0));for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));for(d=a[#a],#b,if(b[d]!=0,a=concat(a,b[d])));a[/CODE]
is this close enough or do I have to find 2 prime as well ?[/QUOTE] yeah I know I didn't optimize it lol. |
[QUOTE=science_man_88;265058]is this close enough[/QUOTE]
It's not the sieve of Eratosthenes. You do [code]b[y]!=0 && b[y]%a[#a]==0[/code] about [TEX]x\pi(\sqrt x)\approx x^{3/2}/\log x[/TEX] times, which is more operations right there than the entire SoE uses. Critically, though, you shouldn't be dividing in the inner loop at all -- you should just be adding and marking 0s. But again, GP isn't the right language for this -- as good as it is on high-level operations it's quite poor at low-level ones like b[y]=0. Consider what it actually does when faced with that statement: 1. Parse the bytecode into machine instructions 2. Go to the location of y and read its codewords. 3. Do some quick bitwise ops to check the type of y. 4. Do some quick bitwise ops to find the location of the length of y. 5. Check that the length of y is at most one machine word. 6. Check that the sign of y is positive. 7. Do some quick bitwise ops to find the value of the first word of the value of y. 8. Read this value into a register, treating it as an unsigned long. 9. Go to the location of b and read its codewords. 10. Do some quick bitwise ops to check the type of b. 11. Do some quick bitwise ops to find the location of the length of b. 12. Check that a is at most the length of b. 13. Do some quick arithmetic to calculate the address pointed to by b[y] 14. Check if the value pointed to by b[y] is at the end of the PARI stack. (?) 15. If so, move the stack pointer avma. (?) 16. Overwrite b[y] with a precomputed value that refers to gen_0. If the value was anything other than one of the special values, a new location would need to be created, the value written, and avma moved. But fortunately I don't think that has to happen in this case. Compare this to what C would have done: 1. Do some quick arithmetic to calculate the address pointed to by b[y] 2. Set this equal to 0. |
[QUOTE=CRGreathouse;265070]It's not the sieve of Eratosthenes. You do
[code]b[y]!=0 && b[y]%a[#a]==0[/code] about [TEX]x\pi(\sqrt x)\approxx^{3/2}/\log x[/TEX] times, which is more operations right there than the entire SoE uses. Critically, though, you shouldn't be dividing in the inner loop at all -- you should just be adding and marking 0s.[/QUOTE] it's not optimal but it does the exact thing the sieve does and gets the correct answer. |
[QUOTE=science_man_88;265071]it's not optimal but it does the exact thing the sieve does and gets the correct answer.[/QUOTE]
but I'll do it over . |
[QUOTE=science_man_88;265071]it's not optimal but it does the exact thing the sieve does and gets the correct answer.[/QUOTE]
No, actually it doesn't. The key operation in your function is the modular division % I quoted, but that operation is barely used at all in the SoE. In fact, it doesn't even need to use it at all! What you need to do is, for each prime p <= sqrt(x) is to go from p^2, p^2 + p, p^2 + 2p, ..., all the way up to x and mark each of these values as 0. By the way, I edited some information into #2248 explaining why GP isn't good at this. But still, I think you'll find that if you make the change I suggested above and replace the vector with a vectorsmall you can get the program to be at least 1000 times faster. :) |
[QUOTE=science_man_88;265073]but I'll do it over .[/QUOTE]
[CODE]erastosthenes(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),forstep(y=a[#a],#b,a[#a],b[y]=0);for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));a[/CODE] is this good enough or did I forget something else ? I know an optimization on this if it's correct though. |
[QUOTE=science_man_88;265076]is this good enough or did I forget something else ?[/QUOTE]
The core of it works properly. Unfortunately your copying between arrays isn't working quite right and it's giving the primes as 2, 3, 3, 3, 3, .... If your optimization involves not copying between arrays (very, very slow in GP) then why don't you do that and see if it fixes that problem. Two birds with one stone and all that. |
[QUOTE=CRGreathouse;265075]No, actually it doesn't. The key operation in your function is the modular division % I quoted, but that operation is barely used at all in the SoE. In fact, it doesn't even need to use it at all!
What you need to do is, for each prime p <= sqrt(x) is to go from p^2, p^2 + p, p^2 + 2p, ..., all the way up to x and mark each of these values as 0. By the way, I edited some information into #2248 explaining why GP isn't good at this. But still, I think you'll find that if you make the change I suggested above and replace the vector with a vectorsmall you can get the program to be at least 1000 times faster. :)[/QUOTE] is that so: [CODE] l(a[#a]>sqrt(x),forstep(y=a[#a],#b,a[#a],b[y]=0);for(z=1,#b,if(b[z]!= = (x)->a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),forstep(y=a[#a], 47)>Eratosthenes(10) * user interrupt after 33,172 ms. 48)>Eratosthenes3(x)= a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),f = (x)->a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),for(y=1,#b,if(b[ 48)>Eratosthenes3(10) = [2, 3, 5, 7] 48)>## * last result computed in 0 ms. 49)>[/CODE] actually my old code gave me a decrease of at least 1/33000 by the looks of these. |
[QUOTE=science_man_88;265081]is that so:[/QUOTE]
Yes, I had to interrupt the code to see what it was doing as it had fallen into an infinite loop. If your third version works, good for you. |
[QUOTE=CRGreathouse;265082]Yes, I had to interrupt the code to see what it was doing as it had fallen into an infinite loop. If your third version works, good for you.[/QUOTE]
no my first version worked better than my last one the loop isn't infinite it just has to recalculate a[#a] twice on the fly. that's basically the major difference I see. |
[QUOTE=science_man_88;265083]no my first version worked better than my last one the loop isn't infinite it just has to recalculate a[#a] twice on the fly. that's basically the major difference I see.[/QUOTE]
The version that you posted in post #2252 is an infinite loop, since each prime it adds to a is 3. For example, if I add [code]print1(b[z]" ");[/code] just before [code]if(b[z]!=0,a=concat(a,b[z]);break)[/code] I get [code]3 3 3 3 3 3 3 3 [3000 more 3s skipped] 3 3 3 3 3 3 3 ^C *** at top-level: Erastosthenes1(100) *** ^------------------- *** in function Erastosthenes1: ...int1(b[z]" ");if(b[z]!=0, *** a=concat(a,b[z]);bre *** ^-------------------- *** user interrupt after 330 ms. *** Break loop: type <Return> to continue; 'break' to go back to GP[/code] |
My browser's having a hard time rendering that uber-long string of 3's. Is it really necessary?
|
[QUOTE=CRGreathouse;265085]The version that you posted in post #2252 is an infinite loop, since each prime it adds to a is 3. For example, if I add
[code]print1(b[z]" ");[/code] just before [code]if(b[z]!=0,a=concat(a,b[z]);break)[/code] I get [code]3 3 3 3 3 3 3 3 3 (skipped 3's ) 3 3 3 ^C *** at top-level: Erastosthenes1(100) *** ^------------------- *** in function Erastosthenes1: ...int1(b[z]" ");if(b[z]!=0, *** a=concat(a,b[z]);bre *** ^-------------------- *** user interrupt after 330 ms. *** Break loop: type <Return> to continue; 'break' to go back to GP[/code][/QUOTE] I don't get any of that on my end except me stopping it. |
What do you get if you add the print statement I suggested?
|
[QUOTE=bsquared;265086]My browser's having a hard time rendering that uber-long string of 3's. Is it really necessary?[/QUOTE]
After my computer lovingly crafted each and every one of those 3s? I'm shocked by your audacity. |
[QUOTE=CRGreathouse;265089]After my computer lovingly crafted each and every one of those 3s? I'm shocked by your audacity.[/QUOTE]
They were a beautiful thing... I'm almost sorry to see them go. :max:But now they're back thanks to sm88's quote. :max: |
[QUOTE=CRGreathouse;265088]What do you get if you add the print statement I suggested?[/QUOTE]
I get a whole bunch of 3's but I don't see how it's possible outside of one thing: for the start of the forstep you need a if statement: [CODE](x)->a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),forstep(y=if(a[#a]==2,a[#a],a[#a]-2),#b,a[#a],b[y]=0);for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));a[/CODE] this happened because I made a have 2 in it to start I think. and because trying to eliminate multiples of three with a n+2 setup in the array made 3 in position 1 not 3 so it never gets eliminated. |
[QUOTE=science_man_88;265091]I get a whole bunch of 3's but I don't see how it's possible outside of one thing:
for the start of the forstep you need a if statement: [CODE](x)->a=[2];b=vector(x-2,n,n+2);until(a[#a]>sqrt(x),forstep(y=if(a[#a]==2,a[#a],a[#a]-2),#b,a[#a],b[y]=0);for(z=1,#b,if(b[z]!=0,a=concat(a,b[z]);break())));a[/CODE] this happened because I made a have 2 in it to start I think. and because trying to eliminate multiples of three with a n+2 setup in the array made 3 in position 1 not 3 so it never gets eliminated.[/QUOTE] now it's missing the part to add in the values greater than sqrt(x). |
Let's do this again. Use only one array and do something like this:
[code]b=vector(x,unused,1); b[1]=0; ptr=2; while(b[ptr] <= sqrt(x), \\ cross off numbers divisible by b[ptr] ptr++; while(b[ptr] == 0, ptr++); ); \\ Now b[n] is 1 if n is prime and 0 otherwise \\ Output the primes numberOfPrimes=sum(i=1,#b,b[i]); n=0; vector(numberOfPrimes,unused, while(b[n]==0,n++); n )[/code] |
[QUOTE=CRGreathouse;265096]Let's do this again. Use only one array and do something like this:
[code]b=vector(x,unused,1); b[1]=0; ptr=2; while(b[ptr] <= sqrt(x), \\ cross off numbers divisible by b[ptr] ptr++; while(b[ptr] == 0, ptr++); ); \\ Now b[n] is 1 if n is prime and 0 otherwise \\ Output the primes numberOfPrimes=sum(i=1,#b,b[i]); n=0; vector(numberOfPrimes,unused, while(b[n]==0,n++); n )[/code][/QUOTE] [CODE], ptr++););numberOfPrimes=sum(i=1,#b,b[i]);n=0;vector (x)->b=vector(x,unused,1);b[1]=0;ptr=2;while(b[ptr]<= >eratosene4(100) hile: array index (101) out of allowed range [1-100].[/CODE] |
I can't say why you got that without seeing your code, but my display code is slightly wrong (worked in my original version, failed when I translated it into a simpler version for you).
Try [code]n=0; vector(numberOfPrimes,unused, n++; while(b[n]==0,n++); n )[/code] |
[QUOTE=CRGreathouse;265121]I can't say why you got that without seeing your code, but my display code is slightly wrong (worked in my original version, failed when I translated it into a simpler version for you).
Try [code]n=0; vector(numberOfPrimes,unused, n++; while(b[n]==0,n++); n )[/code][/QUOTE] still gives the error I got. |
TF code
based on what i see at:
[url]http://www.mersenne.org/various/math.php[/url] I made a TF code a MTF code and a loop involving MTF to check prime exponents. [CODE](17:58)>?TF TF(x,y,B)=a=binary(x);b=0;for(c=1,#a,if(a[c]==1,if(c==1,b=((a[c]^2)*B)%y,b=((b^2)*B)%y),b=(b^2)%y));return(b);if(b==1,return(1),return(0)) (19:56)>?MTF MTF(p,d,e)=for(k=d,e,if((2*k*p+1)%8==7&&TF(p,2*k*p+1,2)==1,return(2*k*p+1);break()))[/CODE] I changed it once I made it after showing myself it could work with other bases ( I gave examples in the thread named that) as well as for higher exponents using lower binary values but a modified values to multiply all but that last 1 by ( which lets you test exponents k*p+1 , which I've been pointed to ( and kinda realized) that most of the supposed use of the altered version ( one situation I haven't coded yet) is useless due to the knowledge of 3 mod 4 and isprime(2p+1)) . the main problem with the alteration is finding a point at which the extra bits in the multiply , square and mod operations, are worth it for less of each. |
Wow. It'll take quite a bit of time to catch up to everything... How long was I out?
And I see familiar faces here. What have you all been working on? |
[QUOTE=3.14159;268161]Wow. It'll take quite a bit of time to catch up to everything... How long was I out?[/QUOTE]
Six hundred years. I didn't want to be the one to tell you this, but everyone you knew is now dead. Wait, wrong cue card. [QUOTE=3.14159;268161]And I see familiar faces here. What have you all been working on?[/QUOTE] In the PARI world I'm trying to get new features added: [url]http://pari.math.u-bordeaux.fr/archives/pari-dev-1107/msg00006.html[/url] Some of these I wrote already, some are just ideas. Outside of that, I'm going to the Rabin conference/party/thing on the 29th-30th. I've been fairly busy with the OEIS, too -- submissions increase faster than the editorial board. :smile: |
Aww, man.. I saw that PARI can't be updated for Windows anymore..
Or, rather, "don't have a reliable Windows machine to make them anymore" [quote]Warning: I no longer have a Windows machine on which to build these reliably. The old packages remain and may or may not cooperate nicely with your own Windows machine (they don't with mines). For the time being, I am no longer able to update them, so the 2.5.* stable branch is not available in this form. Sorry. [i]Karim Belabas.[/i][/quote] |
I'm going to try to make a Windows build for them -- I have access to a Windows machine -- but it's not going to be easy. There are lots of little difficulties (and some big ones).
For now either use the development version 2.4.x or get the SVN version (up to date but no high-res graphics or readline). |
I have a bit of a question..
Suppose I use this for(n=1,20, if(random(10) > 5, print(random(10), " is > 5") And I'd like it only to print #s > 5 (Or better yet, the first random(10))... How exactly would I go about in doing that? |
Store the random number in a variable, test the variable, then print it if it's larger than 5.
|
[QUOTE=CRGreathouse;268340]Store the random number in a variable, test the variable, then print it if it's larger than 5.[/QUOTE]
yeah that's how I did it in : [CODE]randomprime(n)=until(a>10^(n-1)&&isprime(a),a=random(10^(n)-1));return(a)[/CODE] kinda I think. |
[QUOTE=CRGreathouse;268340]Store the random number in a variable, test the variable, then print it if it's larger than 5.[/QUOTE]
That's a nice step.. But I'd need it to also continuously redefine the variable for me. What I'm trying to do is something similar to this: [code]for(n=1,20, print(random(10))) 4 7 8 0 0 6 5 4 3 5 4 7 6 1 3 0 9 4 5 7[/code] It printed random numbers from 0-9 (Why the hell not 1-10?? for random([b]10[/b])) I'd like to do the above, but excluding those numbers less than or equal to 5. If it's impossible in PARI/GP, feel free to point me to where it is possible to do so. Though I'd be very disappointed if it were actually impossible to do. |
[QUOTE=3.14159;268354]That's a nice step.. But I'd need it to also continuously redefine the variable for me.
What I'm trying to do is something similar to this: [code]for(n=1,20, print(random(10))) 4 7 8 0 0 6 5 4 3 5 4 7 6 1 3 0 9 4 5 7[/code] It printed random numbers from 0-9 (Why the hell not 1-10?? for random([b]10[/b])) I'd like to do the above, but excluding those numbers less than or equal to 5. If it's impossible in PARI/GP, feel free to point me to where it is possible to do so. Though I'd be very disappointed if it were actually impossible to do.[/QUOTE] if you allow an upper limit I'd try something like: [code]for(n=1,20,until(a>5,a=random(10));print(a","))[/CODE] [QUOTE]random({N=2^31}): random integer between 0 and N-1.[/QUOTE] that's why |
Thanks for the advice on the random function.
Okay, that's closer, but not quite yet it. What I would like to do is this: random(a)+b - (random(c)+d) Print 20 results, however, all results that are < 0 get the message, "no number here" instead. |
[QUOTE=3.14159;268373]Thanks for the advice on the random function.
Okay, that's closer, but not quite yet it. What I would like to do is this: random(a)+b - (random(c)+d) Print 20 results, however, all results that are < 0 get the message, "no number here" instead.[/QUOTE] so: [CODE]e=0;until(e>20,for(b=1,100,for(d=1,100,a=random();c=random();if((a+b)-(c+d)>0,print1((a+b)-(c+d)",");e=e+1;if(e>20,break(3))))))[/CODE] except more cleaned up ? |
Kind of like this:
[code](21:17) gp > for(n=1,20, print(random(40)+6 - (random(30)+11))) 15 -30 -27 9 7 0 21 2 -12 7 -3 19 -3 18 -34 0 5 -22 8 -22[/code] But all negative numbers/zeros get replaced by "no number here" |
[QUOTE=3.14159;268377]Kind of like this:
[code](21:17) gp > for(n=1,20, print(random(40)+6 - (random(30)+11))) 15 -30 -27 9 7 0 21 2 -12 7 -3 19 -3 18 -34 0 5 -22 8 -22[/code] But all negative numbers/zeros get replaced by "no number here"[/QUOTE] altering yours I got : [CODE]for(n=1,20,a=random(40);b=random(30);print(if(a+6 - (b+11)>0,a+6 - (b+11),"no number here")))[/CODE] |
Yes! I'll do some more testing
This looks exactly on target Thanks a million! :smile: |
[QUOTE=3.14159;268379]Yes! I'll do some more testing
This looks exactly on target Thanks a million! :smile:[/QUOTE] remember you can do a if within a print to choose what to print. |
Okay, thanks for the tip.
Also, I didn't understand how [code]for(n=1,20,a=random(40);b=random(30);print(if(a+6 - (b+11)>0,a+6 - (b+11),"no number here")))[/code] worked. |
It stores a and b, which are random numbers on [0, 39] and [0, 29] respectively. If a+6 - (b+11)>0, that is, a > b+5, it prints the difference a+6 - (b+11). Otherwise, it prints "no number here". This process is repeated a total of 20 times.
I would write this [code]{for(n=1,20, t=random(40)-random(30)-5; print(if(t > 0, t, "no number here")) )}[/code] |
Anywho, in the spirit of somewhat old times..
Primality testing 17025*2^18800+1 [N-1, Brillhart-Lehmer-Selfridge] Running N-1 test using base 7 Special modular reduction using all-complex FFT length 1536 on 17025*2^18800+1 Calling Brillhart-Lehmer-Selfridge with factored part 99.93% 17025*2^18800+1 is prime! (0.4559s+0.0015s) |
[QUOTE=science_man_88;268380]remember you can do a if within a print to choose what to print.[/QUOTE]
.. Sm88 .. Primecare .. See light .. Shut up & drive .. Lol |
[CODE]LLRM(X) = s=Mod(14,2^X-1);for(x=1,(X-3)/2,s=Mod(s^4,2^x-1)-4*Mod(s^2,2^X-1)+2);if(s==0,1,0)[/CODE]
yeah I need a better name than this but I used something CRG talked about to speed it up it originally only beat lucaslehmer(x) under 1279 now at over 9000 + it seems like 10 + times faster than it. |
[QUOTE=science_man_88;268689][CODE]LLRM(X) = s=Mod(14,2^X-1);for(x=1,(X-3)/2,s=Mod(s^4,2^x-1)-4*Mod(s^2,2^X-1)+2);if(s==0,1,0)[/CODE]
yeah I need a better name than this but I used something CRG talked about to speed it up it originally only beat lucaslehmer(x) under 1279 now at over 9000 + it seems like 10 + times faster than it.[/QUOTE] [CODE](16:16)>LLRM(3217) %40 = 1 (16:16)>## *** last result computed in 32 ms. (16:16)>lucaslehmer(3217) %41 = 1 (16:16)>## *** last result computed in 281 ms. (16:16)>LLRM(9689) %42 = 1 (16:17)>## *** last result computed in 485 ms. (16:17)>lucaslehmer(9689) %43 = 1 (16:17)>## *** last result computed in 4,922 ms. (16:17)>LLRM(21701) %44 = 1 (16:23)>## *** last result computed in 4,156 ms. (16:23)>lucaslehmer(21701) %45 = 1 (16:24)>## *** last result computed in 39,562 ms.[/CODE] found a typo in the code. never mind it's twice as slow. |
Pardon the question, but has anyone found the (10[sup]13[/sup])th prime? Or (10[sup]14[/sup])th prime?
I've only found up to 10[sup]12[/sup] Edit: Scratch that, found way more past 10[sup]13[/sup] and 10[sup]14.[/sup] |
prime(10^19) = 465675465116607065549.
|
[QUOTE=CRGreathouse;268994]prime(10^19) = 465675465116607065549.[/QUOTE]
[url]http://factordb.com/index.php?id=1100000000442493681[/url] i)..prp<309> |\| |_| |_ |_ |
..*9<73
[QUOTE=cmd;269010][url]http://factordb.com/index.php?id=1100000000442493681[/url]
i)..prp<309> |\| |_| |_ |_[/QUOTE] [url]http://factordb.com/index.php?id=1100000000442497462[/url] |) |= |\/| () [url]http://factordb.com/index.php?id=1100000000442493555[/url] Demo=andiamo (venice~it)s|_ang |
64-bit
I can't figure which one it is I upgraded to the 64-bit windows OS disk I had ( the 32 and 64 are in the same packaging).
|
[QUOTE=science_man_88;270501]I can't figure which one it is I upgraded to the 64-bit windows OS disk I had ( the 32 and 64 are in the same packaging).[/QUOTE]
Do you mean what version of Windows or what version of PARI? There isn't a released 64-bit Windows version of PARI/GP that has readline and the other nice features yet; there are serious difficulties building it. For Windows, right-click My Computer and choose Properties; this should tell you what version you're using. (Of course this is off-topic here.) |
1 Attachment(s)
a few codes I've sent through PM:
[CODE]ali(n)=a=n-1;b=[];for(x=0,floor(.5*a),if(x==0,if(sigma(a^2)-a^2==n,b=concat(b,a^2)),if(sigma((x)*(a-x))-((x)*(a-x))==n,b=concat(b,(x)*(a-x)))));b=vecsort(b,,8)[/CODE] checks up to floor(.5*a) in matching pairs multiplying them together and if the sigma checks go through puts them into b. sorts the vector at the end. [CODE]for(z=6,6,print(z);for(x=1,#ali(z),print("\t"ali(z)[x]);for(y=1,#ali(ali(z)[x]),print("\t\t"ali(ali(z)[x])[y]);for(h=1,#ali(ali(z)[x])[y],print("\t\t\t"ali(ali(ali(z)[x])[y])[h]);for(i=1,#ali(ali(ali(z)[x])[y])[h],print("\t\t\t\t"ali(ali(ali(ali(z)[x])[y])[h])[i]))))))[/CODE] z can be changed to any group of consecutive numbers this basically does the previous one on top of each other forming generation in later talk I believe I talk of it as aligen(n), please note these codes only find ones that have (n-1) split into 2 divisors. the attachment has been scanned since I want to double check it, it represents most that I have tried with these codes. |
[QUOTE=science_man_88;270856]a few codes I've sent through PM:
[CODE]ali(n)=a=n-1;b=[];for(x=0,floor(.5*a),if(x==0,if(sigma(a^2)-a^2==n,b=concat(b,a^2)),if(sigma((x)*(a-x))-((x)*(a-x))==n,b=concat(b,(x)*(a-x)))));b=vecsort(b,,8)[/CODE] checks up to floor(.5*a) in matching pairs multiplying them together and if the sigma checks go through puts them into b. sorts the vector at the end. [CODE]for(z=6,6,print(z);for(x=1,#ali(z),print("\t"ali(z)[x]);for(y=1,#ali(ali(z)[x]),print("\t\t"ali(ali(z)[x])[y]);for(h=1,#ali(ali(z)[x])[y],print("\t\t\t"ali(ali(ali(z)[x])[y])[h]);for(i=1,#ali(ali(ali(z)[x])[y])[h],print("\t\t\t\t"ali(ali(ali(ali(z)[x])[y])[h])[i]))))))[/CODE] z can be changed to any group of consecutive numbers this basically does the previous one on top of each other forming generation in later talk I believe I talk of it as aligen(n), please note these codes only find ones that have (n-1) split into 2 divisors. the attachment has been scanned since I want to double check it, it represents most that I have tried with these codes.[/QUOTE] made a remake of aligen: [CODE]aligen2(w,s)=c=[];for(z=w,s,for(h=1,#c,if(z==c[h],next(2)));print(z);for(x=1,#ali(z),if(ali(z)[x]<s&&ali(z)[x]>w,c=concat(c,ali(z)[x]));print("\t"ali(z)[x]);for(y=1,#ali(ali(z)[x]),print("\t\t"ali(ali(z)[x])[y]))))[/CODE] this checks the middle generation only ( my original thought did both generations ( I simplified to 2) back but was much slower I think), and puts the values between the limits in a vector that gets checked next z value etc. to check if it can be skipped without loss of information ( except maybe an extra generation, but for some the next one up is in the same limits and still does get checked) . |
connecting the loose ends for aliquot3.txt gives me:
[CODE]3 4 9 15 33 87 249 553 949 1273 247 1205 1673 3029 3893 4313 5129 9353 10229 11993 12629 13193 13529 13973 14453 14873 14933[/CODE] |
[QUOTE=science_man_88;270856]a few codes I've sent through PM:
[CODE]ali(n)=a=n-1;b=[];for(x=0,floor(.5*a),if(x==0,if(sigma(a^2)-a^2==n,b=concat(b,a^2)),if(sigma((x)*(a-x))-((x)*(a-x))==n,b=concat(b,(x)*(a-x)))));b=vecsort(b,,8)[/CODE] checks up to floor(.5*a) in matching pairs multiplying them together and if the sigma checks go through puts them into b. sorts the vector at the end. [CODE]for(z=6,6,print(z);for(x=1,#ali(z),print("\t"ali(z)[x]);for(y=1,#ali(ali(z)[x]),print("\t\t"ali(ali(z)[x])[y]);for(h=1,#ali(ali(z)[x])[y],print("\t\t\t"ali(ali(ali(z)[x])[y])[h]);for(i=1,#ali(ali(ali(z)[x])[y])[h],print("\t\t\t\t"ali(ali(ali(ali(z)[x])[y])[h])[i]))))))[/CODE] z can be changed to any group of consecutive numbers this basically does the previous one on top of each other forming generation in later talk I believe I talk of it as aligen(n), please note these codes only find ones that have (n-1) split into 2 divisors. the attachment has been scanned since I want to double check it, it represents most that I have tried with these codes.[/QUOTE] does anyone else know the way to extend the first code to check for 3 values to multiply that add up to a? I get the sigma check change would be something like: sigma(<new variable>*x*(a-(x+<new variable>)))==n but I want it to act correct so it checks sums of 2 numbers first then three. |
[QUOTE=science_man_88;271229]does anyone else know the way to extend the first code to check for 3 values to multiply that add up to a? I get the sigma check change would be something like: sigma(<new variable>*x*(a-(x+<new variable>)))==n but I want it to act correct so it checks sums of 2 numbers first then three.[/QUOTE]
wasteful ( doesn't solve for 2 as well because if so it throws 0 into sigma) method found! [CODE](n)->a=n-1;for(x=1,a,for(y=1,a,for(z=1,a,if(1+x+y+z==a&&sigma(x*y)-(x*y)==n,print(x*y),if(1+x+y+z==a&&sigma(x*z)-(x*z)==n,print(x*z),if(1+x+y+z==a&&sigma(y*z)-(y*z)==n,print(y*z),if(1+x+y+z==a&&sigma(x*y*z)-(x*y*z)==n,print(x*y*z))))))))[/CODE] |
Asking the pari experts....
Question 1:
I have a file with a million lines. (So what?) ASCII. Each line is a number (no other characters except 0..9) I have (in a separate variable) the number k. (not your business where did I get it from) I want to read in a variable the number on the line k on the file. How the hell I do that in pari/gp??? There is no "seek", "readline" or equivalents. The only possibility I found is to read all the file in a vector with "readvec" and take the component I need. What it would take (beside of "ages") few wheelbarrows of memory. Question 2: How the hell I can take the class "c" in "x=Mod(c,m)" as integer? Beside of the "ugly" construction "component(x,2)", which took me 2 days to find it (no specs in the help). They provide the member function ".mod" to take the m (like x.mod will return m), but no function to take the c. I figured later that I can apply a "component()" function to any internal structure, and from that it did not take me long to find out the the class c is in fact... the second component of the structure (the Mod structure is stored internally "viceversa"). In spite of the fact this question looks somehow tendentious, or picky, the reason is the fact that the only way to do a (somehow fast) "modexp" in pari is "a=Mod(x,m)^y" (which is taking modular reduction at each step, that is why is just "somehow" fast, and not really fast as in yafu), and latter if I want to compare "a==b" this will promote b from t_INT into t_INTMOD, killing all the process further. That is why the "component()" step has to be added. Any idea? (Question 1 is really important, otherwise I have to write a separate application, and use "extern()" to read the line and put it in a separate file, and feed that to pari. Question 2 is just a curiosity.) Thanks in advance. |
[QUOTE=science_man_88;271237]wasteful ( doesn't solve for 2 as well because if so it throws 0 into sigma) method found!
[CODE](n)->a=n-1;for(x=1,a,for(y=1,a,for(z=1,a,if(1+x+y+z==a&&sigma(x*y)-(x*y)==n,print(x*y),if(1+x+y+z==a&&sigma(x*z)-(x*z)==n,print(x*z),if(1+x+y+z==a&&sigma(y*z)-(y*z)==n,print(y*z),if(1+x+y+z==a&&sigma(x*y*z)-(x*y*z)==n,print(x*y*z))))))))[/CODE][/QUOTE] few changes got it working for both 2 and 3 number combinations: [CODE]a=n-1;for(x=0,a,for(y=1,a,for(z=1,a,if(x+y+z==a&&x!=0&&sigma(x*y)-(x*y)==n,print(x*y),if(x+y+z==a&&x!=0&&sigma(x*z)-(x*z)==n,print(x*z),if(x+y+z==a&&sigma(y*z)-(y*z)==n,print(y*z),if(x+y+z==a&&x!=0&&sigma(x*y*z)-(x*y*z)==n,print(x*y*z))))))))[/CODE] okay testing proves I stupidly forgot to print a^2 ( later will turn to a concat) and use a vecsort at the end. |
[QUOTE=LaurV;271261]Question 2:
How the hell I can take the class "c" in "x=Mod(c,m)" as integer? [/QUOTE] Look at lift() and centerlift() |
I'm out of practice, but getting line N out of a file is a one-liner in perl.
|
Pardon the question.. but similar to how (4!+1) = 5^2; (5!+1) = 11^2; (7!+1) = 71^2.. are there any examples of (n!+1) = x^3 (where n and x are integers)
|
[QUOTE=3.14159;275455]Pardon the question.. but similar to how (4!+1) = 5^2; (5!+1) = 11^2; (7!+1) = 71^2.. are there any examples of (n!+1) = x^3 (where n and x are integers)[/QUOTE]
Doubtful, since x would have to be divisible by only cubes of large primes which are very rare (sum converges). Just to be divisible by one prime cube where the prime is over 10,000 is a one-in-two billion chance, and most of that is on the lower end (which doesn't take care of the bulk of the number). |
[QUOTE=CRGreathouse;275461]Doubtful, since x would have to be divisible by only cubes of large primes which are very rare (sum converges). Just to be divisible by one prime cube where the prime is over 10,000 is a one-in-two billion chance, and most of that is on the lower end (which doesn't take care of the bulk of the number).[/QUOTE]
Ah, cool. Thanks for the heads-up. So would it be safe to say that (n!+1) = x^y has no solutions when y > 2? (n, x, and y all being integers, of course.) Lastly, reading through this thread.. can't help but feel a slight bit of embarrassment at what I've posted so far. |
making a new variable name via script ?
[QUOTE](x,y)->eval(Str(x))==y[/QUOTE]
the problem is when I retyped the name with a question mark it worked, But the flaw was copying in the name or other ways of getting the value behind the variable like typing it's name don't seem to work on my end. |
[QUOTE=science_man_88;275749]the problem is when I retyped the name with a question mark it worked, But the flaw was copying in the name or other ways of getting the value behind the variable like typing it's name don't seem to work on my end.[/QUOTE]
from a fresh window of gp.exe: [QUOTE](21:55)>var(x,y) = eval(Str(x))==return(y) %1 = (x,y)->eval(Str(x))==return(y) (21:55)>?weekend weekend: unknown identifier (21:55)>var(weekend,5) %2 = 5 (21:55)>?weekend weekend: user defined variable (21:56)>weekend %3 = weekend[/QUOTE] okay so retyping the equality with one equal sign can add the value after creation |
A PARI script to search for sociable numbers
I've been searching for sociable numbers in the 10^15 area for a while now, using the following PARI script:
[code] pS=readvec(pariStatus) a=pS[1] for(b=a,a+10^9,i=b;for(n=1,200,m=n;i=sigma(i)-i;if(i<b||i>b*20,break()));if(m>199,print(b);write("pariCyclesFound",b));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b))) [/code]The script uses the file pariStatus to keep track of how far along it is, such as: [code] 1000010274050000 [/code]In this example, the script will start searching at 1000010274050000, find a cycle and keep searching: [code] ... 1000010274060000 5.6430000000000000000000000000000000000 seconds/10000 1000010274063430 1000010274070000 5.5030000000000000000000000000000000000 seconds/10000 ... [/code]In truth. it doesn't guarantee it's a cycle, but if the sequence goes for 200 iterations without increasing to b*20, it is at least interesting. In the above example, that is a known amicable number. I do a lot more printing than I probably need to, but it shows activity. I'd prefer overwriting the pariStatus file, but couldn't figure out how to do that with PARI, so you can see the extern("rm pariStatus") command to get rid of the previous file prior to writing the current value. This is a time waster and is definitely hazardous if a power loss is timed just right, but I couldn't quite figure out a better way. For a while I had nine machines set up to search 11*10^14, 12*10^14, ... 19*10^14, respectively. Each was set to start on power-up, so if the AC was lost, they would begin again, within a few numbers (seconds) of where they were interrupted. Occasionally, I would record their progress in case an interruption occurred at a time that caused pariStatus to be corrupted or lost. I'm sure this could be improved upon and it may even have some flaw(s) I've not considered. So, I'm open to all critique... |
[QUOTE=EdH;275857]I've been searching for sociable numbers in the 10^15 area for a while now, using the following PARI script:
[code] pS=readvec(pariStatus) a=pS[1] for(b=a,a+10^9,i=b;for(n=1,200,m=n;i=sigma(i)-i;if(i<b||i>b*20,break()));if(m>199,print(b);write("pariCyclesFound",b));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b))) [/code]The script uses the file pariStatus to keep track of how far along it is, such as: [code] 1000010274050000 [/code]In this example, the script will start searching at 1000010274050000, find a cycle and keep searching: [code] ... 1000010274060000 5.6430000000000000000000000000000000000 seconds/10000 1000010274063430 1000010274070000 5.5030000000000000000000000000000000000 seconds/10000 ... [/code]In truth. it doesn't guarantee it's a cycle, but if the sequence goes for 200 iterations without increasing to b*20, it is at least interesting. In the above example, that is a known amicable number. I do a lot more printing than I probably need to, but it shows activity. I'd prefer overwriting the pariStatus file, but couldn't figure out how to do that with PARI, so you can see the extern("rm pariStatus") command to get rid of the previous file prior to writing the current value. This is a time waster and is definitely hazardous if a power loss is timed just right, but I couldn't quite figure out a better way. For a while I had nine machines set up to search 11*10^14, 12*10^14, ... 19*10^14, respectively. Each was set to start on power-up, so if the AC was lost, they would begin again, within a few numbers (seconds) of where they were interrupted. Occasionally, I would record their progress in case an interruption occurred at a time that caused pariStatus to be corrupted or lost. I'm sure this could be improved upon and it may even have some flaw(s) I've not considered. So, I'm open to all critique...[/QUOTE] add one more layer to this: [QUOTE]for(x=2,1000000,if((sigma(x)-x)!=x && sigma(sigma(x)-x)-(sigma(x)-x)==x ,print(x","sigma(x)-x)))[/QUOTE] and you have sociable numbers. this checks for amicable pairs ( just adding another sigma and a copy of what is there should work for 3). |
[QUOTE=science_man_88;275858]add one more layer to this:
and you have sociable numbers. this checks for amicable pairs ( just adding another sigma and a copy of what is there should work for 3).[/QUOTE] I'm not aware of any order 3 sociable numbers, but there are order 4, 5, 6, 8, 9 and 28, so I'd like my search to include at least that many levels. |
[QUOTE=EdH;275868]I'm not aware of any order 3 sociable numbers, but there are order 4, 5, 6, 8, 9 and 28, so I'd like my search to include at least that many levels.[/QUOTE]
I was thinking of making a loop to change the level of the sigma, but it might have to jump a lot to work. |
[QUOTE=EdH;275868]I'm not aware of any order 3 sociable numbers, but there are order 4, 5, 6, 8, 9 and 28, so I'd like my search to include at least that many levels.[/QUOTE]
sorry for posting again so soon people. the following is what I come to. [CODE]for(x=3,1000000,b=x;for(y=1,100,if(sigma(b)-b!=0,b=sigma(b)-b);if(b==x,print(x","y);break())))[/CODE] |
[QUOTE=science_man_88;275873]sorry for posting again so soon people. the following is what I come to.
[CODE]for(x=3,1000000,b=x;for(y=1,100,if(sigma(b)-b!=0,b=sigma(b)-b);if(b==x,print(x","y);break())))[/CODE][/QUOTE] Shouldn't be a problem with valid users posting "too soon!" I like your approach, but I need to get a better picture of the whole. I think there should be a break if sigma(b)-b falls below x. I'm going to need to add my tracking and check timing. Off to study, but maybe not too much tonight. Unfortunately, a preliminary trial with edits here didn't find the test cycle from the previous example...:sad:, but it ran faster:smile:... Thanks... |
[QUOTE=EdH;275919]Shouldn't be a problem with valid users posting "too soon!"
I like your approach, but I need to get a better picture of the whole. I think there should be a break if sigma(b)-b falls below x. I'm going to need to add my tracking and check timing. Off to study, but maybe not too much tonight. Unfortunately, a preliminary trial with edits here didn't find the test cycle from the previous example...:sad:, but it ran faster:smile:... Thanks...[/QUOTE] how long does it take to cycle because I did checks up to 1000000 in length and it can't find it also I checked by hand that's not the full cycle if it is one: [CODE](09:00)>sigma(1000010274070000)-1000010274070000 %154 = 1421114600505088 (09:00)>sigma(%)-% %155 = 1415563371597376 (09:01)>sigma(%)-% %156 = 1515576001115456 (09:01)>sigma(%)-% %157 = 2075868183012544 (09:01)>sigma(%)-% %158 = 2096174602984256 (09:01)>sigma(%)-% %159 = 2063421874812754 (09:01)>sigma(%)-% %160 = 1043754255742634 (09:01)>sigma(%)-% %161 = 545264138006614 (09:01)>sigma(%)-% %162 = 336474062378666 (09:01)>sigma(%)-% %163 = 195012689531734 (09:01)>sigma(%)-% %164 = 97506344765870 (09:01)>sigma(%)-% %165 = 78039653764690 (09:01)>sigma(%)-% %166 = 67385070876590 (09:01)>sigma(%)-% %167 = 53908056701290 (09:01)>sigma(%)-% %168 = 56557293070742 (09:01)>sigma(%)-% %169 = 31204023763258[/CODE] |
in fact:
[url]http://factordb.com/sequences.php?se=1&eff=2&aq=1000010274060000&action=all&fr=0&to=100[/url] proves it's not a cycle. [CODE](10:15)>for(x=1000010274060000,1000010274060000,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x || isprime(b),print(x","y);break()))) 1000010274060000,228 (10:15)>## *** last result computed in 1,235 ms. (10:15)>for(x=1000010274060000,1000010274060000,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x,print(x","y);break()))) (10:22)>## *** last result computed in 12,859 ms.[/CODE] also finds an end with alteration. |
I think you have confused my checkpoint, 1000010274060000, with the amicable number, [URL="http://factordb.com/sequences.php?se=1&eff=2&aq=1000010274063430&action=all&fr=0&to=100"]1000010274063430[/URL].:smile:
[code][85486] 45 Needham 2004 1000010274063430=2*7^2*5*31*53507*1230371 1161476072988602=2*7^2*17*251*571*727*6691[/code]are known amicable ([URL]http://amicable.adsl.dk/aliquot/c2/c2_16.txt[/URL]): [code] ? sigma(1000010274063430)-1000010274063430 %1 = 1161476072988602 ? sigma(1161476072988602)-1161476072988602 %2 = 1000010274063430 ? [/code]Note: In the following examples, the lines with time are checkpoints and the line without a time is the sociable result. The new code now works, but when I add in all the print stuff, it bogs down. Although a slight improvement, this new method does not result in a significant savings: [code] for(b=1000010274050000,1000010274080000,i=b;for(n=1,100,i=sigma(i)-i;if(i<b||i>b*20,break());if(i==b,print(b);write("pariCyclesFound",b);break()));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b))) [/code]gives: [code] ... 1000010274060000 5.6750000000000000000000000000000000000 seconds/10000 1000010274063430 1000010274070000 5.4970000000000000000000000000000000000 seconds/10000 ... [/code]compared to: [code] for(b=1000010274050000,1000010274080000,i=b;for(n=1,200,m=n;i=sigma(i)-i;if(i<b||i>b*20,break()));if(m>199,print(b);write("pariCyclesFound",b));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b)))[/code]which gives: [code] ... 1000010274060000 5.6780000000000000000000000000000000000 seconds/10000 1000010274063430 1000010274070000 5.5030000000000000000000000000000000000 seconds/10000 ... [/code]Actually, the loss is probably not directly the print statements, but the checks I do for whether to print. I also have received some info via PM that may help with my pariStatus reading/writing. I have to study it a bit. But, I have some other issues to attend to for now. Thanks everyone for the assistance... |
[QUOTE=EdH;275997]I think you have confused my checkpoint, 1000010274060000, with the amicable number, [URL="http://factordb.com/sequences.php?se=1&eff=2&aq=1000010274063430&action=all&fr=0&to=100"]1000010274063430[/URL].:smile:
[/QUOTE] yeah probably but you should still find it with mine: [CODE](10:29)>for(x=1000010274063430,1000010274063430,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x,print(x","y);break()))) 1000010274063430,2[/CODE] |
[QUOTE=EdH;275857]II'd prefer overwriting the pariStatus file, but couldn't figure out how to do that with PARI, so you can see the extern("rm pariStatus") command to get rid of the previous file prior to writing the current value. This is a time waster and is definitely hazardous if a power loss is timed just right, but I couldn't quite figure out a better way.[/QUOTE]
GP doesn't have much in the way of I/O. I suggest [code]extern("rm oldStatus; mv pariStatus oldStatus"); write("pariStatus", ...); extern("rm oldStatus");[/code] |
[QUOTE=CRGreathouse;276052]GP doesn't have much in the way of I/O. I suggest
[code]extern("rm oldStatus; mv pariStatus oldStatus"); write("pariStatus", ...); extern("rm oldStatus");[/code][/QUOTE] Thanks. That would keep me from possibly losing the value. Is the first "rm oldStatus" just to "be sure" before the mv? |
[QUOTE=science_man_88;276004]yeah probably but you should still find it with mine:
[CODE](10:29)>for(x=1000010274063430,1000010274063430,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x,print(x","y);break()))) 1000010274063430,2[/CODE][/QUOTE] In computer terms, yours seems to be taking an eternity to search from 1000010274060000 to 1000010274080000: [code] ? for(x=1000010274060000,1000010274080000,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x,print(x","y);break()))) ^C *** sigma: user interrupt after 7mn, 32,978 ms. [/code]whereas, mine (without the extra prints) is done in a few seconds: [code] ? for(b=100001027406,100001027408,c=b*10000;for(d=c,c+10000,i=d;f=d*20;for(e=1,5,i=sigma(i)-i;if(i<d,break());if(i>f,break());if(i==d,print(d);break())))) 1000010274063430 time = 8,999 ms. [/code] |
[CODE]lucas2(t,z,p,q,x)= v=[t,z];for(y=3,x,v=concat(v,p*v[y-1]-q*v[y-2]));v=vector(x,i,v[i])[/CODE]
though I tried using it like:[CODE]for(z=-1000,1000,for(a=-1000,1000,if(lucas2(2,3,z,a,13)==[2,3,5,7,13,17,19,31,61,89,107,127,521],print(z","a))))[/CODE] and got told "bug in" I'm guessing something like the segmentation fault: please report errors. |
error testing
how can it be done in PARI because if possible it may allow me to test something , I'm currently working on related to Mersenne prime exponents. because laurv talked of code like:[CODE] sqrt(Mod(2,2^61-1)) [/CODE] well so far in testing manually I see a pattern ( gulp) to the result and I might be able to tell just by the first calculation and this kinda had to do with trying the lucas-lehmer test backwards. I think that if p is a mersenne prime exponent the result isn't and error on the second level down by adding 2 ( I am positive it's likely known though).
|
[QUOTE=science_man_88;280037]how can it be done in PARI because if possible it may allow me to test something , I'm currently working on related to Mersenne prime exponents. because laurv talked of code like:[CODE] sqrt(Mod(2,2^61-1)) [/CODE] well so far in testing manually I see a pattern ( gulp) to the result and I might be able to tell just by the first calculation and this kinda had to do with trying the lucas-lehmer test backwards. I think that if p is a mersenne prime exponent the result isn't and error on the second level down by adding 2 ( I am positive it's likely known though).[/QUOTE]
doh I found it I think lol trap() still not sure how to use it. |
The general format is
trap(, code, what to do if the code fails) as in [code]trap(,1/0,"don't divide by 0")[/code] Technically you can put the error type before the first comma to catch only certain types of errors, but usually you don't want to do this. |
[QUOTE=CRGreathouse;280126]The general format is
trap(, code, what to do if the code fails) as in [code]trap(,1/0,"don't divide by 0")[/code] Technically you can put the error type before the first comma to catch only certain types of errors, but usually you don't want to do this.[/QUOTE] thanks I already proved myself wrong without it but I still want to see if it even semi fits: [CODE](17:34)>b=2;forprime(p=2281,2281,for(y=1,p-1,b=sqrt(Mod(b,1<<p-1));trap(,break(),);if(y==p-1,print(p)));b=2) *** sqrt: user interrupt after 13,750 ms. *** break: bug in PARI/GP (Segmentation Fault), please report *** Break loop (type 'break' or Control-d to go back to GP)[/CODE] |
| All times are UTC. The time now is 23:20. |
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.