mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Miscellaneous Math (https://www.mersenneforum.org/forumdisplay.php?f=56)

 a1call 2017-02-23 19:33

[CODE]print("\nBMT-200-A-Alternative-Factorials=Falsorials-PRPs.gp\n")

allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
/*allocatemem()
allocatemem()
allocatemem()*/

for(n= 3,361,{
forstep (i=3,19^6,1,
falsorial=i;
while(falsorial<10^10000,
thePRP=falsorial-1;
if(isprime(thePRP),print("*** Found a PRP: ",thePRP);/*next(19);*/);
thePRP=falsorial+1;
if(isprime(thePRP),print("*** Found a PRP: ",thePRP);/*next(19);*/);
falsorial=falsorial*(falsorial-1);\\print(falsorial);
);
)
})
print("**** End of Run ****")

[/CODE]

[CODE]BMT-200-A-Alternative-Factorials=Falsorials-PRPs.gp

*** Warning: new stack size = 16000000 (15.259 Mbytes).
*** Warning: new stack size = 32000000 (30.518 Mbytes).
*** Warning: new stack size = 64000000 (61.035 Mbytes).
*** Warning: new stack size = 128000000 (122.070 Mbytes).
*** Warning: new stack size = 256000000 (244.141 Mbytes).
*** Warning: new stack size = 512000000 (488.281 Mbytes).
*** Warning: new stack size = 1024000000 (976.563 Mbytes).
*** Warning: new stack size = 2048000000 (1953.125 Mbytes).
*** Found a PRP: 2
*** Found a PRP: 5
*** Found a PRP: 7
*** Found a PRP: 29
*** Found a PRP: 31
*** Found a PRP: 571580604871
*** Found a PRP: 326704387862983487112031
*** Found a PRP: 16845328252499500651236397153977990330504340553067339752454847112098786861695056576556281766791174572528101695441337548927503867619907346753328587052299845752995063651095010143881791492916271631929764191872816126007112807376635349122779799074047977425098804268779836897066162950536070254479907774146254749119086282986975224594605724404812512447585238430350583157170038516624031
*** Found a PRP: 3
*** Found a PRP: 5
*** Found a PRP: 11
*** Found a PRP: 13
*** Found a PRP: 131
*** Found a PRP: 17291
*** Found a PRP: 17293
*** Found a PRP: 298995971
*** Found a PRP: 19
*** Found a PRP: 379
*** Found a PRP: 430214650034342688019
*** Found a PRP: 34256325853337283636265092964542074845419466233928682710028092071120635978227632019
*** Found a PRP: 5
*** Found a PRP: 7
*** Found a PRP: 29
*** Found a PRP: 31
*** Found a PRP: 571580604871
*** Found a PRP: 326704387862983487112031
...
[/CODE]

 rogue 2017-02-23 21:31

Next up, write sieving code. :smile:

 paulunderwood 2017-02-23 21:44

[QUOTE=a1call;453559][CODE]
----- snip
if(isprime(thePRP),print("*** Found a PRP: ",thePRP);/*next(19);*/);
------
[/CODE]
[/quote]
Since you are using [c]isprime()[/c] you can say "Found a prime:". Otherwise you can use the faster PRP test [c]ispseudoprime()[/c] :smile:

 a1call 2017-02-24 00:10

[QUOTE=paulunderwood;453573]Since you are using [c]isprime()[/c] you can say "Found a prime:". Otherwise you can use the faster PRP test [c]ispseudoprime()[/c] :smile:[/QUOTE]

Thank you Paul. I have never used ispseudoprime before. Somehow I was under the impression that isprime was like mathematica's PrimeQ a PRP test, because of how incredibly fast it works. I just reread the isprime specs and it does indeed seem like a deterministic test if the result is 1 and not 0 or 2.:smile:

 a1call 2017-02-24 00:19

[QUOTE=rogue;453572]Next up, write sieving code. :smile:[/QUOTE]
Hi rogue. The most efficient way I can think of for doing that, would require years for a table or primorial-chunks for gcd testing, complete enough to be used for 1k dd (on a single computer).:smile:

 CRGreathouse 2017-02-24 01:05

[QUOTE=a1call;453589]Thank you Paul. I have never used ispseudoprime before. Somehow I was under the impression that isprime was like mathematica's PrimeQ a PRP test, because of how incredibly fast it works. I just reread the isprime specs and it does indeed seem like a deterministic test if the result is 1 and not 0 or 2.:smile:[/QUOTE]

isprime is indeed a primality test. If no parameter is given (or if the parameter is 0), it uses "a combination of algorithms" that do prove (but do not certify) primality. In practice this means APR-CL unless the number is small enough to prove with Jaeschke/SPRP or Feitsma/BPSW. If you want a certificate you need flag=1.

 danaj 2017-02-24 03:14

Perhaps belongs in a different thread, but there is a Pari branch with ECPP now, from one of Enge's grad students. It isn't super fast yet, but it has amazing promise and being able to use Pari's class polynomial support makes it scale well. A. Enge is looking to integrating some of his work including CM into Pari which would make it even faster. Fast ECPP in Pari, Mmmmmmm good.

 a1call 2017-02-24 03:54

Limiting the upper bound (the while loop) seems to give results more quickly.

[CODE]print("\nBMT-210-G-Alternative-Factorials=Falsorials-PRPs.gp\n")

allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
allocatemem()
/*allocatemem()
allocatemem()
allocatemem()*/

for(n= 3,361,{
forstep (i=3,19^5,1,
falsorial=i;
while(falsorial<10^3000,
if(falsorial>10^1000,
thePRP=falsorial-1;
if(ispseudoprime(thePRP),print("\n*** Found a PRP: ",thePRP);print("n=",n," i= ",i," Parity=(-)");/*next(19);*/);
thePRP=falsorial+1;
if(ispseudoprime(thePRP),print("\n*** Found a PRP: ",thePRP);print("n=",n," i= ",i," Parity=(+)");/*next(19);*/);
);
falsorial=falsorial*(falsorial-1);\\print(falsorial);
);
)
})
print("**** End of Run ****")

[/CODE]1st result has 1698 dd:
[CODE]*** Found a PRP: 528543903400998376216462536295840865232212026591522804126753683198052904246211578867881875721391126838417321626215721768106159386475426442888958677299795608408843964054902579437058466689496242399924972171168925408976405066908253208940881347541387095540660018089923854962992208425473191501327714414505710841422273727416623881022995453335887379692702592594200531837571609053521001389964949653598414785407171401335185108635554440829801569987547006242836426475234211396620101564665587235535216720787505110602826010853255699392206997484574446996175954382170199887094525609056236199660872303080004685319677204135448560789529973804238004318106006061627459069145992636910008770955045784560706446622639548952692716439776027956179796176554336062014549151629649587768017601932205591009197732536815808946288291538533697197894398029898373762645573588574475577190912800589712752371729026379381799279816529232412023570736254890148697309106189572156587639012283349942675392645019677308534553319436083282523713364888912080729845712767211135912421624693545587562006539869213292472659046434591149613470140246109072539285467514012966526552754843809683771804131946428505152427680230371708483012987253337695798858432875334609685222960338694664192661492660645418294435290021649128691712700241796243567677685270314955680987765602454720327598916276133452087212452138828200307651534549792933932088660943061833068984586042449051242604874041565694653252247389351941488811612446683207669960189084223989964241087341462416365963017320532800716633979703511319750931846512819864769257580748344774656713190532034045313285035125145232944695150548276769237354636796949101532300229478794094939177106715673443968320915246886628585226831
n=3 i= 46 Parity=(+)
[/CODE]

ETA BTW, I think the easiest way of counting number of decimal digits is to paste the number in Microsoft Word, and let it do the counting. #digits, takes more memory than I have for semi-large integers

 paulunderwood 2017-02-24 04:07

[QUOTE=a1call;453603]
ETA BTW, I think the easiest way of counting number of decimal digits is to paste the number in Microsoft Word, and let it do the counting. #digits, takes more memory than I have for semi-large integers[/QUOTE]

:no: [c]number_of_digits = #digits(thePRP)[/c] should not use a great deal of stack.

 a1call 2017-02-24 04:39

[QUOTE=paulunderwood;453606]:no: [c]number_of_digits = #digits(thePRP)[/c] should not use a great deal of stack.[/QUOTE]
Not for this exercise. But I have had stack problems with it in the past. Going by memory probably around 80k dd or less.

 paulunderwood 2017-02-24 05:00

[QUOTE=a1call;453607]Not for this exercise. But I have had stack problems with it in the past. Going by memory probably around 80k dd or less.[/QUOTE]

[CODE]? allocatemem(100000000)
*** Warning: new stack size = 100000000 (95.367 Mbytes).
? n = 10^5000000;#digits(n)
5000001
[/CODE]

All times are UTC. The time now is 16:16.