|
I gave up on matrices.
[CODE]multiplication(i,j=i)= for(i=1,i,for(j=1,j,print1(i*j"\t"));print(""))[/CODE] |
[QUOTE=science_man_88;227950]1) I couldn't get [] to work for matrices[/QUOTE]
You'll need to define the matrix first. You tried -- but you only made it a 0x0 matrix, you need an i x j. try matrix(i,j,unused1, unused2, 0) or even matrix(i,j,x,y,x*y) which avoids the need for the loops. [QUOTE=science_man_88;227950]2) I wanted the whole table but I couldn't get it all to print[/QUOTE] Just return the matrix. [code]foo(i,j)=matrix(i,j,x,y,x*y)[/code] |
[QUOTE=science_man_88;227951]I gave up on matrices.
[CODE]multiplication(i,j=i)= for(i=1,i,for(j=1,j,print1(i*j"\t"));print(""))[/CODE][/QUOTE] That looks fine, except that your variables are bumping into each other: you have the "max j" and the "loop j" and likewise for i. If you rename them this should work nicely. |
[QUOTE=CRGreathouse;227952]You'll need to define the matrix first. You tried -- but you only made it a 0x0 matrix, you need an i x j. try matrix(i,j,unused1, unused2, 0) or even matrix(i,j,x,y,x*y) which avoids the need for the loops.
Just return the matrix. [code]foo(i,j)=matrix(i,j,x,y,x*y)[/code][/QUOTE] I tried your idea: [CODE][COLOR="Red"]*** too many arguments: ...[/COLOR]plication(i,j=i)=x=Mat(i,j,x,y,x*y);for(i=1,i [COLOR="Red"]^--------------------[/COLOR][/CODE] |
[CODE]mt(i,j=i)= for(x=1,i,for(n=1,j,print1(x*n"\t"));print(""))[/CODE]
is my best idea so far shortens the name and if a variable is missing it knows they are equal so to get to a square we only need one input. the trouble comes when the width is more than one page. seems to print to 20,20 perfectly any higher next thing that shows up is the last one. |
[QUOTE=science_man_88;227955]I tried your idea:
[CODE][COLOR="Red"]*** too many arguments: ...[/COLOR]plication(i,j=i)=x=Mat(i,j,x,y,x*y);for(i=1,i [COLOR="Red"]^--------------------[/COLOR][/CODE][/QUOTE] You used "Mat" instead of "matrix". |
[QUOTE=science_man_88;227956][CODE]mt(i,j=i)= for(x=1,i,for(n=1,j,print1(x*n"\t"));print(""))[/CODE]
is my best idea so far shortens the name and if a variable is missing it knows they are equal so to get to a square we only need one input.[/QUOTE] Yep, works pretty well. [QUOTE=science_man_88;227956]the trouble comes when the width is more than one page. seems to print to 20,20 perfectly any higher next thing that shows up is the last one.[/QUOTE] That's a limitation of the terminal. If you want bigger tables, or better-looking tables, I recommend using write() and write1() to create an HTML file which you can display in a browser. Something like this: [code]table(h,w)={ write("output.html", "<!doctype html>"); write("output.html", "<html><head><meta charset='utf-8'><title>Multiplication</title></head>"); write("output.html", "<body><table>"); for(i=1,h, write1("output.html", "<tr>"); for(j=1,w, write1("output.html", "<td>", i*j, "</td>") ); write("output.html", "</tr>"); ); write("output.html", "</table></body></html>"); };[/code] |
[QUOTE=CRGreathouse;227963]Yep, works pretty well.
That's a limitation of the terminal. If you want bigger tables, or better-looking tables, I recommend using write() and write1() to create an HTML file which you can display in a browser. Something like this: [code]table(h,w)={ write("output.html", "<!doctype html>"); write("output.html", "<html><head><meta charset='utf-8'><title>Multiplication</title></head>"); write("output.html", "<body><table>"); for(i=1,h, write1("output.html", "<tr>"); for(j=1,w, write1("output.html", "<td>", i*j, "</td>") ); write("output.html", "</tr>"); ); write("output.html", "</table></body></html>"); };[/code][/QUOTE] I may make it more efficient I tested it twice in a row and it was going about 452.5 bits/s = baud ? each time. |
The inefficiency comes from opening and closing the file for every write1() statement. You might get more speed by creating a string with the whole row, then writing that; I'm not sure.
|
[quote=CRGreathouse]The inefficiency comes from opening and closing the file for every write1() statement. You might get more speed by creating a string with the whole row, then writing that; I'm not sure.[/quote]
What's sm88 trying to do? |
[CODE]mt(i,j=i)=for(n=1,i,v=vector(j,x,n*x);write1("E:\\output.txt",v);write("E:\\output.txt","\n"))[/CODE]
this seems nice to me though it's a .txt (still able to open in a browser(well dah)) and writes at a rate of about 5881.4 bits/s = baud ? . For someone with so many programming books I sure suck lol, lets see: Dynamic html (useful to remember Pari's escape characters), Advanced MS-DOS (don't know if i have much use for it as it's the 1986 edition I think). Advanced Programming in the Unix environment (didn't check title until later and did nothing),the art of assembly language 2nd edition,and I think "Beginning Programming all in one desk reference for dummies" is the last one. |
[QUOTE=3.14159;227980]What's sm88 trying to do?[/QUOTE]
He's writing a multiplication table -- at first to the terminal, now to a file. I don't know his underlying motivation, but since this teaches valuable tools (matrices, file I/O) I'm not really concerned. |
[QUOTE=science_man_88;227983][CODE]mt(i,j=i)=for(n=1,i,v=vector(j,x,n*x);write1("E:\\output.txt",v);write("E:\\output.txt","\n"))[/CODE]
this seems nice to me though it's a .txt (still able to open in a browser(well dah)) and writes at a rate of about 5881.4 bits/s = baud ?[/QUOTE] Yes, it's faster since it's writing a whole row at a time. Pari's file I/O performance is bad -- actually there's been an outstanding 'todo' on that particular feature for a while (some suggested a "t_FILE" type to avoid this performance penalty). [QUOTE=science_man_88;227983]For someone with so many programming books I sure suck lol, lets see: Dynamic html (useful to remember Pari's escape characters), Advanced MS-DOS (don't know if i have much use for it as it's the 1986 edition I think). Advanced Programming in the Unix environment (didn't check title until later and did nothing),the art of assembly language 2nd edition,and I think "Beginning Programming all in one desk reference for dummies" is the last one.[/QUOTE] You'll learn programming faster than I did when I started, since I was younger than you are now, when I started. Just remember that while books can help, you can [i]only[/i] learn programming by doing. You can read a dozen books about how to read a bike, the physics of bike-riding, etc. without learning to ride. :grin: |
[QUOTE=3.14159;227878]A180362![/QUOTE]
Ok. When [url=http://oeis.org/classic/A180362]A180362[/url] shows up, I'll edit it to add [code]%S A180362 5,13,109,163,257,271,379,433,487,541,769,3329,7681,7937,9473,10753,11777, %T A180362 12289,13313,14081,14593,15361,17921,18433,19457,22273,23041,23297,25601, %U A180362 26113,26881,30977,31489,32257,36097,36353,37501,37633,37889,39937,40193 %H A180362 Charles R Greathouse IV, <a href="b180362.txt">Table of n, a(n) for n = 1..10000</a>[/code] and whatever else it may need. Did you already put the cross-refs in? |
it's up, and doesn't look like he did.
|
[QUOTE=science_man_88]it's up, and doesn't look like he did.
[/QUOTE] Meh, I made the corrections myself. All Charles needs to do is add the b-file. Or, if I think of a working script, I can do this myself! |
[QUOTE=3.14159;228026]Meh, I made the corrections myself. All Charles needs to do is add the b-file. Or, if I think of a working script, I can do this myself![/QUOTE]
pretty sure I got it working if i declare it in 2 parts [CODE](11:08) gp > w(x) = x^x %85 = (x)->x^x (11:09) gp > w(4) %86 = 256 (11:10) gp > w(4) %87 = 256 (11:10) gp > kbb(a, x, m) = { if(x<w(m), for(n=a, x, if(isprime(n*m^m+1), print1(n*m^m+1, ", ")) )); } (11:10) gp > kbb(1,3,4) 257, 769,[/CODE] |
[CODE]kbb(a, x, m) = if(x<m^m, for(n=a, x, if(isprime(n*m^m+1), print1(n*m^m+1, ", ")) ))[/CODE]
this one works on my end to give the same result in one function. |
if you want them in order pass them to to a vector then use vecsort() to sort them oh and if you want to loop until a specific n you could add loops then call this script to print for all x<m^m etc.
|
[QUOTE=science_man_88]if you want them in order pass them to to a vector then use vecsort() to sort them oh and if you want to loop until a specific n you could add loops then call this script to print for all x<m^m etc.
[/QUOTE] Yes, yes, I'm working on the sequence's b-file. Unless Charles decided to hinder my efforts by sending it first. |
[QUOTE=3.14159;228031]Yes, yes, I'm working on the sequence's b-file.[/QUOTE]
I was just giving suggestions as I know what's likely to be asked. |
[QUOTE=science_man_88]I was just giving suggestions as I know what's likely to be asked.
[/QUOTE] I give up. I'll let Charles work on that. I'll work on the sequence, with no restriction. In this case, b > 1. (This means the set of 4n + 1 primes is a subset of this sequence.) |
Here's what I submitted:
[code]%I A180362 %S A180362 5,13,109,163,257,271,379,433,487,541,769,3329,7681,7937,9473,10753,11777, %T A180362 12289,13313,14081,14593,15361,17921,18433,19457,22273,23041,23297,25601, %U A180362 26113,26881,30977,31489,32257,36097,36353,37501,37633,37889,39937,40193 %N A180362 Primes of the form k * n^n + 1 with k < n^n. %C A180362 A result of Heath-Brown shows, on the GRH, that this sequence is infinite; can this be proved unconditionally? The averaged result of Bombieri-Friedlander-Iwaniec does not seem to be strong enough. %H A180362 Charles R Greathouse IV, <a href="b180362.txt">Table of n, a(n) for n = 1..10000</a> %H A180362 D. R. Heath-Brown, "<a href="http://eprints.maths.ox.ac.uk/166/">Zero-free regions for Dirichlet L-functions, and the least prime in an arithmetic progression</a>", Proceedings of the London Mathematical Society 64:3 (1992), pp. 265-338. %F A180362 k * n^n + 1, where k < n^n. %e A180362 a(4) = 109, because 4 * 3^3 + 1 = 109, which is prime, and 4 < 27. %o A180362 (PARI) isA180362(n)=my(b=2);while(b^b<n,if(n%(b^b)==1 && n < b^(2*b), return(isprime(n)));b++);0 %K A180362 nonn,new %O A180362 1,1 %A A180362 Kevin Batista (kevin762401(AT)yahoo.com), Aug 30 2010[/code] |
[QUOTE=3.14159;228031]Yes, yes, I'm working on the sequence's b-file. Unless Charles decided to hinder my efforts by sending it first.[/QUOTE]
I told you in #1246 I already had it... |
[QUOTE=CRGreathouse]I told you in #1246 I already had it...
[/QUOTE] Not to worry, I only submitted the corrections. I never sent the b-file. Also: Made sequence A175768 for the sequence with no restriction on k, meaning the 4n + 1 primes are a subset of the sequence, and my guess on how much of the sequence are 4n + 1 primes = 85 to 99 percent. I have the members of this sequence up to 4k. |
[QUOTE=3.14159;228033]I'll work on the sequence, with no restriction. In this case, b > 1. (This means the set of 4n + 1 primes is a subset of this sequence.)[/QUOTE]
When you submit that sequence, you may wish to include its relative density in the primes (you can probably work this out, or else find my post where I give it) and cross-references (at least to A180362). |
[QUOTE=CRGreathouse]When you submit that sequence, you may wish to include its relative density in the primes (you can probably work this out, or else find my post where I give it) and cross-references (at least to A180362).
[/QUOTE] I made a crossref to A180362. I sent it. It's A175768. |
I assume that when you say "no restrictions" you mean "k > 1".
[QUOTE=3.14159;228037]Made sequence A175768 for the sequence with no restriction on k, meaning the 4n + 1 primes are a subset of the sequence, and my guess on how much of the sequence are 4n + 1 primes = 85 to 99 percent.[/QUOTE] It's 94.7028...%. I prefer to express this as a relative density of 0.527967... in the primes. [QUOTE=3.14159;228037]I have the members of this sequence up to 4k.[/QUOTE] To 10,000 it's just the 4n+1 primes plus 530 exceptions: [code]163,271,379,487,811,919,1459,1567,1783,1999,2539,2647,2971,3079,3187,3511,3727,3943,4051,4159,4483,4591,5023,5347,5563,5779,6211,6427,6967,7507,7723,8263,8803,9127,9343,9883,10099,10531,10639,11071,11287,11503,11719,11827,12043,12583,12799,12907,13339,13879,14419,14851,15391,15607,15823,16363,16903,17011,17443,17551,17659,18199,18307,18523,19387,19603,19819,19927,20143,20359,20899,21871,22303,23059,23167,23599,24247,24571,25111,25219,25759,25867,26083,26407,26731,26839,26947,27271,27487,27919,28027,28351,29863,30187,30403,30727,31051,31159,31267,31699,32563,32779,32887,33211,33427,33751,33967,34183,34939,35803,35911,36343,36451,36559,37423,37747,37963,38287,38611,39043,39367,39799,40123,40231,40771,40879,41203,41851,41959,42283,42391,42499,43579,44119,44983,45307,45523,45631,46171,46279,46819,47143,47251,47791,48871,49411,49627,49843,50383,50599,50707,50923,51031,51679,51787,52543,53299,53407,53623,53731,54163,54919,55243,55351,56431,57727,57943,58699,59023,59239,59671,59779,59887,60103,60427,60859,61291,61507,61723,62047,63127,63559,63667,64747,65071,65179,65287,65719,65827,67231,67339,67447,67987,68311,68743,69067,69499,69931,70039,70687,71011,71119,71443,71551,71983,72091,72307,72739,73063,73387,73819,74143,75223,75979,76303,76519,77167,77383,77491,78031,78139,78571,78787,79111,79867,80191,80407,81163,81703,81919,82351,82567,82891,83431,85159,86131,86239,87103,87211,87427,87643,87751,89371,90019,90127,91099,91423,91639,92179,92503,93151,94447,94771,95203,95311,95419,95527,95959,96823,96931,97039,97579,97687,98011,98227,98443,99523,100279,100927,101359,101467,102547,102763,102871,103087,103843,103951,104059,104383,104491,104707,105031,106219,106543,106759,106867,107839,108271,108379,109567,109891,110323,110431,110647,110863,111187,112807,113023,113131,113779,114319,114643,114859,114967,115183,115399,115831,116047,116371,116803,116911,117127,117883,117991,118423,118747,118751,119179,119503,119611,119827,120691,120907,121123,121447,122203,122527,122743,123499,123931,124147,124363,124471,125119,125551,125659,126199,126307,126631,126739,127711,127819,128467,128683,129223,129439,129763,130087,130303,130411,130843,131059,131251,131707,132247,132679,133327,133543,134191,134731,134839,134947,135271,136027,136351,136999,138079,138403,138511,138727,139267,139483,139591,140779,141319,142183,142939,143263,144451,144667,144883,145207,145423,145531,145963,146719,147151,147583,147799,148123,148339,148663,149419,150067,150607,151471,151579,151687,151903,152443,152767,153523,153739,154279,154387,154927,155251,156007,156979,157303,157411,157519,157627,157951,158923,159463,159571,159787,160651,161407,161731,161839,161947,162703,163027,163243,163351,163567,164431,165079,165511,166807,167023,167779,167887,168211,168643,169399,169831,170047,170263,170371,171559,172423,173827,174259,174367,174583,174799,174907,175447,175663,176419,177283,177823,178039,178903,179119,179659,180307,180847,181063,181387,181603,181711,181927,182467,182899,183439,183763,183871,183979,184087,184627,184843,185167,185491,185599,185707,185923,186247,186679,187003,187111,187219,187651,188299,188407,189271,190027,190243,190783,190891,192187,193051,193751,194239,194671,195103,195319,195427,195751,195967,196291,196831,197371,197479,197803,198127,199207,201151,201907,202231,202339,202879,202987,203311,203419,204067,204931,206251,206551,207199,207307,207523,207847,208279,208387,208927,209359[/code] To 1e5 there are 5313 exceptions, but b-files are usually limited to 10,000 terms. |
[QUOTE=Charles]It's 94.7028...%. I prefer to express this as a relative density of 0.527967... in the primes.
[/QUOTE] Excellent. Slightly less than 95% of them, or about 18 parts in 19 are 4n + 1 primes. How did you manage to work out the percentage? [QUOTE=Charles]To 10,000 it's just the 4n+1 primes plus 530 exceptions: [/QUOTE] How do you sort them in numerical order, anyway? Also: I assume you listed the set of primes that were not 4n + 1s. |
[QUOTE=3.14159;228048]I assume you listed the set of primes that were not 4n + 1s.[/QUOTE]
That is what I wrote, yes. [QUOTE=3.14159;228048]Excellent. Slightly less than 95% of them, or about 18 parts in 19 are 4n + 1 primes. How did you manage to work out the percentage?[/QUOTE] Clearly all primes of the form k * 4^4 + 1 are also of the form k1 * 2^2 + 1, and in general so are all primes of the form k * (2a)^(2a) + 1. Generalizing further, we need only consider primes; numbers of the form k * 25^25 + 1 are also of the form k1 * 5^5 + 1. This avoids the general use of inclusion-exclusion (which requires a great number of terms) since powers of distinct primes are coprime. So you have a simple, rapidly-converging infinite product over primes. The first dozen primes are sufficient to calculate 60 decimal places. |
[QUOTE=Charles]Clearly all primes of the form k * 4^4 + 1 are also of the form k1 * 2^2 + 1, and in general so are all primes of the form k * (2a)^(2a) + 1. Generalizing further, we need only consider primes; numbers of the form k * 25^25 + 1 are also of the form k1 * 5^5 + 1.
[/QUOTE] I think you meant, k * (2[sup]a[/sup]) ^ (2[sup]a[/sup]) + 1. [QUOTE=Charles]This avoids the general use of inclusion-exclusion (which requires a great number of terms) since powers of distinct primes are coprime. So you have a simple, rapidly-converging infinite product over primes. The first dozen primes are sufficient to calculate 60 decimal places. [/QUOTE] The series can be expressed as.. ? |
[QUOTE=3.14159;228055]I think you meant, k * (2[sup]a[/sup]) ^ (2[sup]a[/sup]) + 1.[/QUOTE]
No. If that was all I knew I couldn't avoid terms like 6^6 which would cause a combinatorial explosion in the calculation. [QUOTE=3.14159;228055]The series can be expressed as.. ?[/QUOTE] See if you can come up with it! 1 in 4 are 1 mod 4, 1 in 27 are 1 mod 27, so since 4 and 27 are relatively prime, 1/4 + 1/27 - 1/(4 * 27) are 1 mod 4 or 1 mod 27. |
[QUOTE=Charles]No. If that was all I knew I couldn't avoid terms like 6^6 which would cause a combinatorial explosion in the calculation.
[/QUOTE] Oh, I see my error. Pardon. [QUOTE=Charles]See if you can come up with it! 1 in 4 are 1 mod 4, 1 in 27 are 1 mod 27, so since 4 and 27 are relatively prime, 1/4 + 1/27 - 1/(4 * 27) are 1 mod 4 or 1 mod 27. [/QUOTE] 1/4 + 1/27 + 1/3125 + 1/823543 + 1/285311670611 + 1/302875106592253 + 1/827240261886336764177 + 1/1978419655660313589123979 + 1/20880467999847912034355032910567 + .... - 1/(4 * 27 * 3125 * 823543 * 285311670611 * 302875106592253 * 827240261886336764177 * 1978419655660313589123979 * 20880467999847912034355032910567 * ....) ?? |
That will be too large, since it's only removing one congruence class mod H (the huge product). You removed 1 mod H, but you also need to remove 1 + 4*27, 1 + 2*4*27, ... mod H.
I suggest looking at it as a product rather than a sum, it's easier to express. |
[QUOTE=Charles]I suggest looking at it as a product rather than a sum, it's easier to express.
[/QUOTE] Speaking of the sums.. The sums converge to this: [URL="http://www.research.att.com/~njas/sequences/A094289"]A094289[/URL] Here are the first 300 digits: [code]0.287358251306224179736418045878932206955908802685881709299499368947089329278688939770209124280029090055929603180132199757677833174625274203928613500682866624372279071764951496386358568820464783694988950221338310990369641738444509170337274489547045606825482008978904241753401587644678759939089840746020[/code] When extended to all integers: [code]1.29128599706266354040728259059560054149861936827452231731000244513694453876523445555881704112942970898499507092481543054841048741928486419757916355594791369649697415687802079972917794827300902564923055072096663812846701205368574597870300127789412928825355177022238337531934574925996777964830084954911[/code] |
Factoring each for a bit: Former: 2 ^ 2 * 5 * 7 * 62039 * 4685224417 * c283
Latter: 43 * 199 * 353 * 285101 * 546233 * 72659134783 * 74772299267 * c260 |
[QUOTE=3.14159;228063]Speaking of the sums.. The sums converge to this: [URL="http://www.research.att.com/~njas/sequences/A094289"]A094289[/URL] [/QUOTE]
Right. That's too big, though -- it double-counts 4 residues mod 108, for example, making it too large by at least 4/108 = 3.7...%. |
[QUOTE=Charles]Right. That's too big, though -- it double-counts 4 residues mod 108, for example, making it too large by at least 4/108 = 3.7...%.
[/QUOTE] Right.. [QUOTE=Charles]That will be too large, since it's only removing one congruence class mod H (the huge product). You removed 1 mod H, but you also need to remove 1 + 4*27, 1 + 2*4*27, ... mod H. [/QUOTE] 109, 217, 325, 433, 541, 649, 757, ... ? |
Also: Here is the number, such that n[sup]n[/sup] = 2 :
1.559610469462369349970388768765002993284883511843... Any simple way to compute that? I had to do trial and error to get the first 10 or so digits. |
It's log(2)/W(log(2)). I can calculate it with
[code]\p 1000 log(2)/LambertW(log(2))[/code] using my program [code]LambertW(x)={ my(e,t,w,ep); if(x <= 0, if (x == 0, return (0)); if (x == -exp(-1), return (-1)); if (x < -exp(-1), error("LambertW: "x" is out of range, exiting.")) ); \\ Initial approximation for iteration if (x < 1, ep = sqrt(2*exp(1)*x+2); \\ Using ep as a temporary variable w = ep * (1 - ep * (1/3 + 11/72*ep)) - 1 , w = log(x)); if (x>3, w = w - log(w)); t = 1; ep = eps() * (1 + abs(w)); \\ ep = epsilon while (abs(t) > ep, \\ Main (Halley) loop, cubic convergence e = exp(w); t = w*e - x; t = t/(e*(w+1) - .5*(w+2)*t/(w+1)); w = w - t ); w }; addhelp(LambertW,"Primary branch of Lambert's W function. Finds an L >= -1 such that L * exp(L) = x, where x >= -1/e."); \\ 2 ^ -(decimal_precision * lg(10) - 1) eps()={ precision(2.>>(32 * ceil(precision(1.) * 0.103810252965230073370947482)), 9) };[/code] Be careful, though; you'll want some guard digits. If you want 10,000 digits I'd calculate it with \p 10100 just in case. |
Excellent. I have turned it into the basis of my program, wroot. (Finds x such that x[sup]x[/sup] = n.)
Based on w(n) being n[sup]n[/sup]. |
[QUOTE=3.14159;228083]Excellent. I have turned it into the basis of my program, wroot. (Finds x such that x[sup]x[/sup] = n.)[/QUOTE]
Presumably wroot(n)=log(n)/LambertW(log(n)). This can be calculated directly slightly more quickly by other means, but I'm sure for your purposes very little speed is needed. |
[QUOTE=CRGreathouse]Presumably wroot(n)=log(n)/LambertW(log(n)). This can be calculated directly slightly more quickly by other means, but I'm sure for your purposes very little speed is needed.
[/QUOTE] This is for when I wish to find an n-digit k-b-b, and it will show what the value of b is. |
[QUOTE=3.14159;228086]This is for when I wish to find an n-digit k-b-b, and it will show what the value of b is.[/QUOTE]
I gathered as much. You might also consider writing a script that, given b and a max-k value, estimates the number of primes. :smile: |
[QUOTE=CRGreathouse]You might also consider writing a script that, given b and a max-k value, estimates the number of primes.
[/QUOTE] Step 1. log(x) Step 2. Divide by prime factors. (p)/(p-1) Step 3. Divide result by number of candidates Step 4. These are the amount of primes expected. Step 5. Divide step 2 number by sieve efficiency to learn chances of finding a prime after sieving to x. |
[QUOTE=3.14159;228088]Step 1. log(x)
Step 2. Divide by prime factors. (p)/(p-1) Step 3. Divide result by number of candidates Step 4. These are the amount of primes expected. Step 5. Divide step 2 number by sieve efficiency to learn chances of finding a prime after sieving to x.[/QUOTE] Step 1: log(N) Step 2: Divide step 1 by (1.781*log(x)) to learn chances of finding a prime after sieving to x. |
[QUOTE=Axn]Step 1: log(N)
Step 2: Divide step 1 by (1.781*log(x)) to learn chances of finding a prime after sieving to x.[/QUOTE] Something about something else comes to mind at this point.. Oh, I know! Sieving for k * 28657[sup]28657[/sup] + 1 and k * 2[sup]328750[/sup] + 1. (Approx. 127700 and 98900 digits, respectively.) So, that'll register for items 1 and 15. If I prove them prime using Proth.exe, add item 20 to the list as well. |
Back to the restricted sequence A180362:
Not that it's worth calculating -- checking bases is surely faster for tractable values -- but as a curios, the restriction is [TEX]\left\lceil\frac{\log\sqrt n}{\operatorname{W}(\log\sqrt n)}\right\rceil\le b\le\left\lfloor\frac{\log n}{\operatorname{W}(\log n)}\right\rfloor[/TEX] |
yeah on the kbb(a,x,m) idea I apparently forgot I was doing something in Pari and I didn't realize it even though by the time it gives it should of been at 3 in the morning when I was asleep this started.
a=1;for(m=1,10,for(x=1,m^m-1,kbb(a,x,m))) *** print1: user interrupt after 10h, 40mn, 54,500 ms. (13:11) gp > a=1;for(m=1,10,for(x=1,m^m-1,kbb(a,x,m))) oh and the highest I see in the last few lines before I stopped it: 65390961287 which is over 65 billion stand corrected: 98418329759 over 98 billion found. |
I'm currently looking for a record k-b-b. (b = 28657), and that will also be the largest prime I ever found,. :smile: (Expected digits: 127737)
|
I have now commenced testing for k * 2[sup]328750[/sup] + 1 and k * 28657[sup]28657[/sup] + 1.
Here are the odds for each: For the base 2 search: 1 in 4283. For the base 28657 search: 1 in 5812. Time taken for each test: 3 minutes for base 2. Expected time for base 2 completion: Approx. 9 days. Expected time for base 28657 completion: .. 81 days. |
[QUOTE=3.14159;228275]I have now commenced testing for k * 2[sup]328750[/sup] + 1 and k * 28657[sup]28657[/sup] + 1.
Here are the odds for each: For the base 2 search: 1 in 4283. For the base 28657 search: 1 in 5812. Time taken for each test: 3 minutes for base 2. Expected time for base 2 completion: Approx. 9 days. Expected time for base 28657 completion: .. 81 days.[/QUOTE] well you're k range for your expected length is 10^114477-10^114478 how many candidates in that range ? |
[QUOTE=science_man_88]well you're k range for your expected length is 10^114477-10^114478 how many candidates in that range ?
[/QUOTE] For b = 2, there are approximately 4880 candidates. For b = 28657, there are approx. 9289 candidates. Therefore, I will expect only one prime for each, in 9 and 81 days, respectively. |
I'll report back in November for b = 28657.
|
[CODE](13:45) gp > ifprime(x,p)=if(isprime(x),p)
%11 = (x,p)->if(isprime(x),p) (13:47) gp > ifprime(100,print(100)) 100[/CODE] I was trying something like this but I can't get it to work. why ? is my best question. |
[QUOTE=science_man_88;228683]
I was trying something like this but I can't get it to work. why ? is my best question.[/QUOTE] [CODE]ifprime(p, s)=if(isprime(p), eval(s)); ifprime(100, "print(p)"); ifprime(101, "print(p)"); 101 [/CODE] |
[QUOTE=axn;228687][CODE]ifprime(p, s)=if(isprime(p), eval(s));
ifprime(100, "print(p)"); ifprime(101, "print(p)"); 101 [/CODE][/QUOTE] I should of known this by now lol thank you though in a way this is useless but it is useful in the fact that it could create shorter code when typing it in. as by the look of it it saves 10 characters I was also thinking of one where you can do a && part in the if I think I'll try that next. |
[CODE](14:28) gp > ifprime(p,u,s)=if(isprime(p) && eval(u), eval(s))
%42 = (p,u,s)->if(isprime(p)&&eval(u),eval(s)) (14:29) gp > ifprime(100,"p%9==1","print(100)") (14:29) gp > ifprime(101,"p%9==1","print(100)") (14:29) gp > ifprime(101,"p%9==2","print(100)") 100 (14:29) gp > ifprime(101,"p%9==2","print(p)") 101[/CODE] thank you axn |
[QUOTE=science_man_88;228690]I should of known this by now lol thank you though in a way this is useless but it is useful in the fact that it could create shorter code when typing it in. as by the look of it it saves 10 characters I was also thinking of one where you can do a && part in the if I think I'll try that next.[/QUOTE]
okay maybe not but it's at least as small as what it equals. |
Primes p such that p*10[sup]2[/sup] + 1 is also prime: 701, 1301, 1901, 3701, 6101, 6701, 7901, 10301, 13901, 15101, 16301, 19301, 21101, 22901, 27701, 42101, 46301, 52301, 54101, 60101, 64301, 70901, 72701, 81101, 82301, 87701, 88301, 93701, 102101, 112901, 115301, 117101, 123701, 132701, 138101, 139901, 144701, 148301, 157901, 159701, 165701, 166301, 174101, 178301, 186701, 193301, 198701, 201101, 201701, 217901, 220301, 226901, 237701, 243701, 247301, 250301, 252101, 253901, 270701, 285101, 291701, etc...
Submitted to OEIS. |
[QUOTE=3.14159;228701]Primes p such that p*10[sup]2[/sup] + 1 is also prime: 701, 1301, 1901, 3701, 6101, 6701, 7901, 10301, 13901, 15101, 16301, 19301, 21101, 22901, 27701, 42101, 46301, 52301, 54101, 60101, 64301, 70901, 72701, 81101, 82301, 87701, 88301, 93701, 102101, 112901, 115301, 117101, 123701, 132701, 138101, 139901, 144701, 148301, 157901, 159701, 165701, 166301, 174101, 178301, 186701, 193301, 198701, 201101, 201701, 217901, 220301, 226901, 237701, 243701, 247301, 250301, 252101, 253901, 270701, 285101, 291701, etc...
There should only be two sequences: One for p * 10[sup]n[/sup] + 1, when n is any integer, and for 100.[/QUOTE] so in other words p*10[SUP]n[/SUP]+1 and p*10[SUP]100[/SUP] +1 ? or p*10[SUP]n[/SUP]+1 and p*100[SUP]n[/SUP]+1 ? |
Former.
|
[QUOTE=3.14159;227558]There are no redundancies.
Not every square number is a fourth power. Fourth powers are a subset of squares. Nevermind. This is getting me nowhere. Fine, you win. You can keep the smug expression on your face.[/QUOTE] what CRG was saying is: [COLOR="Red"]1[/COLOR],4,9,[COLOR="Red"]16[/COLOR],25,36,49,64,[COLOR="Red"]81[/COLOR],100,121,144,169,196,225,[COLOR="Red"]256[/COLOR],.... [COLOR="Red"]1,16,81,256,[/COLOR]........ 1) n[SUP]2[/SUP] 2) n[SUP]4[/SUP] Red) the members in both sequences see how all the n[SUP]4[/SUP] are Red ? since it includes all fourth powers only then it can be stated: within the set of fourth powers, instead of in the set of squares and the set of fourth powers, so the redundancy stays. |
Sent the b-file for the first 20000 terms. The 20000th term is 284727701.
|
[QUOTE=science_man_88;228683][CODE](13:45) gp > ifprime(x,p)=if(isprime(x),p)
%11 = (x,p)->if(isprime(x),p) (13:47) gp > ifprime(100,print(100)) 100[/CODE] I was trying something like this but I can't get it to work. why ? is my best question.[/QUOTE] [QUOTE=axn;228687][CODE]ifprime(p, s)=if(isprime(p), eval(s)); ifprime(100, "print(p)"); ifprime(101, "print(p)"); 101 [/CODE][/QUOTE] [QUOTE=science_man_88;228690]I should of known this by now lol thank you though in a way this is useless but it is useful in the fact that it could create shorter code when typing it in. as by the look of it it saves 10 characters I was also thinking of one where you can do a && part in the if I think I'll try that next.[/QUOTE] Hats off to axn for understanding the request! Another way to do it -- the way I'd prefer -- is to pass a closure rather than a string. So [code]ifprime(p, s)=if(isprime(p), s(p)); ifprime(100, p -> print(p)); ifprime(101, p -> print(p));[/code] or even [code]ifprime(p, s)=if(isprime(p), s(p)); printme(n)=print(n); ifprime(100, printme); ifprime(101, printme);[/code] |
[QUOTE=3.14159;228701]Primes p such that p*10[sup]2[/sup] + 1 is also prime: 701, 1301, 1901, 3701, 6101, 6701, 7901, 10301, 13901, 15101, 16301, 19301, 21101, 22901, 27701, 42101, 46301, 52301, 54101, 60101, 64301, 70901, 72701, 81101, 82301, 87701, 88301, 93701, 102101, 112901, 115301, 117101, 123701, 132701, 138101, 139901, 144701, 148301, 157901, 159701, 165701, 166301, 174101, 178301, 186701, 193301, 198701, 201101, 201701, 217901, 220301, 226901, 237701, 243701, 247301, 250301, 252101, 253901, 270701, 285101, 291701, etc...
Submitted to OEIS.[/QUOTE] "Primes p such that p*10[sup]2[/sup] + 1 is also prime" are 7,13,19,37,61,67,79,103,.... "Primes p such that (p*10)[sup]2[/sup] + 1 is also prime" are 2,11,13,17,23,43,47,89,101,.... I'm not sure what your sequence is! You should probably discuss your sequences here before submitting so we can iron out the wording and so forth. |
[QUOTE=CRGreathouse;228717]"Primes p such that p*10[sup]2[/sup] + 1 is also prime" are 7,13,19,37,61,67,79,103,.... "Primes p such that (p*10)[sup]2[/sup] + 1 is also prime" are 2,11,13,17,23,43,47,89,101,.... I'm not sure what your sequence is!
You should probably discuss your sequences here before submitting so we can iron out the wording and so forth.[/QUOTE] he's using p*10[SUP]n[/SUP]+1 I believe for one thing he's checking and p*10[SUP]100[/SUP] +1 for another I'm unsure about what he's posted as well for once. |
No, he gave the values for p*10^2+1, but as stated "Primes p such..." should list the p-values and not the p*10^2+1 results!
The sequence as submitted is false! |
[QUOTE=kar_bon;228722]No, he gave the values for p*10^2+1, but as stated "Primes p such..." should list the p-values and not the p*10^2+1 results!
The sequence as submitted is false![/QUOTE] so in other terms: [CODE]forprime(p=1,1500,if(isprime(p*10^2+1),print(p*10^2+1)))[/CODE] |
instead of :
[CODE]forprime(p=1,1500,if(isprime(p*10^2+1),print(p)))[/CODE] |
[QUOTE=science_man_88;228723]so in other terms:
[CODE]forprime(p=1,1500,if(isprime(p*10^2+1),print(p*10^2+1)))[/CODE][/QUOTE] NO! "print(p)" and that will give the correct answers! |
Title is corrected.
|
his sequence he gave is A180469 (already up by him) and it's not the p values it's the full values he listed.
|
The title says "Primes such that p * 100 + 1 is also prime."
But to which value is p*100+1 [b]also[/b] prime? The word "Primes" and "p*100+1 also prime" refer to different values! "also prime" to the value p obviously meant, so p has to be listed here, not p*100+1! PS: And the given PARI code is also false [code] p10n(a, x) = { for(n=a, x, if(isprime(n)&isprime(n*10^2+1), print(10^2+1)) )}; [/code] This will ouput only 101 evertime a value is found! Correct [url=http://factordb.com/new/index.php?query=n*100%2B1&use=n&n=1&VP=on&PR=on&perpage=20&format=1&sent=Show]like this[/url]: "Primes p such that p*100+1 also prime" "Values: 7, 13, 19, 37, 61, 67, ..." |
[QUOTE=kar_bon;228730]
PS: And the given PARI code is also false [code] p10n(a, x) = { for(n=a, x, if(isprime(n)&isprime(n*10^2+1), print(10^2+1)) )}; [/code] This will ouput only 101 evertime a value is found![/QUOTE] yeah I know I've been experimenting with it lol |
[QUOTE=kar_bon;228730]
Correct [url=http://factordb.com/new/index.php?query=n*100%2B1&use=n&n=1&VP=on&PR=on&perpage=20&format=1&sent=Show]like this[/url]: "Primes p such that p*100+1 also prime" "Values: 7, 13, 19, 37, 61, 67, ..."[/QUOTE] also I could see using another if to branch if(p<default(primelimit),code1,code2) type idea. |
[QUOTE=science_man_88;228732]also I could see using another if to branch if(p<default(primelimit),code1,code2) type idea.[/QUOTE]
How about this, for your usual preoccupation with looping over primes that 'may' be huge: [code]forallprime(start,stop,ff)={ my(d=default(primelimit)); if (type(ff) != "T_CLOSURE", error("ff must be a function! For example, forallprime(10,100,n->print(n))")); forprime(p=start,min(stop,d), ff(p)); if(stop > d, print("\nWarning, "stop" > primelimit, calculation may be slow"); for(n=precprime(d)+1, stop, if(isprime(n), ff(n)) ) ) }; addhelp(forallprime, "forallprime(start, stop, ff): Executes the closure (function) ff for all the primes between start and stop.");[/code] |
I'm useless as this code appears to do nothing useful no matter what I do if i try to print the result it give me an error oddly if i use next() as well it doesn't give me that error but it only prints once and only exactly what I type even if it's a variable you use, plus the added bonus of not throwing a out of limits error if i go over primelimit.
|
Sorry, typo... replace T_CLOSURE with t_CLOSURE. (That'll teach me to post without checking my programs!)
|
[QUOTE=CRGreathouse;228748]Sorry, typo... replace T_CLOSURE with t_CLOSURE. (That'll teach me to post without checking my programs!)[/QUOTE]
yeah I got it working under that lol. wonder how many primes I could fit on my drives lol using my ipod and a few other removables + main drive I get 211,474,722,816 249,515,175,936 008,002,691,072 [U]006,624,665,600[/U] 475,617,255,424 bytes if i did the math right apparently doing math visually on here isn't a good way although I got nearly accurate I get 372 as the ending by calculator. |
Depends on what size primes you store and how you store them, naturally.
|
[QUOTE=CRGreathouse;228940]Depends on what size primes you store and how you store them, naturally.[/QUOTE]
well if we encode it as characters for each digit it gets inefficient real quick (almost like BCD in ASM I think), I don't know the best way to encode them efficiently unless maybe we convert them to hex to take less room than the decimal and hex or another power of 2 base is easier for me to convert to binary(yes I've learned that part lol). |
[QUOTE=science_man_88;228942]well if we encode it as characters for each digit it gets inefficient real quick (almost like BCD in ASM I think), I don't know the best way to encode them efficiently unless maybe we convert them to hex to take less room than the decimal and hex or another power of 2 base is easier for me to convert to binary(yes I've learned that part lol).[/QUOTE]
So store them natively in binary, taking k bits to store (naturally!) a k-bit prime. You can even do better: store all but the last bit, unless the prime is 2, in which case store it as 0. You'll need some way to delimit the primes, though. |
[QUOTE=CRGreathouse;228946]So store them natively in binary, taking k bits to store (naturally!) a k-bit prime. You can even do better: store all but the last bit, unless the prime is 2, in which case store it as 0.
You'll need some way to delimit the primes, though.[/QUOTE] let me guess the last bit won't matter as it will always be one except when the decimal representation of the prime is 2. yeah what to separate them in storage with. |
[QUOTE=science_man_88;228948]let me guess the last bit won't matter as it will always be one except when the decimal representation of the prime is 2. yeah what to separate them in storage with.[/QUOTE]
One way would be to used fixed width. Store only 65-bit primes, so every 64 bits is a new prime. There are 412246861431389469 65-bit primes, so you shouldn't run out... unless you have more than 206123 TB. You could even left-pad with 0s and include primes of up to 65 bits. This gives you 837903145466607212 to work with (and moves the storage for all of them to 418951 TB). |
[CODE]a=0;forallprime(1,2^100,p->a=a+1;if(ispower(p+1),print(ispower(p+1)","p","a)))[/CODE]
this can't be efficient |
There are many more primes than (proper) powers, and as it happens ispower is a lot slower than isprime. So if you're able, set up the loop the other way: loop over powers and check the the primality of n-1.
This way you're checking something like 1.126 * 10[SUP]15[/SUP] numbers instead of 1.267 * 10[SUP]30[/SUP]. Edit: 1125910730446094 instead of 1.856(1) * 10[SUP]28[/SUP], assuming (wrongly) that you've precalculated primes up to 2[SUP]100[/SUP]. |
[QUOTE=CRGreathouse;228962]There are many more primes than (proper) powers, and as it happens ispower is a lot slower than isprime. So if you're able, set up the loop the other way: loop over powers and check the the primality of n-1.[/QUOTE]
Two words. Algebraic factors. Zero clock cycles. |
[QUOTE=axn;228965]Two words. Algebraic factors. Zero clock cycles.[/QUOTE]
I was actually more interested in the general issue of looping over primes and powers than the specific one of finding primes one less than a power. |
[QUOTE=axn;228965]Two words. Algebraic factors. Zero clock cycles.[/QUOTE]
I can't even remember the think other than 2kp+1 that anyone told me lol. |
The suggestion was that b - 1 is a factor of b^k - 1, so unless b = 2 this cannot be a prime. If b = 2 then you're just looking for Mersenne primes.
|
best I can do of summing up the LL test is:
s=2*y*(2*k*p+1) s=4*y*k*p+(2*y) but this is useless obviously |
[QUOTE=CRGreathouse;226913]
[code]substring(string_to_search, string_to_find)={ my(s1=Vec(string_to_search),s2=Vec(string_to_find),good); for(i=0,#s1-#s2, good=1; for(j=1,#s2, if(s2[j] != s1[j+i], good=0;break)); if(good, return(1)) ); 0 };[/code][/QUOTE] is there a way to create 2 Vec() on the fly to read the parts around the string into and then put out the result of those 2 and the replaced word in between into another that you can return later the only thing with this idea I don't see happening is skipping ahead the length of the new word inserted and making the loop realize to go until the new end of the string. I know I'm annoying lol. |
got it partly in place now I need a way to get the string after s2 is found and I may have something closer to what is needed.
|
I can't understand you, but I'm glad my script seems to have been useful, either as inspiration or as a tool.
|
[QUOTE=CRGreathouse;229045]I can't understand you, but I'm glad my script seems to have been useful, either as inspiration or as a tool.[/QUOTE]
yeah I found a few shortcuts I think that I didn't know before to help me with long winded scripts so I don't have to hold down left/right ctrl +a is one I found helpful. [CODE]substring(string_to_search,string_to_find)=my(s1=Vec(string_to_search),s2=Vec(string_to_find),good);c[COLOR="Red"]=Vec()[/COLOR];for(i=0,#s1-#s2,good=1;for(j=1,#s2,if(s2[j]!=s1[j+i],good=0;c=concat(c,s1[i+1]);if(good,return(c)));0;eturn(c)));0;[/CODE] still not perfect in one sense and I'd want another vector to take from the end of the s2 found in s1 to the end to create a new s1 in one sense with the word/phrase replaced with something understandable. changed c=Vec() to c="" and works better still. |
my computer is weird and i can't find a virus weird anyways just insert a few things into your script for the vector and i can show you how to have fun with it lol.
|
| All times are UTC. The time now is 23:20. |
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.