trizen
/
experimental-projects


			
				
					
						
						
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164
							#!/usr/bin/perl

# Author: Daniel "Trizen" Șuteu
# Date: 28 January 2019
# https://github.com/trizen

# A new algorithm for generating Fermat superpseudoprimes to any given base.

# See also:
#   https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html

use 5.020;
use warnings;
use experimental qw(signatures);

use Math::AnyNum qw(prod);
use ntheory qw(:all);
use List::Util qw(min);

sub compute_signature ($base, $p) {

    my $pm1 = $p-1;
    my $valuation = valuation($pm1, 2);
    my $d = $p >> $valuation;

    if (powmod($base, $d, $p) == 1) {
        return 1;
    }

    for my $k (0..1000) {
        if (powmod($base, $d << $k, $p) == $pm1) {
            return -$k;
        }
    }

    die "error for p = $p";
}

sub compute_signature_2($base, $p) {
    valuation(znorder($base, $p), 2);
}

sub fermat_superpseudoprimes ($base, $callback) {

    my %common_divisors;

    #for (my $p = 2; $p <= $prime_limit; $p = next_prime($p)) {

    say "# :: Sieving...";

    my @primes = @{primes(nth_prime(9))};

    #~ my @plist = qw(
        #~ 274177 456037 743731 767131 822571 1072999 3192253 3718651 3835651 4112851 5364991 6693571 6904171 7403131 17575219 18697477 22532959 87876091 549915631 2201771401 3038846581 6077693161 8807085601 33544853431 34517128141 69034256281 91408708819 1223875355821 1843845920281 2447750711641 3687691840561 67280421310721 27778299663977101 55556599327954201 220633985108812507 8163457449026062723 9045993389461312747 142445387161415482404826365418175962266689133006163 5840260873618034778597880982145214452934254453252643 14386984103302963722887462907235772188935602433622363 73938149834061418521192073314311208786743496108043 8355010931248940292894704284517166592902015060208747 12791299921292625404166228683375839120106624826691267
    #~ );

    #~ my @plist = qw(
        #~ 274177 456037 743731 767131 822571 1072999 3192253 3718651 3835651 4112851 5364991 6693571 6904171 7403131 17575219 18697477 22532959 87876091 549915631 2201771401 3038846581 6077693161 8807085601 33544853431 34517128141 69034256281 91408708819 1223875355821 1843845920281 2447750711641 3687691840561 67280421310721 27778299663977101 55556599327954201 220633985108812507 8163457449026062723 9045993389461312747 73938149834061418521192073314311208786743496108043 142445387161415482404826365418175962266689133006163 5840260873618034778597880982145214452934254453252643 8355010931248940292894704284517166592902015060208747 12791299921292625404166228683375839120106624826691267 14386984103302963722887462907235772188935602433622363 123282949929736752510916282638002560328626287433241467864741378859343760091491850110380631340085554443 184423640821245254099510263989726185275555264354779963459730437466872086414580403681648321609499543203 20839871412800713713244659830839058936137744872090135870949539433756545764847585616026260341873448381827 28724927333628663335043493854654596556569924971945262012484741274227096101317601075718687102239934184987 29711190933066557355130824115758617039198935271411193755402672305101846182049535876601732152960618620523 31905289862075428959215275670222630052671060733376933678533365681768870949722409836925159638443420973947 29674495668685510550154174642905332730771991799853043350995075531276838753171770199594238596428121188033664754218345562493168782883 9288117144298564802198256663229369144731633433354002568861458641289650529742764072472996680682001931854537068070342161060361829042067 10475096971045985224204423648945582453962513105348124302901261662540724079869634880456766224539126779375883658239075983560088580357347 222116839510593004474253307523543930668605023776267914871251681339939729863340244743625878412836415901800725569032731901731279834484776044761998664069644323493578946903234375700665946172079826882759862218316218817048571800819527667 142376894126290115867996370122591659558575820240587733432472327738901366842401096880664188062628142593054265089749981149009750373904741444692441143668642011359384104964973234824126871496303169031849071681940696261728134524325317233907 150373100348671464029069489193439241062645601096533378367837388267139197117481345691434719685490253565519091210235159497472076447946193382303873095575149207005152947053489672349350845558498042799628426721800080139141883109154820229883 14944654354308364545703605962105292139354124383675054870648192523437530482472401251725664068684814949588779644083289776272581318830344624110152295134965367218621910337465901747035296579855521540699486195373540899785897644915651031264255203023455677284343404880972742346387613171333560275622453586156561794632202184300886820826624562851321402814660588726418255372307360994387346859948919643309365459081627814846050589608870993523 1900809651964297932845396573668937798734368301116723662488816465653850457169635536238329413823488311896159501537950387704882660206924124702319767792112295081202063239825692483258274851577725555804284352026976128284526935861630031791154642239525653834047162490008249561015133998516792689132376286139252916898534542631569361562095241253775712906006443172754937344955794702245521682959506315285400771587080099118585628535260790657982364170344198148111227942260024011356635573591350816373257142371334822961251413464688009750724419940054051790800253090303365833973258551947851236063024914180757185378753321451274470828789328966010596890419576693866086998173674081322567916936236887943395514342216042478259443 5004831813621996457181929178470313224067591736840333403333053754066588253727650366915521346597244725222587967549423370826956044324831220341207948596631672948805032510461048308419037684204151388432680698887028145773159422123671873706110173016671046545046178836191721094152847818094715150485546761404652930193841450748922128992996770221191452081514964873863750029268607451012458591232380128146460231588781900979235959933341661802467564860516273723976863171970643221902021465266026699510786055863724588856974971652523529673657397702162318365177066386768762240851589767278692304553944599037933669102257495381205681692202303167505901612474745434949407066191283856122321325293111725954960389263054839845257110787 11178661563202036143063777249747023194356819978867451859096729634510294538614626588617615282695934762261314028544686230092414924676920777374342554385412407372549333913414897494041914402128603993684996274270646610441302909802246216963780451010650370198031362603738515668330003045277257804787504938784946404280281645216259415346682113813454967600223892298971786525685028643905913017484856640193441937703618062916402081415868709859594283685794229309042131528431201210788373808290734151091125254285820093835119562585830185344010313667457878581696288424074094469596733544005313119286649520297033007212448283454945162944110043649108320312557530536626457636259377272258021919502009137995109019846572545814643778403
    #~ );

    #~ my @plist = qw(
        #~ 274177 743731 767131 822571 1072999 3718651 3835651 4112851 5364991 6693571 6904171 7403131 17575219 22532959 87876091 549915631 3038846581 4540612081 6077693161 9081224161 33544853431 91408708819 399165290221 424665351661 798330580441 849330703321 1287836182261 2575672364521 22754930352733 54786377365501 67280421310721 68264791058197 109572754731001 172157429516701 344314859033401 531099297693901 1062198595387801 27778299663977101 55556599327954201 123282949929736752510916282638002560328626287433241467864741378859343760091491850110380631340085554443 184423640821245254099510263989726185275555264354779963459730437466872086414580403681648321609499543203 28724927333628663335043493854654596556569924971945262012484741274227096101317601075718687102239934184987 29711190933066557355130824115758617039198935271411193755402672305101846182049535876601732152960618620523 29674495668685510550154174642905332730771991799853043350995075531276838753171770199594238596428121188033664754218345562493168782883 9288117144298564802198256663229369144731633433354002568861458641289650529742764072472996680682001931854537068070342161060361829042067 10475096971045985224204423648945582453962513105348124302901261662540724079869634880456766224539126779375883658239075983560088580357347 222116839510593004474253307523543930668605023776267914871251681339939729863340244743625878412836415901800725569032731901731279834484776044761998664069644323493578946903234375700665946172079826882759862218316218817048571800819527667 14944654354308364545703605962105292139354124383675054870648192523437530482472401251725664068684814949588779644083289776272581318830344624110152295134965367218621910337465901747035296579855521540699486195373540899785897644915651031264255203023455677284343404880972742346387613171333560275622453586156561794632202184300886820826624562851321402814660588726418255372307360994387346859948919643309365459081627814846050589608870993523 15138934860914373284797752839612660937165728000662830583966619026242218378744542467998097701577717543933433779456372543364124875975139104223584274971719916992463995171852958469746755435393643320728579515913396931483114314299554494670690520662760601089039869144425387996890652142560896559205545482776597097962420812696798349497370682168388581051251176379861692692147356687314382369128255598672387210049688976439049247273786316437787 30681375389395072412329503040202164762094017359684887649440739250617250080515839769792788333009925091505764609302993910687609447558697513298142661912083898899830781922817496286663463878443385723056045159101879467260447865011831567185515931807154505464757010220637040037133769840747799245852897212379421364379911084369720643157060227533762839978498188655336678279347012121477223103475132027714127287494581903878941860467012149700667 1900809651964297932845396573668937798734368301116723662488816465653850457169635536238329413823488311896159501537950387704882660206924124702319767792112295081202063239825692483258274851577725555804284352026976128284526935861630031791154642239525653834047162490008249561015133998516792689132376286139252916898534542631569361562095241253775712906006443172754937344955794702245521682959506315285400771587080099118585628535260790657982364170344198148111227942260024011356635573591350816373257142371334822961251413464688009750724419940054051790800253090303365833973258551947851236063024914180757185378753321451274470828789328966010596890419576693866086998173674081322567916936236887943395514342216042478259443 5004831813621996457181929178470313224067591736840333403333053754066588253727650366915521346597244725222587967549423370826956044324831220341207948596631672948805032510461048308419037684204151388432680698887028145773159422123671873706110173016671046545046178836191721094152847818094715150485546761404652930193841450748922128992996770221191452081514964873863750029268607451012458591232380128146460231588781900979235959933341661802467564860516273723976863171970643221902021465266026699510786055863724588856974971652523529673657397702162318365177066386768762240851589767278692304553944599037933669102257495381205681692202303167505901612474745434949407066191283856122321325293111725954960389263054839845257110787 11178661563202036143063777249747023194356819978867451859096729634510294538614626588617615282695934762261314028544686230092414924676920777374342554385412407372549333913414897494041914402128603993684996274270646610441302909802246216963780451010650370198031362603738515668330003045277257804787504938784946404280281645216259415346682113813454967600223892298971786525685028643905913017484856640193441937703618062916402081415868709859594283685794229309042131528431201210788373808290734151091125254285820093835119562585830185344010313667457878581696288424074094469596733544005313119286649520297033007212448283454945162944110043649108320312557530536626457636259377272258021919502009137995109019846572545814643778403
    #~ );

    #~ my @plist = qw(
        #~ 3 5 7 11 13 17 19 29 31 37 53 61 67 73 97 101 107 113 127 137 181 193 211 229 241 257 281 313 331 379 409 421 457 461 521 523 617 673 701 761 769 883 991 1009 1171 1289 1471 1489 1801 2011 2861 4481 5101 5209 5237 5521 6203 11161 11971 12689 22051 22543 92569 157321 162289 494761
        #~ 211 457 2857 3169 23371 49257589 98515177 1468383841 10688896597 123282949929736752510916282638002560328626287433241467864741378859343760091491850110380631340085554443 28724927333628663335043493854654596556569924971945262012484741274227096101317601075718687102239934184987 29711190933066557355130824115758617039198935271411193755402672305101846182049535876601732152960618620523 29674495668685510550154174642905332730771991799853043350995075531276838753171770199594238596428121188033664754218345562493168782883 9288117144298564802198256663229369144731633433354002568861458641289650529742764072472996680682001931854537068070342161060361829042067 10475096971045985224204423648945582453962513105348124302901261662540724079869634880456766224539126779375883658239075983560088580357347 222116839510593004474253307523543930668605023776267914871251681339939729863340244743625878412836415901800725569032731901731279834484776044761998664069644323493578946903234375700665946172079826882759862218316218817048571800819527667 142376894126290115867996370122591659558575820240587733432472327738901366842401096880664188062628142593054265089749981149009750373904741444692441143668642011359384104964973234824126871496303169031849071681940696261728134524325317233907 150373100348671464029069489193439241062645601096533378367837388267139197117481345691434719685490253565519091210235159497472076447946193382303873095575149207005152947053489672349350845558498042799628426721800080139141883109154820229883 14944654354308364545703605962105292139354124383675054870648192523437530482472401251725664068684814949588779644083289776272581318830344624110152295134965367218621910337465901747035296579855521540699486195373540899785897644915651031264255203023455677284343404880972742346387613171333560275622453586156561794632202184300886820826624562851321402814660588726418255372307360994387346859948919643309365459081627814846050589608870993523 15138934860914373284797752839612660937165728000662830583966619026242218378744542467998097701577717543933433779456372543364124875975139104223584274971719916992463995171852958469746755435393643320728579515913396931483114314299554494670690520662760601089039869144425387996890652142560896559205545482776597097962420812696798349497370682168388581051251176379861692692147356687314382369128255598672387210049688976439049247273786316437787 30681375389395072412329503040202164762094017359684887649440739250617250080515839769792788333009925091505764609302993910687609447558697513298142661912083898899830781922817496286663463878443385723056045159101879467260447865011831567185515931807154505464757010220637040037133769840747799245852897212379421364379911084369720643157060227533762839978498188655336678279347012121477223103475132027714127287494581903878941860467012149700667 1900809651964297932845396573668937798734368301116723662488816465653850457169635536238329413823488311896159501537950387704882660206924124702319767792112295081202063239825692483258274851577725555804284352026976128284526935861630031791154642239525653834047162490008249561015133998516792689132376286139252916898534542631569361562095241253775712906006443172754937344955794702245521682959506315285400771587080099118585628535260790657982364170344198148111227942260024011356635573591350816373257142371334822961251413464688009750724419940054051790800253090303365833973258551947851236063024914180757185378753321451274470828789328966010596890419576693866086998173674081322567916936236887943395514342216042478259443 5004831813621996457181929178470313224067591736840333403333053754066588253727650366915521346597244725222587967549423370826956044324831220341207948596631672948805032510461048308419037684204151388432680698887028145773159422123671873706110173016671046545046178836191721094152847818094715150485546761404652930193841450748922128992996770221191452081514964873863750029268607451012458591232380128146460231588781900979235959933341661802467564860516273723976863171970643221902021465266026699510786055863724588856974971652523529673657397702162318365177066386768762240851589767278692304553944599037933669102257495381205681692202303167505901612474745434949407066191283856122321325293111725954960389263054839845257110787 11178661563202036143063777249747023194356819978867451859096729634510294538614626588617615282695934762261314028544686230092414924676920777374342554385412407372549333913414897494041914402128603993684996274270646610441302909802246216963780451010650370198031362603738515668330003045277257804787504938784946404280281645216259415346682113813454967600223892298971786525685028643905913017484856640193441937703618062916402081415868709859594283685794229309042131528431201210788373808290734151091125254285820093835119562585830185344010313667457878581696288424074094469596733544005313119286649520297033007212448283454945162944110043649108320312557530536626457636259377272258021919502009137995109019846572545814643778403
    #~ );

    #~ my @plist = qw(
        #~ 274177 2147569201 3037001161 3037002877 3037217041 3037419001 3038555521 3039036001 3040228801 3041910481 3047259601 3050261761 3050679241 3082383361 3170665741 3347875621 4044647521 6074002321 6074005753 6074434081 6074838001 6077111041 6078072001 6080457601 6083820961 6094519201 6100523521 6101358481 6164766721 6341331481 6695751241 8089295041 8590276801 52456994521 102835822741 104913989041 205671645481 1515590649721 3031181299441 20432432259121 40864864518241 67280421310721 15875058966304285441 31750117932608570881 866343646137052296813601 1732687292274104593627201 896052973183478531313669116341 1792105946366957062627338232681 118042515934659926299185815952813181 236085031869319852598371631905626361 46111359278910368495062042946635121401 92222718557820736990124085893270242801 594836534697943753586300354011593066061 1189673069395887507172600708023186132121 105719513418757801848628745863750342833781 112239659620795727953830518736404548999741 211439026837515603697257491727500685667561 224479319241591455907661037472809097999481 5310489179590158108408394656979335298870147414321 10620978359180316216816789313958670597740294828641 13149782730413724839868405817282163597202269787841 26299565460827449679736811634564327194404539575681 1356893205494566231913921125270803256186309213732741 2713786410989132463827842250541606512372618427465481 921347491855125486210849448111794485277702788094439844281 1842694983710250972421698896223588970555405576188879688561 94183565626585432123476118743718656564586398285465035430418381 188367131253170864246952237487437313129172796570930070860836761 2571703890109274439924538995759915900728161270960305578655806689261 5143407780218548879849077991519831801456322541920611157311613378521 344735176623482245886386745292606125652318048757024280998914929874661 689470353246964491772773490585212251304636097514048561997829859749321 39315217293808751773881000102811516475231400264390467443082957421571561 56640567287690574589489576419304727125333373262257453095966972556501401 78630434587617503547762000205623032950462800528780934886165914843143121 113281134575381149178979152838609454250666746524514906191933945113002801 7454898194482803567293113366952607467231377804370650076601535358539522501 14909796388965607134586226733905214934462755608741300153203070717079045001 1867481817310720224821194070875095980971329065506065326788991013455584544341 3734963634621440449642388141750191961942658131012130653577982026911169088681 183498979262803081512921838212750038037755044997327193584312475362609720631421 366997958525606163025843676425500076075510089994654387168624950725219441262841 10743968325151313008582423583445398294547657746978337143459014709152351170902581 21487936650302626017164847166890796589095315493956674286918029418304702341805161 2351735288949787402989708273265270515584320640171924908068250997447792422964231401 4703470577899574805979416546530541031168641280343849816136501994895584845928462801 227118130009739033507510179578111520461401090528411596419218919603221318910043758185581 454236260019478067015020359156223040922802181056823192838437839206442637820087516371161 691229091333988362848944024802948105752090275521252684754144537922847492334915785782201 1382458182667976725697888049605896211504180551042505369508289075845694984669831571564401
    #~ );

    my @plist = qw(
        23 89 151 751 829 1303 1657 2251 6763 10627 11251 16927 28351 29947 149491 157543 747451 10670053 32010157 34233211
        274177 424771 743731 767131 2123851 3718651 3835651 6693571 6904171 20813731 549915631 4540612081 9081224161 33544853431 399165290221 798330580441 1287836182261 2575672364521 54786377365501 67280421310721 109572754731001 172157429516701 344314859033401 531099297693901 1062198595387801 27778299663977101 55556599327954201 123282949929736752510916282638002560328626287433241467864741378859343760091491850110380631340085554443 28724927333628663335043493854654596556569924971945262012484741274227096101317601075718687102239934184987 29711190933066557355130824115758617039198935271411193755402672305101846182049535876601732152960618620523 14944654354308364545703605962105292139354124383675054870648192523437530482472401251725664068684814949588779644083289776272581318830344624110152295134965367218621910337465901747035296579855521540699486195373540899785897644915651031264255203023455677284343404880972742346387613171333560275622453586156561794632202184300886820826624562851321402814660588726418255372307360994387346859948919643309365459081627814846050589608870993523 15138934860914373284797752839612660937165728000662830583966619026242218378744542467998097701577717543933433779456372543364124875975139104223584274971719916992463995171852958469746755435393643320728579515913396931483114314299554494670690520662760601089039869144425387996890652142560896559205545482776597097962420812696798349497370682168388581051251176379861692692147356687314382369128255598672387210049688976439049247273786316437787 30681375389395072412329503040202164762094017359684887649440739250617250080515839769792788333009925091505764609302993910687609447558697513298142661912083898899830781922817496286663463878443385723056045159101879467260447865011831567185515931807154505464757010220637040037133769840747799245852897212379421364379911084369720643157060227533762839978498188655336678279347012121477223103475132027714127287494581903878941860467012149700667 1900809651964297932845396573668937798734368301116723662488816465653850457169635536238329413823488311896159501537950387704882660206924124702319767792112295081202063239825692483258274851577725555804284352026976128284526935861630031791154642239525653834047162490008249561015133998516792689132376286139252916898534542631569361562095241253775712906006443172754937344955794702245521682959506315285400771587080099118585628535260790657982364170344198148111227942260024011356635573591350816373257142371334822961251413464688009750724419940054051790800253090303365833973258551947851236063024914180757185378753321451274470828789328966010596890419576693866086998173674081322567916936236887943395514342216042478259443
    );

    my %seen_p;
    #while (<>) {
    #    my $p = (split(' ', $_))[-1];
    foreach my $p(@plist) {

        $p || next;
        $p =~ /^[0-9]+\z/ or next;
        $p > $primes[-1] or next;
        is_prime($p) || next;
        next if $seen_p{$p}++;

        if ($p > ~0) {
            $p = Math::GMPz->new($p);
        }

        #my $z = znorder($base, $p) // next;
        #my @p_orders = map { znorder($_, $p) } @primes;
        #my $p_sig = join(' ', map { valuation($_, 2) } @p_orders);

        my $p_sig = join(' ', map { compute_signature($_, $p) } @primes);
        #my $p_sig = join(' ', map { compute_signature_2($_, $p) } @primes);

        #$z < 1e10 or next;

        my $pm1 = $p-1;
        #my $pm1 = lcm(@p_orders);
        my $p_added = 0;

        for my $k (2..6) {
            my @list = (addint(mulint($pm1, $k), 1));

            if ($pm1 % $k == 0) {
                push @list, addint(divint($pm1, $k), 1);
            }

            foreach my $q (@list) {

                $q > $primes[-1] or next;

                if (is_prime($q) and $p != $q) {
                #my @q_orders = map { znorder($_, $q) } @primes;
                #my $q_sig = join(' ', map {valuation($_, 2)} @q_orders);
                my $q_sig = join(' ', map { compute_signature($_, $q) } @primes);
                #my $q_sig = join(' ', map { compute_signature_2($_, $q) } @primes);

                if ($p_sig eq $q_sig) {
                    if (not $p_added) {
                        push @{$common_divisors{$p_sig}}, $p;
                        $p_added = 1;
                    }
                    push @{$common_divisors{$p_sig}}, $q;
                }
            }
        }
        }
    }

    say "# :: Creating combinations...";

    my %seen;

    foreach my $arr (values %common_divisors) {

        my $l = scalar(@$arr);

        foreach my $k (2 .. min($l, 3)) {
        #foreach my $k (2 .. $l) {
            forcomb {
                my $n = Math::Prime::Util::GMP::vecprod(@{$arr}[@_]);
                $callback->($n) if !$seen{$n}++;
            } $l, $k;
        }
    }
}

my $base = 2;

fermat_superpseudoprimes(
    $base,           # base
    sub ($n) {
        if ($n > ~0 and  Math::Prime::Util::GMP::is_pseudoprime($n, $base)) {
            say $n;
        }
    }
);