123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556 |
- #!/usr/bin/perl
- =begin
- This script factorizes a natural number given as a command line
- parameter into its prime factors. It first attempts to use trial
- division to find very small factors, then uses other special-purpose
- factorization methods to find slightly larger factors. If any large
- factors remain, it uses the Self-Initializing Quadratic Sieve (SIQS) [2]
- to factorize those.
- [2] Contini, Scott Patrick. 'Factoring integers with the self-
- initializing quadratic sieve.' (1997).
- =cut
- use 5.020;
- use strict;
- use warnings;
- use Math::GMPz;
- use POSIX qw(ULONG_MAX);
- use experimental qw(signatures);
- use Math::Sidef qw();
- use Math::AnyNum qw(is_smooth);
- $Sidef::Types::Number::Number::USE_YAFU = 1;
- $Sidef::Types::Number::Number::USE_FACTORDB = 0;
- $Sidef::Types::Number::Number::VERBOSE = 1;
- $Sidef::Types::Number::Number::SPECIAL_FACTORS = 1;
- use ntheory qw(
- urandomm valuation sqrtmod invmod random_prime factor_exp vecmin vecall
- );
- use Math::Prime::Util::GMP qw(
- is_power powmod vecprod sqrtint rootint logint is_prime trial_factor
- gcd sieve_primes consecutive_integer_lcm lucas_sequence is_pseudoprime
- );
- my $ZERO = Math::GMPz->new(0);
- my $ONE = Math::GMPz->new(1);
- local $| = 1;
- # Tuning parameters
- use constant {
- FACTOR_MIN_LEN => 25, # factor directly numbers with <= this many digits
- PREFACTORIZATION => 0, # try to find very small prime factors of n
- FACTORIZE_WITH_SIDEF => 1, # true to factorize large numbers with Math::Sidef
- MASK_LIMIT => 200, # show Cn if n > MASK_LIMIT, where n ~ log_10(N)
- LOOK_FOR_SMALL_FACTORS => 1,
- LOOK_FOR_SPECIAL_FORMS => 0,
- SCRAPE_FACTORDB => 0, # factorize with factorDB by scrapping instead of API
- MAX_SIZE => 1e4, # ignore numbers with length() greater than this
- TRIAL_DIVISION_LIMIT => 1_000_000,
- PHI_FINDER_ITERATIONS => 100_000,
- FERMAT_ITERATIONS => 100_000,
- NEAR_POWER_ITERATIONS => 1_000,
- CFRAC_ITERATIONS => 50_000,
- HOLF_ITERATIONS => 100_000,
- MILLER_RABIN_ITERATIONS => 100,
- LUCAS_MILLER_ITERATIONS => 50,
- SIQS_TRIAL_DIVISION_EPS => 25,
- SIQS_MIN_PRIME_POLYNOMIAL => 400,
- SIQS_MAX_PRIME_POLYNOMIAL => 4000,
- };
- my @small_primes = sieve_primes(2, TRIAL_DIVISION_LIMIT);
- package Polynomial {
- sub new ($class, $coeff, $A = undef, $B = undef) {
- bless {
- a => $A,
- b => $B,
- coeff => $coeff,
- }, $class;
- }
- sub eval ($self, $x) {
- my $res = $ZERO;
- foreach my $k (@{$self->{coeff}}) {
- $res *= $x;
- $res += $k;
- }
- return $res;
- }
- }
- package FactorBasePrime {
- sub new ($class, $p, $t, $lp) {
- bless {
- p => $p,
- soln1 => undef,
- soln2 => undef,
- t => $t,
- lp => $lp,
- ainv => undef,
- }, $class;
- }
- }
- sub siqs_factor_base_primes ($n, $nf) {
- my @factor_base;
- foreach my $p (@small_primes) {
- my $t = sqrtmod($n, $p) // next;
- my $lp = sprintf('%0.f', log($p) / log(2));
- push @factor_base, FactorBasePrime->new($p, $t, $lp);
- if (scalar(@factor_base) >= $nf) {
- last;
- }
- }
- return \@factor_base;
- }
- sub siqs_create_poly ($A, $B, $n, $factor_base, $first) {
- my $B_orig = $B;
- if (($B << 1) > $A) {
- $B = $A - $B;
- }
- # 0 < $B or die 'error';
- # 2 * $B <= $A or die 'error';
- # ($B * $B - $n) % $A == 0 or die 'error';
- my $g = Polynomial->new([$A * $A, ($A * $B) << 1, $B * $B - $n], $A, $B_orig);
- my $h = Polynomial->new([$A, $B]);
- foreach my $fb (@$factor_base) {
- next if Math::GMPz::Rmpz_divisible_ui_p($A, $fb->{p});
- #<<<
- $fb->{ainv} = int(invmod($A, $fb->{p})) if $first;
- $fb->{soln1} = int(($fb->{ainv} * ( $fb->{t} - $B)) % $fb->{p});
- $fb->{soln2} = int(($fb->{ainv} * (-$fb->{t} - $B)) % $fb->{p});
- #>>>
- }
- return ($g, $h);
- }
- sub siqs_find_first_poly ($n, $m, $factor_base) {
- my $p_min_i;
- my $p_max_i;
- foreach my $i (0 .. $#{$factor_base}) {
- my $fb = $factor_base->[$i];
- if (not defined($p_min_i) and $fb->{p} >= SIQS_MIN_PRIME_POLYNOMIAL) {
- $p_min_i = $i;
- }
- if (not defined($p_max_i) and $fb->{p} > SIQS_MAX_PRIME_POLYNOMIAL) {
- $p_max_i = $i - 1;
- last;
- }
- }
- # The following may happen if the factor base is small
- if (not defined($p_max_i)) {
- $p_max_i = $#{$factor_base};
- }
- if (not defined($p_min_i)) {
- $p_min_i = 5;
- }
- if ($p_max_i - $p_min_i < 20) {
- $p_min_i = vecmin($p_min_i, 5);
- }
- my $target0 = (log("$n") + log(2)) / 2 - log("$m");
- my $target1 = $target0 - log(($factor_base->[$p_min_i]{p} + $factor_base->[$p_max_i]{p}) / 2) / 2;
- # find q such that the product of factor_base[q_i] is approximately
- # sqrt(2 * n) / m; try a few different sets to find a good one
- my ($best_q, $best_a, $best_ratio);
- for (1 .. 30) {
- my $A = $ONE;
- my $log_A = 0;
- my %Q;
- while ($log_A < $target1) {
- my $p_i = 0;
- while ($p_i == 0 or exists $Q{$p_i}) {
- $p_i = $p_min_i + urandomm($p_max_i - $p_min_i + 1);
- }
- my $fb = $factor_base->[$p_i];
- $A *= $fb->{p};
- $log_A += log($fb->{p});
- $Q{$p_i} = $fb;
- }
- my $ratio = exp($log_A - $target0);
- # ratio too small seems to be not good
- if ( !defined($best_ratio)
- or ($ratio >= 0.9 and $ratio < $best_ratio)
- or ($best_ratio < 0.9 and $ratio > $best_ratio)) {
- $best_q = \%Q;
- $best_a = $A;
- $best_ratio = $ratio;
- }
- }
- my $A = $best_a;
- my $B = $ZERO;
- my @arr;
- foreach my $fb (values %$best_q) {
- my $p = $fb->{p};
- #($A % $p == 0) or die 'error';
- my $r = $A / $p;
- #$fb->{t} // die 'error';
- #gcd($r, $p) == 1 or die 'error';
- my $gamma = ($fb->{t} * int(invmod($r, $p))) % $p;
- if ($gamma > ($p >> 1)) {
- $gamma = $p - $gamma;
- }
- my $t = $r * $gamma;
- $B += $t;
- push @arr, $t;
- }
- my ($g, $h) = siqs_create_poly($A, $B, $n, $factor_base, 1);
- return ($g, $h, \@arr);
- }
- sub siqs_find_next_poly ($n, $factor_base, $i, $g, $arr) {
- # Compute the (i+1)-th polynomials for the Self-Initializing
- # Quadratic Sieve, given that g is the i-th polynomial.
- my $v = valuation($i, 2);
- my $z = ((($i >> ($v + 1)) & 1) == 0) ? -1 : 1;
- my $A = $g->{a};
- my $B = ($g->{b} + 2 * $z * $arr->[$v]) % $A;
- return siqs_create_poly($A, $B, $n, $factor_base, 0);
- }
- sub siqs_sieve ($factor_base, $m) {
- # Perform the sieving step of the SIQS. Return the sieve array.
- my @sieve_array = (0) x (2 * $m + 1);
- foreach my $fb (@$factor_base) {
- $fb->{p} > 100 or next;
- $fb->{soln1} // next;
- my $p = $fb->{p};
- my $lp = $fb->{lp};
- my $end = 2 * $m;
- my $i_start_1 = -int(($m + $fb->{soln1}) / $p);
- my $a_start_1 = int($fb->{soln1} + $i_start_1 * $p);
- for (my $i = $a_start_1 + $m ; $i <= $end ; $i += $p) {
- $sieve_array[$i] += $lp;
- }
- my $i_start_2 = -int(($m + $fb->{soln2}) / $p);
- my $a_start_2 = int($fb->{soln2} + $i_start_2 * $p);
- for (my $i = $a_start_2 + $m ; $i <= $end ; $i += $p) {
- $sieve_array[$i] += $lp;
- }
- }
- return \@sieve_array;
- }
- sub siqs_trial_divide ($n, $factor_base_info) {
- # Determine whether the given number can be fully factorized into
- # primes from the factors base. If so, return the indices of the
- # factors from the factor base. If not, return undef.
- my $factor_prod = $factor_base_info->{prod};
- state $g = Math::GMPz::Rmpz_init_nobless();
- state $t = Math::GMPz::Rmpz_init_nobless();
- Math::GMPz::Rmpz_set($t, $n);
- Math::GMPz::Rmpz_gcd($g, $t, $factor_prod);
- while (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
- Math::GMPz::Rmpz_remove($t, $t, $g);
- if (Math::GMPz::Rmpz_cmp_ui($t, 1) == 0) {
- my $factor_index = $factor_base_info->{index};
- return [map { [$factor_index->{$_->[0]}, $_->[1]] } factor_exp($n)];
- }
- Math::GMPz::Rmpz_gcd($g, $t, $g);
- }
- return undef;
- }
- sub siqs_trial_division ($n, $sieve_array, $factor_base_info, $smooth_relations, $g, $h, $m, $req_relations) {
- # Perform the trial division step of the Self-Initializing Quadratic Sieve.
- my $limit = (log("$m") + log("$n") / 2) / log(2) - SIQS_TRIAL_DIVISION_EPS;
- foreach my $i (0 .. $#{$sieve_array}) {
- next if ((my $sa = $sieve_array->[$i]) < $limit);
- my $x = $i - $m;
- my $gx = abs($g->eval($x));
- my $divisors_idx = siqs_trial_divide($gx, $factor_base_info) // next;
- my $u = $h->eval($x);
- my $v = $gx;
- #(($u * $u) % $n == ($v % $n)) or die 'error';
- push @$smooth_relations, [$u, $v, $divisors_idx];
- if (scalar(@$smooth_relations) >= $req_relations) {
- return 1;
- }
- }
- return 0;
- }
- sub siqs_build_matrix ($factor_base, $smooth_relations) {
- # Build the matrix for the linear algebra step of the Quadratic Sieve.
- my $fb = scalar(@$factor_base);
- my @matrix;
- foreach my $sr (@$smooth_relations) {
- my @row = (0) x $fb;
- foreach my $pair (@{$sr->[2]}) {
- $row[$pair->[0]] = $pair->[1] % 2;
- }
- push @matrix, \@row;
- }
- return \@matrix;
- }
- sub siqs_build_matrix_opt ($M) {
- # Convert the given matrix M of 0s and 1s into a list of numbers m
- # that correspond to the columns of the matrix.
- # The j-th number encodes the j-th column of matrix M in binary:
- # The i-th bit of m[i] is equal to M[i][j].
- my $m = scalar(@{$M->[0]});
- my @cols_binary = ("") x $m;
- foreach my $mi (@$M) {
- foreach my $j (0 .. $#{$mi}) {
- $cols_binary[$j] .= $mi->[$j];
- }
- }
- #<<<
- return ([map {
- Math::GMPz::Rmpz_init_set_str(scalar reverse($_), 2)
- } @cols_binary], scalar(@$M), $m);
- #>>>
- }
- sub find_pivot_column_opt ($M, $j) {
- # For a matrix produced by siqs_build_matrix_opt, return the row of
- # the first non-zero entry in column j, or None if no such row exists.
- my $v = $M->[$j];
- if ($v == 0) {
- return undef;
- }
- return valuation($v, 2);
- }
- sub siqs_solve_matrix_opt ($M, $n, $m) {
- # Perform the linear algebra step of the SIQS. Perform fast
- # Gaussian elimination to determine pairs of perfect squares mod n.
- # Use the optimizations described in [1].
- # [1] Koç, Çetin K., and Sarath N. Arachchige. 'A Fast Algorithm for
- # Gaussian Elimination over GF (2) and its Implementation on the
- # GAPP.' Journal of Parallel and Distributed Computing 13.1
- # (1991): 118-122.
- my @row_is_marked = (0) x $n;
- my @pivots = (-1) x $m;
- foreach my $j (0 .. $m - 1) {
- my $i = find_pivot_column_opt($M, $j) // next;
- $pivots[$j] = $i;
- $row_is_marked[$i] = 1;
- foreach my $k (0 .. $m - 1) {
- if ($k != $j and Math::GMPz::Rmpz_tstbit($M->[$k], $i)) {
- Math::GMPz::Rmpz_xor($M->[$k], $M->[$k], $M->[$j]);
- }
- }
- }
- my @perf_squares;
- foreach my $i (0 .. $n - 1) {
- if (not $row_is_marked[$i]) {
- my @perfect_sq_indices = ($i);
- foreach my $j (0 .. $m - 1) {
- if (Math::GMPz::Rmpz_tstbit($M->[$j], $i)) {
- push @perfect_sq_indices, $pivots[$j];
- }
- }
- push @perf_squares, \@perfect_sq_indices;
- }
- }
- return \@perf_squares;
- }
- sub siqs_calc_sqrts ($n, $square_indices, $smooth_relations) {
- # Given on of the solutions returned by siqs_solve_matrix_opt and
- # the corresponding smooth relations, calculate the pair [a, b], such
- # that a^2 = b^2 (mod n).
- my $r1 = $ONE;
- my $r2 = $ONE;
- foreach my $i (@$square_indices) {
- ($r1 *= $smooth_relations->[$i][0]) %= $n;
- ($r2 *= $smooth_relations->[$i][1]);
- }
- $r2 = Math::GMPz->new(sqrtint($r2));
- return ($r1, $r2);
- }
- sub siqs_factor_from_square ($n, $square_indices, $smooth_relations) {
- # Given one of the solutions returned by siqs_solve_matrix_opt,
- # return the factor f determined by f = gcd(a - b, n), where
- # a, b are calculated from the solution such that a*a = b*b (mod n).
- # Return f, a factor of n (possibly a trivial one).
- my ($sqrt1, $sqrt2) = siqs_calc_sqrts($n, $square_indices, $smooth_relations);
- #(($sqrt1 * $sqrt1) % $n == ($sqrt2 * $sqrt2) % $n) or die 'error';
- return Math::GMPz->new(gcd($sqrt1 - $sqrt2, $n));
- }
- sub siqs_find_more_factors_gcd (@numbers) {
- my %res;
- foreach my $i (0 .. $#numbers) {
- my $n = $numbers[$i];
- $res{$n} = $n;
- foreach my $k ($i + 1 .. $#numbers) {
- my $m = $numbers[$k];
- my $fact = Math::GMPz->new(gcd($n, $m));
- if ($fact != 1 and $fact != $n and $fact != $m) {
- if (not exists($res{$fact})) {
- say "SIQS: GCD found non-trivial factor: $fact";
- $res{$fact} = $fact;
- }
- my $t1 = $n / $fact;
- my $t2 = $m / $fact;
- $res{$t1} = $t1;
- $res{$t2} = $t2;
- }
- }
- }
- return (values %res);
- }
- sub siqs_find_factors ($n, $perfect_squares, $smooth_relations) {
- # Perform the last step of the Self-Initializing Quadratic Field.
- # Given the solutions returned by siqs_solve_matrix_opt, attempt to
- # identify a number of (not necessarily prime) factors of n, and
- # return them.
- my @factors;
- my $rem = $n;
- my %non_prime_factors;
- my %prime_factors;
- foreach my $square_indices (@$perfect_squares) {
- my $fact = siqs_factor_from_square($n, $square_indices, $smooth_relations);
- if ($fact > 1 and $fact < $rem) {
- if (is_prime($fact)) {
- if (not exists $prime_factors{$fact}) {
- say "SIQS: Prime factor found: $fact";
- $prime_factors{$fact} = $fact;
- }
- $rem = check_factor($rem, $fact, \@factors);
- if ($rem == 1) {
- last;
- }
- if (is_prime($rem)) {
- push @factors, $rem;
- $rem = 1;
- last;
- }
- if (defined(my $root = check_perfect_power($rem))) {
- say "SIQS: Perfect power detected with root: $root";
- push @factors, $root;
- $rem = 1;
- last;
- }
- }
- else {
- if (not exists $non_prime_factors{$fact}) {
- say "SIQS: Composite factor found: $fact";
- $non_prime_factors{$fact} = $fact;
- }
- }
- }
- }
- if ($rem != 1 and keys(%non_prime_factors)) {
- $non_prime_factors{$rem} = $rem;
- my @primes;
- my @composites;
- foreach my $fact (siqs_find_more_factors_gcd(values %non_prime_factors)) {
- if (is_prime($fact)) {
- push @primes, $fact;
- }
- elsif ($fact > 1) {
- push @composites, $fact;
- }
- }
- foreach my $fact (@primes, @composites) {
- if ($fact != $rem and $rem % $fact == 0) {
- say "SIQS: Using non-trivial factor from GCD: $fact";
- $rem = check_factor($rem, $fact, \@factors);
- }
- if ($rem == 1 or is_prime($rem)) {
- last;
- }
- }
- }
- if ($rem != 1) {
- push @factors, $rem;
- }
- return @factors;
- }
- sub siqs_choose_range ($n) {
- # Choose m for sieving in [-m, m].
- $n = "$n";
- return sprintf('%.0f', exp(sqrt(log($n) * log(log($n))) / 2));
- }
- sub siqs_choose_nf ($n) {
- # Choose parameters nf (sieve of factor base)
- $n = "$n";
- return sprintf('%.0f', exp(sqrt(log($n) * log(log($n))))**(sqrt(2) / 4));
- }
- sub siqs_choose_nf2 ($n) {
- # Choose parameters nf (sieve of factor base)
- $n = "$n";
- return sprintf('%.0f', exp(sqrt(log($n) * log(log($n))) / 2));
- }
- sub siqs_factorize ($n, $nf) {
- # Use the Self-Initializing Quadratic Sieve algorithm to identify
- # one or more non-trivial factors of the given number n. Return the
- # factors as a list.
- my $m = siqs_choose_range($n);
- my @factors;
- my $factor_base = siqs_factor_base_primes($n, $nf);
- my $factor_prod = Math::GMPz->new(vecprod(map { $_->{p} } @$factor_base));
- my %factor_base_index;
- @factor_base_index{map { $_->{p} } @{$factor_base}} = 0 .. $#{$factor_base};
- my $factor_base_info = {
- base => $factor_base,
- prod => $factor_prod,
- index => \%factor_base_index,
- };
- my $smooth_relations = [];
- my $required_relations_ratio = 1;
- my $success = 0;
- my $prev_cnt = 0;
- my $i_poly = 0;
- my ($g, $h, $arr);
- while (not $success) {
- say "*** Step 1/2: Finding smooth relations ***";
- say "SIQS sieving range: [-$m, $m]";
- my $required_relations = sprintf('%.0f', (scalar(@$factor_base) + 1) * $required_relations_ratio);
- say "Target: $required_relations relations.";
- my $enough_relations = 0;
- while (not $enough_relations) {
- if ($i_poly == 0) {
- ($g, $h, $arr) = siqs_find_first_poly($n, $m, $factor_base);
- }
- else {
- ($g, $h) = siqs_find_next_poly($n, $factor_base, $i_poly, $g, $arr);
- }
- if (++$i_poly >= (1 << $#{$arr})) {
- $i_poly = 0;
- }
- my $sieve_array = siqs_sieve($factor_base, $m);
- $enough_relations = siqs_trial_division($n, $sieve_array, $factor_base_info, $smooth_relations, $g, $h, $m, $required_relations);
- if ( scalar(@$smooth_relations) >= $required_relations
- or scalar(@$smooth_relations) > $prev_cnt) {
- printf("Progress: %d/%d relations.\r", scalar(@$smooth_relations), $required_relations);
- $prev_cnt = scalar(@$smooth_relations);
- }
- }
- say "\n\n*** Step 2/2: Linear Algebra ***";
- say "Building matrix for linear algebra step...";
- my $M = siqs_build_matrix($factor_base, $smooth_relations);
- my ($M_opt, $M_n, $M_m) = siqs_build_matrix_opt($M);
- say "Finding perfect squares using Gaussian elimination...";
- my $perfect_squares = siqs_solve_matrix_opt($M_opt, $M_n, $M_m);
- say "Finding factors from congruences of squares...\n";
- @factors = siqs_find_factors($n, $perfect_squares, $smooth_relations);
- if (scalar(@factors) > 1) {
- $success = 1;
- }
- else {
- say "Failed to find a solution. Finding more relations...";
- $required_relations_ratio += 0.05;
- }
- }
- return @factors;
- }
- sub check_factor ($n, $i, $factors) {
- while ($n % $i == 0) {
- $n /= $i;
- push @$factors, $i;
- if (is_prime($n)) {
- push @$factors, $n;
- return 1;
- }
- }
- return $n;
- }
- sub trial_division_small_primes ($n) {
- # Perform trial division on the given number n using all primes up
- # to upper_bound. Initialize the global variable small_primes with a
- # list of all primes <= upper_bound. Return (factors, rem), where
- # factors is the list of identified prime factors of n, and rem is the
- # remaining factor. If rem = 1, the function terminates early, without
- # fully initializing small_primes.
- say "[*] Trial division...";
- my $factors = [];
- my $rem = $n;
- foreach my $p (@small_primes) {
- if (Math::GMPz::Rmpz_divisible_ui_p($rem, $p)) {
- $rem = check_factor($rem, $p, $factors);
- last if ($rem == 1);
- }
- }
- return ($factors, $rem);
- }
- sub fast_fibonacci_factor ($n, $upto) {
- my $g = Math::GMPz::Rmpz_init();
- my ($P, $Q) = (3, 1);
- my $U0 = Math::GMPz::Rmpz_init_set_ui(0);
- my $U1 = Math::GMPz::Rmpz_init_set_ui(1);
- my $V0 = Math::GMPz::Rmpz_init_set_ui(2);
- my $V1 = Math::GMPz::Rmpz_init_set_ui($P);
- foreach my $k (2 .. $upto) {
- # my ($U, $V) = Math::Prime::Util::GMP::lucas_sequence($n, $P, $Q, $k);
- Math::GMPz::Rmpz_set($g, $U1);
- Math::GMPz::Rmpz_mul_ui($U1, $U1, $P);
- Math::GMPz::Rmpz_submul_ui($U1, $U0, $Q);
- Math::GMPz::Rmpz_mod($U1, $U1, $n);
- Math::GMPz::Rmpz_set($U0, $g);
- Math::GMPz::Rmpz_set($g, $V1);
- Math::GMPz::Rmpz_mul_ui($V1, $V1, $P);
- Math::GMPz::Rmpz_submul_ui($V1, $V0, $Q);
- Math::GMPz::Rmpz_mod($V1, $V1, $n);
- Math::GMPz::Rmpz_set($V0, $g);
- foreach my $param ([$U1, 0], [$V1, -$P, 0]) {
- my ($t, @deltas) = @$param;
- foreach my $delta (@deltas) {
- ($delta >= 0)
- ? Math::GMPz::Rmpz_add_ui($g, $t, $delta)
- : Math::GMPz::Rmpz_sub_ui($g, $t, -$delta);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($g, 12345) > 0
- and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- return $g;
- }
- }
- }
- }
- return undef;
- }
- sub fast_power_check ($n, $upto) {
- state $t = Math::GMPz::Rmpz_init_nobless();
- state $g = Math::GMPz::Rmpz_init_nobless();
- my $base_limit = vecmin(logint($n, 2), 150);
- foreach my $base (2 .. $base_limit) {
- Math::GMPz::Rmpz_set_ui($t, $base);
- foreach my $exp (2 .. $upto) {
- Math::GMPz::Rmpz_mul_ui($t, $t, $base);
- foreach my $k ($base <= 10 ? (1 .. ($base_limit >> 1)) : 1) {
- Math::GMPz::Rmpz_mul_ui($g, $t, $k);
- Math::GMPz::Rmpz_sub_ui($g, $g, 1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0 and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- my $r = Math::GMPz::Rmpz_init_set($g);
- return $r if ($r > 12345);
- }
- Math::GMPz::Rmpz_mul_ui($g, $t, $k);
- Math::GMPz::Rmpz_add_ui($g, $g, 1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0 and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- my $r = Math::GMPz::Rmpz_init_set($g);
- return $r if ($r > 12345);
- }
- }
- }
- }
- return undef;
- }
- sub cyclotomic_polynomial ($n, $x) {
- my @d;
- foreach my $p (map { $_->[0] } factor_exp($n)) {
- push @d, map { [$_->[0] * $p, $_->[1] + 1] } @d;
- push @d, [$p, 1];
- }
- push @d, [1, 0];
- my $t = Math::GMPz::Rmpz_init();
- my $u = Math::GMPz::Rmpz_init_set_ui(1);
- my $v = Math::GMPz::Rmpz_init_set_ui(1);
- foreach my $pp (@d) {
- Math::GMPz::Rmpz_ui_pow_ui($t, $x, int($n / $pp->[0]));
- Math::GMPz::Rmpz_sub_ui($t, $t, 1);
- if ($pp->[1] % 2 == 0) {
- Math::GMPz::Rmpz_mul($u, $u, $t);
- }
- else {
- Math::GMPz::Rmpz_mul($v, $v, $t);
- }
- }
- Math::GMPz::Rmpz_divexact($u, $u, $v);
- return $u;
- }
- sub cyclotomic_factorization ($n) {
- my $g = Math::GMPz::Rmpz_init();
- my $base_limit = vecmin(1 + logint($n, 2), 500);
- for (my $base = $base_limit ; $base >= 2 ; $base -= 1) {
- my $lim = 1 + logint($n, $base);
- foreach my $k (1 .. $lim) {
- my $c = cyclotomic_polynomial($k, $base);
- Math::GMPz::Rmpz_gcd($g, $n, $c);
- if ( Math::GMPz::Rmpz_cmp_ui($g, 1) > 0
- and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- return $g if ($g > 1e7);
- }
- }
- }
- return undef;
- }
- sub fast_lucasVmod ($P, $n, $m) { # assumes Q = 1
- my ($V1, $V2) = (Math::GMPz::Rmpz_init_set_ui(2), Math::GMPz::Rmpz_init_set($P));
- my ($Q1, $Q2) = (Math::GMPz::Rmpz_init_set_ui(1), Math::GMPz::Rmpz_init_set_ui(1));
- foreach my $bit (ntheory::todigits($n, 2)) {
- Math::GMPz::Rmpz_mul($Q1, $Q1, $Q2);
- Math::GMPz::Rmpz_mod($Q1, $Q1, $m);
- if ($bit) {
- Math::GMPz::Rmpz_mul($V1, $V1, $V2);
- Math::GMPz::Rmpz_powm_ui($V2, $V2, 2, $m);
- Math::GMPz::Rmpz_submul($V1, $Q1, $P);
- Math::GMPz::Rmpz_submul_ui($V2, $Q2, 2);
- Math::GMPz::Rmpz_mod($V1, $V1, $m);
- }
- else {
- Math::GMPz::Rmpz_set($Q2, $Q1);
- Math::GMPz::Rmpz_mul($V2, $V2, $V1);
- Math::GMPz::Rmpz_powm_ui($V1, $V1, 2, $m);
- Math::GMPz::Rmpz_submul($V2, $Q1, $P);
- Math::GMPz::Rmpz_submul_ui($V1, $Q2, 2);
- Math::GMPz::Rmpz_mod($V2, $V2, $m);
- }
- }
- Math::GMPz::Rmpz_mod($V1, $V1, $m);
- return $V1;
- }
- sub chebyshev_factorization ($n, $B, $A = 127) {
- # The Chebyshev factorization method, taking
- # advantage of the smoothness of p-1 or p+1.
- my $x = Math::GMPz::Rmpz_init_set_ui($A);
- my $i = Math::GMPz::Rmpz_init_set_ui(2);
- Math::GMPz::Rmpz_invert($i, $i, $n);
- my sub chebyshevTmod ($A, $x) {
- Math::GMPz::Rmpz_mul_2exp($x, $x, 1);
- Math::GMPz::Rmpz_set($x, fast_lucasVmod($x, $A, $n));
- Math::GMPz::Rmpz_mul($x, $x, $i);
- Math::GMPz::Rmpz_mod($x, $x, $n);
- }
- my $g = Math::GMPz::Rmpz_init();
- my $lnB = log($B);
- foreach my $p (sieve_primes(2, $B)) {
- chebyshevTmod($p**int($lnB / log($p)), $x); # T_k(x) (mod n)
- Math::GMPz::Rmpz_sub_ui($g, $x, 1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
- return undef if (Math::GMPz::Rmpz_cmp($g, $n) == 0);
- return $g;
- }
- }
- return undef;
- }
- sub fibonacci_factorization ($n, $bound) {
- # The Fibonacci factorization method, taking
- # advantage of the smoothness of `p - legendre(p, 5)`.
- my ($P, $Q) = (1, 0);
- for (my $k = 2 ; ; ++$k) {
- my $D = (-1)**$k * (2 * $k + 1);
- if (Math::GMPz::Rmpz_si_kronecker($D, $n) == -1) {
- $Q = (1 - $D) / 4;
- last;
- }
- }
- state %cache;
- my $g = Math::GMPz::Rmpz_init();
- for (; ;) {
- return undef if $bound <= 1;
- my $d = ($cache{$bound} //= consecutive_integer_lcm($bound));
- my ($U, $V) = map { Math::GMPz::Rmpz_init_set_str($_, 10) } lucas_sequence($n, $P, $Q, $d);
- foreach my $t ($U, $V - 2, $V, $V + 2) {
- Math::GMPz::Rmpz_gcd($g, $t, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($g, 1) > 0
- and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- return $g;
- }
- }
- if ($U == 0) {
- say ":: p±1 seems to be $bound-smooth...";
- $bound >>= 1;
- next;
- }
- return undef;
- #say "=> Lucas p±1...";
- #return lucas_factorization($n, Math::GMPz::Rmpz_init_set_str($d, 10));
- }
- }
- sub lucas_factorization ($n, $d) {
- # The Lucas factorization method, taking
- # advantage of the smoothness of p-1 or p+1.
- my $Q;
- for (my $k = 2 ; ; ++$k) {
- my $D = (-1)**$k * (2 * $k + 1);
- if (Math::GMPz::Rmpz_si_kronecker($D, $n) == -1) {
- $Q = (1 - $D) / 4;
- last;
- }
- }
- my $s = Math::GMPz::Rmpz_scan1($d, 0);
- my $U1 = Math::GMPz::Rmpz_init_set_ui(1);
- my ($V1, $V2) = (Math::GMPz::Rmpz_init_set_ui(2), Math::GMPz::Rmpz_init_set_ui(1));
- my ($Q1, $Q2) = (Math::GMPz::Rmpz_init_set_ui(1), Math::GMPz::Rmpz_init_set_ui(1));
- foreach my $bit (split(//, substr(Math::GMPz::Rmpz_get_str($d, 2), 0, -$s - 1))) {
- Math::GMPz::Rmpz_mul($Q1, $Q1, $Q2);
- Math::GMPz::Rmpz_mod($Q1, $Q1, $n);
- if ($bit) {
- Math::GMPz::Rmpz_mul_si($Q2, $Q1, $Q);
- Math::GMPz::Rmpz_mul($U1, $U1, $V2);
- Math::GMPz::Rmpz_mul($V1, $V1, $V2);
- Math::GMPz::Rmpz_powm_ui($V2, $V2, 2, $n);
- Math::GMPz::Rmpz_sub($V1, $V1, $Q1);
- Math::GMPz::Rmpz_submul_ui($V2, $Q2, 2);
- Math::GMPz::Rmpz_mod($V1, $V1, $n);
- Math::GMPz::Rmpz_mod($U1, $U1, $n);
- }
- else {
- Math::GMPz::Rmpz_set($Q2, $Q1);
- Math::GMPz::Rmpz_mul($U1, $U1, $V1);
- Math::GMPz::Rmpz_mul($V2, $V2, $V1);
- Math::GMPz::Rmpz_sub($U1, $U1, $Q1);
- Math::GMPz::Rmpz_powm_ui($V1, $V1, 2, $n);
- Math::GMPz::Rmpz_sub($V2, $V2, $Q1);
- Math::GMPz::Rmpz_submul_ui($V1, $Q2, 2);
- Math::GMPz::Rmpz_mod($V2, $V2, $n);
- Math::GMPz::Rmpz_mod($U1, $U1, $n);
- }
- }
- Math::GMPz::Rmpz_mul($Q1, $Q1, $Q2);
- Math::GMPz::Rmpz_mul_si($Q2, $Q1, $Q);
- Math::GMPz::Rmpz_mul($U1, $U1, $V1);
- Math::GMPz::Rmpz_mul($V1, $V1, $V2);
- Math::GMPz::Rmpz_sub($U1, $U1, $Q1);
- Math::GMPz::Rmpz_sub($V1, $V1, $Q1);
- Math::GMPz::Rmpz_mul($Q1, $Q1, $Q2);
- my $t = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_gcd($t, $U1, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($t, 1) > 0
- and Math::GMPz::Rmpz_cmp($t, $n) < 0) {
- return $t;
- }
- Math::GMPz::Rmpz_gcd($t, $V1, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($t, 1) > 0
- and Math::GMPz::Rmpz_cmp($t, $n) < 0) {
- return $t;
- }
- for (1 .. $s) {
- Math::GMPz::Rmpz_mul($U1, $U1, $V1);
- Math::GMPz::Rmpz_mod($U1, $U1, $n);
- Math::GMPz::Rmpz_powm_ui($V1, $V1, 2, $n);
- Math::GMPz::Rmpz_submul_ui($V1, $Q1, 2);
- Math::GMPz::Rmpz_powm_ui($Q1, $Q1, 2, $n);
- Math::GMPz::Rmpz_gcd($t, $U1, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($t, 1) > 0
- and Math::GMPz::Rmpz_cmp($t, $n) < 0) {
- return $t;
- }
- Math::GMPz::Rmpz_gcd($t, $V1, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($t, 1) > 0
- and Math::GMPz::Rmpz_cmp($t, $n) < 0) {
- return $t;
- }
- }
- return undef;
- }
- sub pollard_pm1_lcm_find_factor ($n, $bound) {
- # Pollard p-1 method (LCM).
- my $g = Math::GMPz::Rmpz_init();
- my $t = Math::GMPz::Rmpz_init_set_ui(random_prime(1e6));
- foreach my $p (sieve_primes(2, $bound)) {
- Math::GMPz::Rmpz_powm_ui($t, $t, $p**int(log(ULONG_MAX >> 32) / log($p)), $n);
- Math::GMPz::Rmpz_sub_ui($g, $t, 1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
- return undef if ($g == $n);
- return $g;
- }
- }
- return undef;
- }
- sub pollard_pm1_factorial_find_factor ($n, $bound2) {
- # Pollard p-1 method (factorial).
- my $bound1 = 1e5;
- state %cache;
- my $g = Math::GMPz::Rmpz_init();
- my $t = Math::GMPz::Rmpz_init_set_ui(random_prime(1e6));
- if (exists $cache{$n}) {
- $t = $cache{$n}{value};
- $bound1 = $cache{$n}{bound};
- }
- else {
- foreach my $k (2 .. $bound1) {
- Math::GMPz::Rmpz_powm_ui($t, $t, $k, $n);
- Math::GMPz::Rmpz_sub_ui($g, $t, 1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
- return undef if ($g == $n);
- return $g;
- }
- }
- }
- while ($bound1 >= $bound2) {
- $bound2 *= 2;
- }
- foreach my $p (sieve_primes($bound1, $bound2)) {
- Math::GMPz::Rmpz_powm_ui($t, $t, $p, $n);
- Math::GMPz::Rmpz_sub_ui($g, $t, 1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
- return undef if ($g == $n);
- return $g;
- }
- }
- $cache{$n}{value} = $t;
- $cache{$n}{bound} = $bound2 + 1;
- return undef;
- }
- sub pollard_rho_find_factor ($n, $max_iter) {
- # Pollard rho method, using the polynomial:
- # f(x) = x^2 - 1, with x_0 = 1+floor(log_2(n)).
- state %cache;
- my $u = logint($n, 2) + 1;
- my $x = Math::GMPz::Rmpz_init_set_ui($u);
- my $y = Math::GMPz::Rmpz_init_set_ui($u * $u - 1);
- if (exists $cache{$n}) {
- $x = $cache{$n}{x};
- $y = $cache{$n}{y};
- }
- my $g = Math::GMPz::Rmpz_init();
- for (1 .. $max_iter) {
- # f(x) = x^2 - 1
- Math::GMPz::Rmpz_powm_ui($x, $x, 2, $n);
- Math::GMPz::Rmpz_sub_ui($x, $x, 1);
- # f(f(x)) = (x^2 - 1)^2 - 1 = (x^2 - 2) * x^2
- Math::GMPz::Rmpz_powm_ui($g, $y, 2, $n);
- Math::GMPz::Rmpz_sub_ui($y, $g, 2);
- Math::GMPz::Rmpz_mul($y, $y, $g);
- Math::GMPz::Rmpz_sub($g, $x, $y);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) != 0) {
- return undef if ($g == $n);
- return $g;
- }
- }
- $cache{$n}{x} = $x;
- $cache{$n}{y} = $y;
- return undef;
- }
- sub pollard_pm1_ntheory_factor ($n, $max_iter) {
- my ($p, $q) = Math::Prime::Util::GMP::pminus1_factor($n, $max_iter);
- return $p if defined($q);
- return undef;
- #return pollard_pm1_factorial_find_factor($n, $max_iter);
- }
- sub williams_pp1_ntheory_factor ($n, $max_iter) {
- my ($p, $q) = Math::Prime::Util::GMP::pplus1_factor($n, $max_iter);
- return $p if defined($q);
- return undef;
- }
- sub pollard_rho_ntheory_factor ($n, $max_iter) {
- my ($p, $q) =
- (rand(1) < 0.5)
- ? (Math::Prime::Util::GMP::prho_factor($n, $max_iter))
- : (Math::Prime::Util::GMP::pbrent_factor($n, $max_iter));
- return $p if defined($q);
- return undef;
- #return pollard_rho_find_factor($n, $max_iter >> 1);
- }
- sub pollard_rho_sqrt_find_factor ($n, $max_iter) {
- # Pollard rho method, using the polynomial:
- # f(x) = x^2 + c
- #
- # where
- # c = floor(sqrt(n)) - (floor(sqrt(n))^2 - n)
- # c = n + s - s^2, with s = floor(sqrt(n))
- #
- # and
- # x_0 = 3^2 + c
- my $s = Math::GMPz->new(sqrtint($n));
- my $c = $n + $s - $s * $s;
- my $a0 = 3;
- my $a1 = ($a0 * $a0 + $c);
- my $a2 = ($a1 * $a1 + $c);
- my $g = Math::GMPz::Rmpz_init();
- for (1 .. $max_iter) {
- Math::GMPz::Rmpz_sub($g, $a2, $a1);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) != 0) {
- return undef if ($g == $n);
- return $g;
- }
- Math::GMPz::Rmpz_powm_ui($a1, $a1, 2, $n);
- Math::GMPz::Rmpz_add($a1, $a1, $c);
- Math::GMPz::Rmpz_powm_ui($a2, $a2, 2, $n);
- Math::GMPz::Rmpz_add($a2, $a2, $c);
- Math::GMPz::Rmpz_powm_ui($a2, $a2, 2, $n);
- Math::GMPz::Rmpz_add($a2, $a2, $c);
- }
- return undef;
- }
- sub pollard_rho_exp_find_factor ($n, $max_iter) {
- my $B = logint($n, 5)**2;
- if ($B > 50_000) {
- $B = 50_000;
- }
- my $e = Math::GMPz::Rmpz_init_set_str(consecutive_integer_lcm($B), 10);
- my $c = 2 * $e - 1;
- my $x = Math::GMPz::Rmpz_init_set_ui(1);
- my $y = Math::GMPz::Rmpz_init();
- my $g = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_powm($x, $x, $e, $n);
- Math::GMPz::Rmpz_add($x, $x, $c);
- Math::GMPz::Rmpz_mod($x, $x, $n);
- Math::GMPz::Rmpz_powm($y, $x, $e, $n);
- Math::GMPz::Rmpz_add($y, $y, $c);
- Math::GMPz::Rmpz_mod($y, $y, $n);
- for (1 .. $max_iter) {
- Math::GMPz::Rmpz_powm($x, $x, $e, $n);
- Math::GMPz::Rmpz_add($x, $x, $c);
- Math::GMPz::Rmpz_mod($x, $x, $n);
- Math::GMPz::Rmpz_powm($y, $y, $e, $n);
- Math::GMPz::Rmpz_add($y, $y, $c);
- Math::GMPz::Rmpz_mod($y, $y, $n);
- Math::GMPz::Rmpz_powm($y, $y, $e, $n);
- Math::GMPz::Rmpz_add($y, $y, $c);
- Math::GMPz::Rmpz_mod($y, $y, $n);
- Math::GMPz::Rmpz_sub($g, $x, $y);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) != 0) {
- return undef if ($g == $n);
- return $g;
- }
- }
- return undef;
- }
- sub phi_finder_factor ($n, $max_iter) {
- # Phi-finder algorithm for semiprimes, due to Kyle Kloster (2010)
- my $E = $n - 2 * Math::GMPz->new(sqrtint($n)) + 1;
- my $E0 = Math::GMPz->new(powmod(2, -$E, $n));
- my $L = logint($n, 2);
- my $i = 0;
- while (Math::GMPz::Rmpz_scan1($E0, 0) < Math::GMPz::Rmpz_sizeinbase($E0, 2) - 1) {
- Math::GMPz::Rmpz_mul_2exp($E0, $E0, $L);
- Math::GMPz::Rmpz_mod($E0, $E0, $n);
- ++$i;
- return undef if ($i > $max_iter);
- }
- my $t = 0;
- foreach my $k (0 .. $L) {
- if (Math::GMPz->new(powmod(2, $k, $n)) == $E0) {
- $t = $k;
- last;
- }
- }
- my $phi = abs($i * $L - $E - $t);
- my $q = ($n - $phi + 1);
- my $p = ($q + Math::GMPz->new(sqrtint(abs($q * $q - 4 * $n)))) >> 1;
- (($n % $p) == 0) ? $p : undef;
- }
- sub fermat_find_factor ($n, $max_iter) {
- # Fermat's factorization method, trying to represent `n` as a difference of two squares:
- # n = a^2 - b^2, where n = (a-b) * (a+b).
- my $p = Math::GMPz::Rmpz_init(); # p = floor(sqrt(n))
- my $q = Math::GMPz::Rmpz_init(); # q = p^2 - n
- Math::GMPz::Rmpz_sqrtrem($p, $q, $n);
- Math::GMPz::Rmpz_neg($q, $q);
- for (my $j = 1 ; $j <= $max_iter ; ++$j) {
- Math::GMPz::Rmpz_addmul_ui($q, $p, 2);
- Math::GMPz::Rmpz_add_ui($q, $q, 1);
- Math::GMPz::Rmpz_add_ui($p, $p, 1);
- if (Math::GMPz::Rmpz_perfect_square_p($q)) {
- Math::GMPz::Rmpz_sqrt($q, $q);
- my $r = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_sub($r, $p, $q);
- return $r if ($r > 12345);
- }
- }
- return undef;
- }
- sub holf_find_factor ($n, $max_iter) {
- # Hart’s One-Line Factoring Algorithm
- my $m = Math::GMPz::Rmpz_init();
- my $s = Math::GMPz::Rmpz_init();
- foreach my $i (1 .. $max_iter) {
- Math::GMPz::Rmpz_mul_ui($s, $n, 4 * $i);
- Math::GMPz::Rmpz_sqrt($s, $s);
- Math::GMPz::Rmpz_add_ui($s, $s, 1);
- Math::GMPz::Rmpz_mul($m, $s, $s);
- Math::GMPz::Rmpz_mod($m, $m, $n);
- if (Math::GMPz::Rmpz_perfect_square_p($m)) {
- Math::GMPz::Rmpz_sqrt($m, $m);
- Math::GMPz::Rmpz_sub($m, $s, $m);
- Math::GMPz::Rmpz_gcd($m, $m, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($m, 1) > 0
- and Math::GMPz::Rmpz_cmp($m, $n) < 0) {
- return $m if ($m > 12345);
- }
- }
- }
- return undef;
- }
- sub miller_rabin_factor ($n, $tries) {
- # Miller-Rabin factorization method.
- # https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
- my $D = $n - 1;
- my $s = Math::GMPz::Rmpz_scan1($D, 0);
- my $r = $s - 1;
- my $d = $D >> $s;
- if ($s > 20 and $tries > 10) {
- $tries = 10;
- }
- my $x = Math::GMPz::Rmpz_init();
- my $g = Math::GMPz::Rmpz_init();
- for (1 .. $tries) {
- my $p = random_prime(1e7);
- Math::GMPz::Rmpz_powm($x, Math::GMPz::Rmpz_init_set_ui($p), $d, $n);
- foreach my $k (0 .. $r) {
- last if (Math::GMPz::Rmpz_cmp_ui($x, 1) == 0);
- last if (Math::GMPz::Rmpz_cmp($x, $D) == 0);
- foreach my $i (1, -1) {
- Math::GMPz::Rmpz_gcd($g, $x + $i, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($g, 1) > 0
- and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- if ($g > 12345) {
- return $g;
- }
- }
- }
- Math::GMPz::Rmpz_powm_ui($x, $x, 2, $n);
- }
- }
- return undef;
- }
- sub lucas_miller_factor ($n, $I, $tries) {
- # Lucas-Miller factorization method.
- my $D = $n + $I;
- my $s = Math::GMPz::Rmpz_scan1($D, 0);
- my $r = $s;
- my $d = $D >> $s;
- if ($s > 10 and $tries > 5) {
- $tries = 5;
- }
- my $x = Math::GMPz::Rmpz_init();
- my $g = Math::GMPz::Rmpz_init();
- foreach my $i (1 .. $tries) {
- my $P = 1 + int(rand(1e6));
- my $Q = 1 + int(rand(1e6));
- $Q *= -1 if (rand(1) < 0.5);
- next if ntheory::is_square($P * $P - 4 * $Q);
- Math::GMPz::Rmpz_set($x, $d);
- foreach my $k (0 .. $r) {
- my ($U, $V) = lucas_sequence($n, $P, $Q, $x);
- foreach my $t ($U, $V) {
- Math::GMPz::Rmpz_set_str($g, $t, 10);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($g, 1) > 0
- and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- return $g if ($g > 12345);
- }
- }
- Math::GMPz::Rmpz_set_str($g, $V, 10);
- Math::GMPz::Rmpz_sub_ui($g, $g, $P);
- Math::GMPz::Rmpz_gcd($g, $g, $n);
- if ( Math::GMPz::Rmpz_cmp_ui($g, 1) > 0
- and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
- return $g if ($g > 12345);
- }
- Math::GMPz::Rmpz_mul_2exp($x, $x, 1);
- }
- }
- return undef;
- }
- sub simple_cfrac_find_factor ($n, $max_iter) {
- # Simple version of the continued-fraction factorization method.
- # Efficient for numbers that have factors relatively close to sqrt(n)
- my $x = Math::GMPz::Rmpz_init();
- my $y = Math::GMPz::Rmpz_init();
- my $z = Math::GMPz::Rmpz_init_set_ui(1);
- my $t = Math::GMPz::Rmpz_init();
- my $w = Math::GMPz::Rmpz_init();
- my $r = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_sqrt($x, $n);
- Math::GMPz::Rmpz_set($y, $x);
- Math::GMPz::Rmpz_add($w, $x, $x);
- Math::GMPz::Rmpz_set($r, $w);
- my $f2 = Math::GMPz::Rmpz_init_set($x);
- my $f1 = Math::GMPz::Rmpz_init_set_ui(1);
- foreach (1 .. $max_iter) {
- # y = r*z - y
- Math::GMPz::Rmpz_mul($t, $r, $z);
- Math::GMPz::Rmpz_sub($y, $t, $y);
- # z = (n - y*y) / z
- Math::GMPz::Rmpz_mul($t, $y, $y);
- Math::GMPz::Rmpz_sub($t, $n, $t);
- Math::GMPz::Rmpz_divexact($z, $t, $z);
- # r = (x + y) / z
- Math::GMPz::Rmpz_add($t, $x, $y);
- Math::GMPz::Rmpz_div($r, $t, $z);
- # f1 = (f1 + r*f2) % n
- Math::GMPz::Rmpz_addmul($f1, $f2, $r);
- Math::GMPz::Rmpz_mod($f1, $f1, $n);
- # swap f1 with f2
- ($f1, $f2) = ($f2, $f1);
- if (Math::GMPz::Rmpz_perfect_square_p($z)) {
- my $g = Math::GMPz->new(gcd($f1 - Math::GMPz->new(sqrtint($z)), $n));
- if ($g > 1 and $g < $n) {
- return $g if ($g > 12345);
- }
- }
- last if ($z == 1);
- }
- return undef;
- }
- sub store_factor ($rem, $f, $factors) {
- $f || return;
- if (ref($f) ne 'Math::GMPz') {
- $f =~ /^[0-9]+\z/ or return;
- $f = Math::GMPz->new($f);
- }
- $f < $$rem or return;
- $$rem % $f == 0 or die 'error';
- if (is_prime($f)) {
- say("`-> prime factor: ", $f);
- $$rem = check_factor($$rem, $f, $factors);
- }
- else {
- say("`-> composite factor: ", $f);
- $$rem /= $f;
- # Try to find a small factor of f
- my $f_rem = find_small_factors($f, $factors);
- if ($f_rem < $f) {
- $$rem *= $f_rem;
- }
- else {
- # Use SIQS to factorize f
- find_prime_factors($f, $factors);
- foreach my $p (@$factors) {
- if ($$rem % $p == 0) {
- $$rem = check_factor($$rem, $p, $factors);
- last if $$rem == 1;
- }
- }
- }
- }
- return 1;
- }
- sub find_small_factors ($rem, $factors) {
- # Some special-purpose factorization methods to attempt to find small prime factors.
- # Collect the identified prime factors in the `$factors` array and return 1 if all
- # prime factors were found, or otherwise the remaining factor.
- my %state = (
- cyclotomic_check => 1,
- fast_power_check => 1,
- fast_fibonacci_check => 1,
- );
- my $len = length($rem);
- if (0) {
- foreach my $k (1 .. 10000) {
- #use Math::Sidef qw(cyclotomic);
- #my $c = cyclotomic($k, 2);
- my $c = Math::GMPz->new(2)**$k - 1;
- if (gcd($rem, $c) == $rem and $c != $rem) {
- my @factors = factorize($c);
- foreach my $f (@factors) {
- store_factor(\$rem, $f, $factors);
- }
- if (is_prime($rem)) {
- push(@$factors, $rem);
- $rem = 1;
- }
- if ($rem == 1) {
- return 1;
- }
- }
- }
- }
- my @factorization_methods = (
- sub {
- say "=> Miller-Rabin method...";
- miller_rabin_factor($rem, ($len > 1000) ? 15 : MILLER_RABIN_ITERATIONS);
- },
- sub {
- if ($len < 3000) {
- say "=> Lucas-Miller method (n+1)...";
- lucas_miller_factor($rem, +1, ($len > 1000) ? 15 : LUCAS_MILLER_ITERATIONS);
- }
- },
- sub {
- if ($len < 3000) {
- say "=> Lucas-Miller method (n-1)...";
- lucas_miller_factor($rem, -1, ($len > 1000) ? 15 : LUCAS_MILLER_ITERATIONS);
- }
- },
- #~ sub {
- #~ say "=> Phi finder method...";
- #~ phi_finder_factor($rem, PHI_FINDER_ITERATIONS);
- #~ },
- #~ sub {
- #~ say "=> Fermat's method...";
- #~ fermat_find_factor($rem, FERMAT_ITERATIONS);
- #~ },
- sub {
- say "=> HOLF method...";
- holf_find_factor($rem, HOLF_ITERATIONS);
- },
- #~ sub {
- #~ say "=> CFRAC simple...";
- #~ simple_cfrac_find_factor($rem, CFRAC_ITERATIONS);
- #~ },
- sub {
- $state{fast_fibonacci_check} || return undef;
- say "=> Fast Fibonacci check...";
- my $f = fast_fibonacci_factor($rem, 5000);
- $f // do { $state{fast_fibonacci_check} = 0 };
- $f;
- },
- #~ sub {
- #~ say "=> Pollard rho (11M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e11));
- #~ },
- sub {
- say "=> Pollard p-1 (500K)...";
- pollard_pm1_ntheory_factor($rem, 500_000);
- },
- #~ sub {
- #~ say "=> Williams p±1 (500K)...";
- #~ williams_pp1_ntheory_factor($rem, 500_000);
- #~ },
- #~ sub {
- #~ if ($len < 1000) {
- #~ say "=> Chebyshev p±1 (500K)...";
- #~ chebyshev_factorization($rem, 500_000, int(rand(1e6)) + 2);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len < 1000) {
- #~ say "=> Chebyshev p±1 (1M)...";
- #~ chebyshev_factorization($rem, 1_000_000, int(rand(1e6)) + 2);
- #~ }
- #~ },
- #~ sub {
- #~ say "=> Williams p±1 (1M)...";
- #~ williams_pp1_ntheory_factor($rem, 1_000_000);
- #~ },
- sub {
- if ($len < 500) {
- say "=> Fibonacci p±1 (500K)...";
- fibonacci_factorization($rem, 500_000);
- }
- },
- sub {
- $state{cyclotomic_check} || return undef;
- say "=> Fast cyclotomic check...";
- my $f = cyclotomic_factorization($rem);
- $f // do { $state{cyclotomic_check} = 0 };
- $f;
- },
- sub {
- $state{fast_power_check} || return undef;
- say "=> Fast power check...";
- my $f = fast_power_check($rem, 500);
- $f // do { $state{fast_power_check} = 0 };
- $f;
- },
- #~ sub {
- #~ say "=> Pollard rho (12M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e12));
- #~ },
- #~ sub {
- #~ say "=> Pollard p-1 (5M)...";
- #~ pollard_pm1_factorial_find_factor($rem, 5_000_000);
- #~ },
- #~ sub {
- #~ say "=> Williams p±1 (3M)...";
- #~ williams_pp1_ntheory_factor($rem, 3_000_000);
- #~ },
- #~ sub {
- #~ say "=> Pollard rho (13M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e13));
- #~ },
- #~ sub {
- #~ say "=> Williams p±1 (5M)...";
- #~ williams_pp1_ntheory_factor($rem, 5_000_000);
- #~ },
- #~ sub {
- #~ if ($len > 40) {
- #~ say "=> Pollard rho (14M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e14));
- #~ }
- #~ },
- #~ sub {
- #~ say "=> Pollard p-1 (8M)...";
- #~ pollard_pm1_ntheory_factor($rem, 8_000_000);
- #~ },
- #~ sub {
- #~ if ($len < 150) {
- #~ say "=> Pollard rho-exp...";
- #~ pollard_rho_exp_find_factor($rem, ($len > 50 ? 2 : 1) * 200);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 50) {
- #~ say "=> Pollard p-1 (10M)...";
- #~ pollard_pm1_factorial_find_factor($rem, 10_000_000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 50) {
- #~ say "=> Williams p±1 (10M)...";
- #~ williams_pp1_ntheory_factor($rem, 10_000_000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard rho (15M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e15));
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard p-1 (20M)...";
- #~ pollard_pm1_factorial_find_factor($rem, 20_000_000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Williams p±1 (20M)...";
- #~ williams_pp1_ntheory_factor($rem, 20_000_000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard rho-exp...";
- #~ pollard_rho_exp_find_factor($rem, 1000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard rho (16M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e16));
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard p-1 (50M)...";
- #~ pollard_pm1_factorial_find_factor($rem, 50_000_000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard p+1 (50M)...";
- #~ williams_pp1_ntheory_factor($rem, 50_000_000);
- #~ }
- #~ },
- #~ sub {
- #~ if ($len > 70) {
- #~ say "=> Pollard rho (16M)...";
- #~ pollard_rho_ntheory_factor($rem, int sqrt(1e16));
- #~ }
- #~ },
- );
- MAIN_LOOP: for (; ;) {
- if ($rem <= 1) {
- last;
- }
- if (is_prime($rem)) {
- push @$factors, $rem;
- $rem = 1;
- last;
- }
- $len = length($rem);
- if ($len <= FACTOR_MIN_LEN) { # factorize with SIQS directly
- return $rem;
- }
- if ( (Math::GMPz::Rmpz_popcount($rem + 1) == 1)
- or (Math::GMPz::Rmpz_popcount($rem - 1) == 1)
- or (Math::GMPz::Rmpz_popcount(($rem * 3) - 1) == 1)
- or (Math::GMPz::Rmpz_popcount(($rem * 3) + 1) == 1)) {
- return $rem;
- }
- printf("\n[*] Factoring %s (%s digits)...\n\n", ($len > MASK_LIMIT ? "C$len" : $rem), $len);
- say "=> Perfect power check...";
- if (defined(my $f = check_perfect_power($rem))) {
- my $exp = 1;
- for (my $t = $f ; $t < $rem ; ++$exp) {
- $t *= $f;
- }
- my @r = (is_prime($f) ? $f : factorize($f));
- push(@$factors, (@r) x $exp);
- return 1;
- }
- foreach my $i (0 .. $#factorization_methods) {
- my $code = $factorization_methods[$i];
- my $f = $code->();
- if (store_factor(\$rem, $f, $factors)) {
- # Move the successful factorization method at the top
- unshift(@factorization_methods, splice(@factorization_methods, $i, 1));
- next MAIN_LOOP;
- }
- }
- last;
- }
- return $rem;
- }
- sub check_perfect_power ($n) {
- # Check whether n is a perfect and return its perfect root.
- # Returns undef otherwise.
- if ((my $exp = is_power($n)) > 1) {
- my $root = Math::GMPz->new(rootint($n, $exp));
- say "`-> perfect power: $root^$exp";
- return $root;
- }
- return undef;
- }
- sub find_prime_factors ($n, $factors) {
- # Return one or more prime factors of the given number n. Assume
- # that n is not a prime and does not have very small factors.
- if (length("$n") > 65) {
- find_all_prime_factors($n, $factors);
- return 1;
- }
- push @$factors, map { Math::GMPz->new("$_") } Math::Prime::Util::GMP::factor($n);
- }
- sub find_all_prime_factors ($n, $factors) {
- if (!ref($n)) {
- $n = Math::GMPz->new($n);
- }
- if (length("$n") < 45) {
- say ":: Factoring with ntheory (<45 len)...";
- push @$factors, map { Math::GMPz->new("$_") } ntheory::factor($n);
- return 1;
- }
- state $t = do {
- my $z = Math::GMPz::Rmpz_init_nobless();
- Math::GMPz::Rmpz_primorial_ui($z, 1e7);
- $z;
- };
- my $g = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_gcd($g, $n, $t);
- my $rem = $n / $g;
- if ($g > $rem or ($g > ~0 and length("$rem") < 80)) {
- say ":: Factoring with ntheory (small factors)...";
- push @$factors, map { Math::GMPz->new("$_") } ntheory::factor($g);
- push @$factors, map { Math::GMPz->new("$_") } ntheory::factor($rem);
- return 1;
- }
- if (length("$n") > 50) {
- say ":: Factoring with FactorDB...";
- my @f;
- if (SCRAPE_FACTORDB) {
- chomp(@f = `$^X /home/swampyx/Other/Programare/experimental-projects/factordb/scrape_factordb.pl $n`);
- }
- else {
- chomp(@f = `$^X /home/swampyx/Other/Programare/experimental-projects/factordb/get_factordb.pl $n`);
- }
- if (Math::GMPz->new(vecprod(@f)) == $n and vecall { is_prime($_) } @f) {
- push @$factors, @f;
- return 1;
- }
- say ":: FactorDB failed...";
- }
- if (length("$n") < 70) {
- say ":: Factoring with ntheory (<70 len)...";
- push @$factors, map { Math::GMPz->new("$_") } Math::Prime::Util::GMP::factor($n);
- return 1;
- }
- say ":: Too large to factor...";
- return 0;
- }
- sub special_form_factorization ($n) {
- my %seen_divisor;
- my @near_power_params;
- my @diff_powers_params;
- my @cong_powers_params;
- my @sophie_params;
- #
- ## Close to a perfect power
- #
- my $near_power = sub ($r, $e, $k) {
- my @factors;
- foreach my $d (ntheory::divisors($e)) {
- my $x = $r**$d;
- foreach my $j (1, -1) {
- my $t = $x - $k * $j;
- my $g = Math::GMPz->new(gcd($t, $n));
- if ($g > TRIAL_DIVISION_LIMIT and $g < $n and !$seen_divisor{$g}++) {
- push @factors, $g;
- }
- }
- }
- @factors;
- };
- foreach my $j (1 .. NEAR_POWER_ITERATIONS) {
- foreach my $k (1, -1) {
- my $u = $k * $j * $j;
- if ($n + $u > 0) {
- if (my $e = is_power($n + $u)) {
- my $r = Math::GMPz->new(rootint($n + $u, $e));
- say "[*] Near power detected: $r^$e ", sprintf("%s %s", ($k == 1) ? ('-', $u) : ('+', -$u));
- push @near_power_params, [$r, $e, $j];
- }
- }
- }
- }
- #
- ## Difference of powers
- #
- my $diff_powers = sub ($r1, $e1, $r2, $e2) {
- my @factors;
- my @d1 = ntheory::divisors($e1);
- my @d2 = ntheory::divisors($e2);
- foreach my $d1 (@d1) {
- my $x = $r1**$d1;
- foreach my $d2 (@d2) {
- my $y = $r2**$d2;
- foreach my $j (1, -1) {
- my $t = $x - $j * $y;
- my $g = Math::GMPz->new(gcd($t, $n));
- if ($g > TRIAL_DIVISION_LIMIT and $g < $n and !$seen_divisor{$g}++) {
- push @factors, $g;
- }
- }
- }
- }
- @factors;
- };
- my $diff_power_check = sub ($r1, $e1) {
- my $u = $r1**$e1;
- my $dx = abs($u - $n);
- if (Math::GMPz::Rmpz_perfect_power_p($dx)) {
- my $e2 = is_power($dx) || 1;
- my $r2 = Math::GMPz->new(rootint($dx, $e2));
- if ($u > $n) {
- say "[*] Difference of powers detected: ", sprintf("%s^%s - %s^%s", $r1, $e1, $r2, $e2);
- }
- else {
- say "[*] Sum of powers detected: ", sprintf("%s^%s + %s^%s", $r1, $e1, $r2, $e2);
- # Sophie Germain's identity:
- # n^4 + 4^(2k+1) = n^4 + 4*(4^(2k)) = n^4 + 4*((2^k)^4)
- if ($r1 == 4 and ($e1 % 2 == 1) and $e2 == 4) { # n = r1^(2k+1) + r2^4
- push @sophie_params, [$r2, Math::GMPz->new(rootint($r1**($e1 - 1), 4))];
- }
- if ($r2 == 4 and ($e2 % 2 == 1) and $e1 == 4) { # n = r2^(2k+1) + r1^4
- push @sophie_params, [$r1, Math::GMPz->new(rootint($r2**($e2 - 1), 4))];
- }
- }
- push @diff_powers_params, [$r1, $e1, $r2, $e2];
- }
- };
- # Sum and difference of powers of the form a^k ± b^k, where a and b are small.
- foreach my $r1 (reverse 2 .. logint($n, 2)) {
- my $t = logint($n, $r1);
- $diff_power_check->(Math::GMPz->new($r1), $t); # sum of powers
- $diff_power_check->(Math::GMPz->new($r1), $t + 1); # difference of powers
- }
- # Sum and difference of powers of the form a^k ± b^k, where a and b are large.
- foreach my $e1 (reverse 2 .. logint($n, 2)) {
- my $t = Math::GMPz->new(rootint($n, $e1));
- $diff_power_check->($t, $e1); # sum of powers
- $diff_power_check->($t + 1, $e1); # difference of powers
- }
- #
- ## Congruence of powers
- #
- my $cong_powers = sub ($r, $e1, $k, $e2) {
- my @factors;
- my @divs1 = ntheory::divisors($e1);
- my @divs2 = ntheory::divisors($e2);
- foreach my $d1 (@divs1) {
- my $x = $r**$d1;
- foreach my $d2 (@divs2) {
- my $y = $k**$d2;
- foreach my $j (-1, 1) {
- my $t = $x - $j * $y;
- my $g = Math::GMPz->new(gcd($t, $n));
- if ($g > TRIAL_DIVISION_LIMIT and $g < $n and !$seen_divisor{$g}++) {
- if ($r == $k) {
- say "[*] Congruence of powers: a^$d1 == b^$d2 (mod n) -> $g";
- }
- else {
- say "[*] Congruence of powers: $r^$d1 == $k^$d2 (mod n) -> $g";
- }
- push @factors, $g;
- }
- }
- }
- }
- @factors;
- };
- my @congrunce_range = reverse(2 .. 64);
- my $process_congruence = sub ($root, $e) {
- for my $j (1, 0) {
- my $k = $root + $j;
- my $u = Math::GMPz::Rmpz_init();
- ref($k)
- ? Math::GMPz::Rmpz_set($u, $k)
- : Math::GMPz::Rmpz_set_ui($u, $k);
- Math::GMPz::Rmpz_powm_ui($u, $u, $e, $n);
- foreach my $z ($u, $n - $u) {
- if (Math::GMPz::Rmpz_perfect_power_p($z)) {
- my $t = is_power($z) || 1;
- my $r1 = rootint($z, $t);
- my $r2 = rootint($z, $e);
- push @cong_powers_params, [Math::GMPz->new($r1), $t, Math::GMPz->new($k), $e];
- push @cong_powers_params, [Math::GMPz->new($r2), $e, Math::GMPz->new($k), $e];
- }
- }
- }
- };
- for my $e (@congrunce_range) {
- my $root = Math::GMPz->new(rootint($n, $e));
- $process_congruence->($root, $e);
- }
- for my $root (@congrunce_range) {
- my $e = Math::GMPz->new(logint($n, $root));
- $process_congruence->($root, $e);
- }
- # Sophie Germain's identity
- # x^4 + 4y^4 = (x^2 + 2xy + 2y^2) * (x^2 - 2xy + 2y^2)
- my $sophie = sub ($A, $B) {
- my @factors;
- foreach my $f (
- #<<<
- $A*$A + (($B*$B)<<1) - (($A*$B<<1)),
- $A*$A + (($B*$B)<<1) + (($A*$B)<<1),
- #>>>
- ) {
- my $g = Math::GMPz->new(gcd($f, $n));
- if ($g > TRIAL_DIVISION_LIMIT and $g < $n and !$seen_divisor{$g}++) {
- push @factors, $g;
- }
- }
- @factors;
- };
- my $sophie_check_root = sub ($r1) {
- {
- my $x = 4 * $r1**4;
- my $dx = $n - $x;
- if (is_power($dx, 4)) {
- my $r2 = Math::GMPz->new(rootint($dx, 4));
- say "[*] Sophie Germain special form detected: $r2^4 + 4*$r1^4";
- push @sophie_params, [$r2, $r1];
- }
- }
- {
- my $y = $r1**4;
- my $dy = $n - $y;
- if (($dy % 4 == 0) and is_power($dy >> 2, 4)) {
- my $r2 = Math::GMPz->new(rootint($dy >> 2, 4));
- say "[*] Sophie Germain special form detected: $r1^4 + 4*$r2^4";
- push @sophie_params, [$r1, $r2];
- }
- }
- };
- # Try to find n = x^4 + 4*y^4, for x or y small.
- foreach my $r1 (map { Math::GMPz->new($_) } 2 .. logint($n, 2)) {
- $sophie_check_root->($r1);
- }
- { # Try to find n = x^4 + 4*y^4 for x and y close to floor(n/5)^(1/4).
- my $k = Math::GMPz->new(rootint($n / 5, 4));
- for my $j (0 .. 1000) {
- $sophie_check_root->($k + $j);
- }
- }
- my @divisors;
- foreach my $args (@near_power_params) {
- push @divisors, $near_power->(@$args);
- }
- foreach my $args (@diff_powers_params) {
- push @divisors, $diff_powers->(@$args);
- }
- foreach my $args (@cong_powers_params) {
- push @divisors, $cong_powers->(@$args);
- }
- foreach my $args (@sophie_params) {
- push @divisors, $sophie->(@$args);
- }
- @divisors = sort { $a <=> $b } @divisors;
- my @factors;
- foreach my $d (@divisors) {
- my $g = Math::GMPz->new(gcd($n, $d));
- if ($g > TRIAL_DIVISION_LIMIT and $g < $n) {
- while ($n % $g == 0) {
- $n /= $g;
- push @factors, $g;
- }
- }
- }
- return sort { $a <=> $b } @factors;
- }
- sub verify_prime_factors ($n, $factors) {
- Math::GMPz->new(vecprod(@$factors)) == $n or die 'product of factors != n';
- foreach my $p (@$factors) {
- is_prime($p) or die "not prime detected: $p";
- }
- sort { $a <=> $b } @$factors;
- }
- sub factorize ($n) {
- # Factorize the given integer n >= 1 into its prime factors.
- if ($n < 1) {
- die "Number needs to be an integer >= 1";
- }
- my $len = length($n);
- printf("\n[*] Factoring %s (%d digits)...\n", $n, $len);
- return () if ($n <= 1);
- return ($n) if is_prime($n);
- if (my $e = is_power($n)) {
- my $root = Math::GMPz->new(rootint($n, $e));
- say "[*] Perfect power detected: $root^$e";
- my @factors = (is_prime($root) ? $root : factorize($root));
- return verify_prime_factors($n, [(@factors) x $e]);
- }
- my @divisors = ((LOOK_FOR_SPECIAL_FORMS and $len > 65) ? special_form_factorization($n) : ());
- if (@divisors) {
- say "[*] Divisors found so far: ", join(', ', @divisors);
- my @composite;
- my @factors;
- foreach my $d (@divisors) {
- $d > 1 or next;
- if (is_prime($d)) {
- push @factors, $d;
- }
- else {
- push @composite, $d;
- }
- }
- push @factors, map { factorize($_) } reverse @composite;
- my $rem = $n / Math::GMPz->new(vecprod(@factors));
- if ($rem > 1) {
- push @factors, factorize($rem);
- }
- return verify_prime_factors($n, \@factors);
- }
- my ($factors, $rem) = ([], $n);
- if (@$factors) {
- say "[*] Prime factors found so far: ", join(', ', @$factors);
- }
- else {
- say "[*] No small factors found...";
- }
- if ($rem != 1) {
- if (LOOK_FOR_SMALL_FACTORS) {
- say "[*] Trying to find more small factors...";
- $rem = find_small_factors($rem, $factors);
- }
- else {
- say "[*] Skipping the search for more small factors...";
- }
- if ($rem > 1) {
- find_all_prime_factors($rem, $factors);
- }
- }
- return verify_prime_factors($n, $factors);
- }
- sub prefactorization ($n, $factors) {
- if (length("$n") <= 45) {
- say ":: Factoring $n with ntheory (<45 len)...";
- push @$factors, Math::Prime::Util::GMP::factor($n);
- return 1;
- }
- state $t = do {
- my $z = Math::GMPz::Rmpz_init_nobless();
- Math::GMPz::Rmpz_primorial_ui($z, 1e7);
- $z;
- };
- my $g = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_gcd($g, $n, $t);
- my $rem = $n / $g;
- if ($g > $rem or ($g > ~0 and length("$rem") < 80) or length("$rem") <= 50) {
- say ":: Factoring $n with ntheory (small factors)...";
- push @$factors, ntheory::factor($g);
- push @$factors, ntheory::factor($rem);
- @$factors = map { Math::GMPz->new("$_") } @$factors;
- @$factors = sort { $a <=> $b } @$factors;
- return 1;
- }
- return 0;
- }
- eval { require GDBM_File };
- my $cache_db = "factors.db";
- dbmopen(my %cache_db, $cache_db, 0666)
- or die "Can't create/access database <<$cache_db>>: $!";
- while (<>) {
- next if /^\h*#/;
- /\S/ or next;
- my $n = (split(' ', $_))[-1];
- $n || next;
- $n =~ /^[0-9]+\z/ or next;
- #~ next if ($n < ~0);
- #~ next if length($n) >= 35;
- next if ($n < ~0);
- next if length($n) > MAX_SIZE;
- #next if length($n) > 1000;
- my $len = length($n);
- next if exists($cache_db{$n});
- #is_pseudoprime($n, 2) && next;
- my @factors;
- if ($len <= FACTOR_MIN_LEN) {
- @factors = Math::Prime::Util::GMP::factor($n);
- }
- else {
- $n = Math::GMPz::Rmpz_init_set_str($n, 10);
- }
- #~ if ( (Math::GMPz::Rmpz_popcount($n + 1) == 1)
- #~ or (Math::GMPz::Rmpz_popcount($n - 1) == 1)
- #~ or (Math::GMPz::Rmpz_popcount(($n * 3) - 1) == 1)
- #~ or (Math::GMPz::Rmpz_popcount(($n * 3) + 1) == 1)) {
- #~ next;
- #~ }
- #next if is_smooth($n, 1e7);
- if (PREFACTORIZATION and not @factors) {
- prefactorization($n, \@factors);
- }
- if (not @factors) {
- if (FACTORIZE_WITH_SIDEF) {
- @factors = map { Math::Sidef::factor($_) } Math::Sidef::special_factors($n, 0);
- }
- else {
- @factors = eval { factorize($n) }
- }
- }
- if (scalar(@factors) == 1) {
- Math::Prime::Util::GMP::is_provable_prime($n) || die "BPSW counter-example(?): $n";
- }
- if ( scalar(@factors) > 1
- and (vecall { is_prime($_) } @factors)
- and Math::GMPz->new(vecprod(@factors)) == Math::GMPz->new($n)) {
- say ":: Successfully factorized!";
- $cache_db{"$n"} = join(' ', @factors);
- }
- }
- dbmclose(%cache_db);
|