mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Miscellaneous Math

Reply
 
Thread Tools
Old 2017-02-23, 19:33   #12
a1call
 
a1call's Avatar
 
"Rashid Naimi"
Oct 2015
Remote to Here/There

111100111002 Posts
Default

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:
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
...
a1call is offline   Reply With Quote
Old 2017-02-23, 21:31   #13
rogue
 
rogue's Avatar
 
"Mark"
Apr 2003
Between here and the

2·3,023 Posts
Default

Next up, write sieving code.

Last fiddled with by rogue on 2017-02-23 at 21:32
rogue is online now   Reply With Quote
Old 2017-02-23, 21:44   #14
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

350510 Posts
Default

Quote:
Originally Posted by a1call View Post
Code:
----- snip
      if(isprime(thePRP),print("*** Found a PRP: ",thePRP);/*next(19);*/);
------
Since you are using isprime() you can say "Found a prime:". Otherwise you can use the faster PRP test ispseudoprime()

Last fiddled with by paulunderwood on 2017-02-23 at 21:44
paulunderwood is offline   Reply With Quote
Old 2017-02-24, 00:10   #15
a1call
 
a1call's Avatar
 
"Rashid Naimi"
Oct 2015
Remote to Here/There

22×487 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
Since you are using isprime() you can say "Found a prime:". Otherwise you can use the faster PRP test ispseudoprime()
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.
a1call is offline   Reply With Quote
Old 2017-02-24, 00:19   #16
a1call
 
a1call's Avatar
 
"Rashid Naimi"
Oct 2015
Remote to Here/There

22×487 Posts
Default

Quote:
Originally Posted by rogue View Post
Next up, write sieving code.
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).

Last fiddled with by a1call on 2017-02-24 at 00:19
a1call is offline   Reply With Quote
Old 2017-02-24, 01:05   #17
CRGreathouse
 
CRGreathouse's Avatar
 
Aug 2006

134648 Posts
Default

Quote:
Originally Posted by a1call View Post
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.
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.
CRGreathouse is offline   Reply With Quote
Old 2017-02-24, 03:14   #18
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

2·3·151 Posts
Default

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.
danaj is offline   Reply With Quote
Old 2017-02-24, 03:54   #19
a1call
 
a1call's Avatar
 
"Rashid Naimi"
Oct 2015
Remote to Here/There

36348 Posts
Default

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 ****")
1st result has 1698 dd:
Code:
*** Found a PRP: 528543903400998376216462536295840865232212026591522804126753683198052904246211578867881875721391126838417321626215721768106159386475426442888958677299795608408843964054902579437058466689496242399924972171168925408976405066908253208940881347541387095540660018089923854962992208425473191501327714414505710841422273727416623881022995453335887379692702592594200531837571609053521001389964949653598414785407171401335185108635554440829801569987547006242836426475234211396620101564665587235535216720787505110602826010853255699392206997484574446996175954382170199887094525609056236199660872303080004685319677204135448560789529973804238004318106006061627459069145992636910008770955045784560706446622639548952692716439776027956179796176554336062014549151629649587768017601932205591009197732536815808946288291538533697197894398029898373762645573588574475577190912800589712752371729026379381799279816529232412023570736254890148697309106189572156587639012283349942675392645019677308534553319436083282523713364888912080729845712767211135912421624693545587562006539869213292472659046434591149613470140246109072539285467514012966526552754843809683771804131946428505152427680230371708483012987253337695798858432875334609685222960338694664192661492660645418294435290021649128691712700241796243567677685270314955680987765602454720327598916276133452087212452138828200307651534549792933932088660943061833068984586042449051242604874041565694653252247389351941488811612446683207669960189084223989964241087341462416365963017320532800716633979703511319750931846512819864769257580748344774656713190532034045313285035125145232944695150548276769237354636796949101532300229478794094939177106715673443968320915246886628585226831
n=3 i= 46 Parity=(+)
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

Last fiddled with by a1call on 2017-02-24 at 03:57
a1call is offline   Reply With Quote
Old 2017-02-24, 04:07   #20
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

1101101100012 Posts
Default

Quote:
Originally Posted by a1call View Post
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
number_of_digits = #digits(thePRP) should not use a great deal of stack.

Last fiddled with by paulunderwood on 2017-02-24 at 04:22
paulunderwood is offline   Reply With Quote
Old 2017-02-24, 04:39   #21
a1call
 
a1call's Avatar
 
"Rashid Naimi"
Oct 2015
Remote to Here/There

22×487 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
number_of_digits = #digits(thePRP) should not use a great deal of stack.
Not for this exercise. But I have had stack problems with it in the past. Going by memory probably around 80k dd or less.
a1call is offline   Reply With Quote
Old 2017-02-24, 05:00   #22
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

66618 Posts
Default

Quote:
Originally Posted by a1call View Post
Not for this exercise. But I have had stack problems with it in the past. Going by memory probably around 80k dd or less.
Code:
? allocatemem(100000000)
  ***   Warning: new stack size = 100000000 (95.367 Mbytes).
? n = 10^5000000;#digits(n)
5000001
paulunderwood is offline   Reply With Quote
Reply

Thread Tools


All times are UTC. The time now is 01:27.

Sun Dec 6 01:27:29 UTC 2020 up 2 days, 21:38, 0 users, load averages: 2.44, 2.78, 2.78

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.