mersenneforum.org Linear Congruence order 4 sequence not testing in PFGW?
 User Name Remember Me? Password
 Register FAQ Search Today's Posts Mark Forums Read

 2017-07-07, 22:22 #1 carpetpool     "Sam" Nov 2016 13716 Posts Linear Congruence order 4 sequence not testing in PFGW? Is there a way to make pfgw test recurrence relations of order 3 or 4 at all? The Linear() doesn't work for these types of sequences, can anyone find a different way please? H. C. Williams, R. K. Guy, Some fourth-order linear divisibility sequences, Intl. J. Number Theory 7 (5) (2011) 1255-1277, P1=4, P2=3, Q=2. What is in bold requires three P, P, Q parameters, not two which is why the Linear() function doesn't work. I would like to run pfgw to test larger terms of A215458. So far (let a(n) be the nth term of this sequence) a(n) is prime for n = 3, 5, 7, 11, 17, 19, 23, 101, 107, 109, 113, 163, 283, 311, 331, 347, 359, 701, 1153, 1597, 1621... (up to 2000). At least does anyone care to verify these? Thanks! Last fiddled with by carpetpool on 2017-07-07 at 22:24
 2017-07-08, 03:47 #2 CRGreathouse     Aug 2006 10111001000112 Posts I can verify those terms (at least as generating probable primes), and can find the next few as 2063, 2437, 2909, 3319, 6011, and 12829. I don't know of a good way to do this in pfgw. It would be really convenient to have a general linear recurrence operator, not just the current restricted version. (I recognize it wouldn't have nice computational and divisibility properties.) I hope someone else knows more than me? A bad solution would be to generate terms in gp and then run them through pfgw, like so: Code: a(n)=([0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; -4, 8, -7, 4]^n*[0; 1; 4; 7])[1, 1] P=prod(i=1,3e5,prime(i)); forprime(n=2e4,1e5, t=a(n); if(gcd(t,P)==1, write("delme.pfgw", t))) in gp then Code: pfgw64 -f100 -s300000 delme.pfgw on the command line. Last fiddled with by CRGreathouse on 2017-07-08 at 03:57
 2017-07-08, 04:35 #3 carpetpool     "Sam" Nov 2016 13716 Posts And to clarify this is a strong-divisiblity sequence so if a(n) is prime, so is n. I think the Linear() function is still a problem for some recurrence relations, but as for this specific sequence has the closed form: a(n) = (2^n - (1/2 - (i sqrt(7))/2)^n - (1/2 + (i sqrt(7))/2)^n + 1)/2 So if PFGW can handle square roots, and imaginary numbers, this will take care of this sequence only (and many others too). But last time I tried these on PFGW, it evaluated these expressions differently and gave different results. Anyone know more? Last fiddled with by carpetpool on 2017-07-08 at 04:36
 2017-07-08, 04:39 #4 CRGreathouse     Aug 2006 5,923 Posts OK, I ran a quick test of the (bad) code I presented, except I bumped the trial division to 1e7. Interestingly (at least to me) it found 46147 and 46471 but no others to 60k.
 2017-07-08, 04:55 #5 carpetpool     "Sam" Nov 2016 31110 Posts Thanks for those! Are they in the correct order up to n = 60k? a(n) is prime for n = 3, 5, 7, 11, 17, 19, 23, 101, 107, 109, 113, 163, 283, 311, 331, 347, 359, 701, 1153, 1597, 1621, 2063, 2437, 2909, 3319, 6011, 12829, 46147, 46471...
2017-07-08, 14:58   #6
rogue

"Mark"
Apr 2003
Between here and the

2·2,953 Posts

Quote:
 Originally Posted by carpetpool Is there a way to make pfgw test recurrence relations of order 3 or 4 at all? The Linear() doesn't work for these types of sequences, can anyone find a different way please? H. C. Williams, R. K. Guy, Some fourth-order linear divisibility sequences, Intl. J. Number Theory 7 (5) (2011) 1255-1277, P1=4, P2=3, Q=2. What is in bold requires three P, P, Q parameters, not two which is why the Linear() function doesn't work. I would like to run pfgw to test larger terms of A215458. So far (let a(n) be the nth term of this sequence) a(n) is prime for n = 3, 5, 7, 11, 17, 19, 23, 101, 107, 109, 113, 163, 283, 311, 331, 347, 359, 701, 1153, 1597, 1621... (up to 2000). At least does anyone care to verify these? Thanks!
Yes. Use pfgw's scripting language. It isn't too difficult once you get familiar with it. I used it recently for a similar task. If interested, send me a PM and I'll give you an example.

 2017-07-08, 15:01 #7 carpetpool     "Sam" Nov 2016 311 Posts I didn't see your post there but how easy is it to transcribe PARI/GP to Pfgw language? Last fiddled with by carpetpool on 2017-07-08 at 15:02
2017-07-08, 15:08   #8
CRGreathouse

Aug 2006

5,923 Posts

Quote:
 Originally Posted by carpetpool Thanks for those! Are they in the correct order up to n = 60k? a(n) is prime for n = 3, 5, 7, 11, 17, 19, 23, 101, 107, 109, 113, 163, 283, 311, 331, 347, 359, 701, 1153, 1597, 1621, 2063, 2437, 2909, 3319, 6011, 12829, 46147, 46471...
Yes. Just one more up to 100k: 3, 5, 7, 11, 17, 19, 23, 101, 107, 109, 113, 163, 283, 311, 331, 347, 359, 701, 1153, 1597, 1621, 2063, 2437, 2909, 3319, 6011, 12829, 46147, 46471, 74219.

 2017-07-09, 15:05 #9 carpetpool     "Sam" Nov 2016 311 Posts So the code that was used to find these primes, if there was a way to automatically generating it to input one number in the file at a time, and have pfgw test the file (automatically of course) and create a new one once the original file has been processed. Then, of course the process repeats again. Any idea as to how this might go? I am using the code @CRGreathouse gave me: Code: binaryprod(v)=if(#v<4, factorback(v), my(t=#v\2); binaryprod(v[1..t])*binaryprod(v[t+1..#v])); primorial(x)=binaryprod(primes(primepi(x))); P=primorial(3e5); a(n)=([0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; -4, 8, -7, 4]^n*[0; 1; 4; 7])[1, 1] P=prod(i=1,500,prime(i)); for(n=20000,1e5, t=a(n); if(gcd(t,P)==1, write("delme.pfgw", t))) pfgw64 -f100 delme.pfgw Also mention would like to pick up testing at n = 100k to 500k. Thanks!
 2017-07-13, 01:52 #10 CRGreathouse     Aug 2006 5,923 Posts OK, that got a little jumbled -- you compute the primorial using the new code, then using the old code. Code: binaryprod(v)=if(#v<4, factorback(v), my(t=#v\2); binaryprod(v[1..t])*binaryprod(v[t+1..#v])); primorial(x)=binaryprod(primes(primepi(x))); P=primorial(1e7); \\ go big or go home a(n)=([0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; -4, 8, -7, 4]^n*[0; 1; 4; 7])[1, 1] for(n=10^5,10^6, t=a(n); if(gcd(t,P)==1, write("delme.pfgw", t))) You could try something like this if you wanted to run them one-by-one in pfgw, though it may be slow: Code: binaryprod(v)=if(#v<4, factorback(v), my(t=#v\2); binaryprod(v[1..t])*binaryprod(v[t+1..#v])); primorial(x)=binaryprod(primes(primepi(x))); P=primorial(1e7); a(n)=([0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; -4, 8, -7, 4]^n*[0; 1; 4; 7])[1, 1] for(n=10^5,10^6, t=a(n); if(gcd(t,P)==1, system("del delme.pfgw"); write("delme.pfgw", t); system("pfgw64 -s10000000 delme.pfgw"))) This code assumes you're on Windows; if you're on linux, you probably know enough to make the required changes.
 2017-07-13, 05:38 #11 carpetpool     "Sam" Nov 2016 311 Posts That code works, meanwhile I formed one that works with ALL linear recurrence relations. The only issue is that unlike the code (above) this one does not have a minimum value. As you will notice, the maximum value is approximately the same size as the maximum term. For this specific sequence involved in the code, a(n) is approximate to 2^(n-1), so to test terms up to a(n), use the maximum limit as 2^n. Any suggestions to the script? Code: SCRIPT DIM anm1, 8 DIM anm2, 7 DIM anm3, 4 DIM anm4, 1 DIM anmax, 2^100 DIM an OPENFILEAPP prp_file,prp.txt OPENFILEAPP prime_file,prime.txt LABEL start SET an,4*anm1-7*anm2+8*anm3-4*anm4 IF (an > anmax) THEN GOTO the_end IF (an==0) THEN GOTO next_iteration IF (an>0) THEN FACTORIZE an ELSE FACTORIZE -an IF (FACTORFOUND > 1) THEN GOTO next_iteration PRP an IF (ISPRP) THEN WRITE prp_file,an LABEL next_iteration SET anm4, anm3 SET anm3, anm2 SET anm2, anm1 SET anm1, an GOTO start LABEL the_end CLOSEFILE prp_file CLOSEFILE prime_file END

 Thread Tools

 Similar Threads Thread Thread Starter Forum Replies Last Post CRGreathouse Software 2 2017-06-14 16:05 Trejack Information & Answers 2 2016-04-30 05:30 storm5510 Math 27 2009-09-22 23:14 Unregistered Information & Answers 4 2006-10-04 22:38 abiessuunreg Miscellaneous Math 3 2005-03-07 21:03

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

Thu Sep 24 05:36:40 UTC 2020 up 14 days, 2:47, 0 users, load averages: 1.09, 1.26, 1.25

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

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.