PARI's commands

[3.14159]

Wait; How do you code Pollard's Rho?

[science_man_88]

Originally Posted by 3.14159

Wait; How do you code Pollard's Rho?

http://en.wikipedia.org/wiki/Pollard's_rho_algorithm may help.

[science_man_88]

If I knew a few more things maybe I could make it work,However I don't so I won't.

[CRGreathouse]

Originally Posted by Batalov

Thanks! That's what I was looking for.
(and yes, subst() works, too, but I was looking for a hint to best practices)

Yes, thats certainly substpol. Honestly, I haven't tested this vs. subst for performance or anything, but I imagine there are reasons to split the special case out.

[Batalov]

Here's the script by D. Broadhurst in all its beauty minus some lines that I couldn't understand or that were not entirely by the book and cost nothing to be by the book. The script's result appears to tell us that this
27810-digit number is prime by Konyagin-Pomerance N-1 test with 30.2% factored N-1.

I was looking for a way not to put the 7334-digit prime there literally.

If this is correct by the book, I could proceed adding it to PFGW. Does it look proper to you, Charles? Thanks.

[Batalov]

P.S. That previous number was actually 32.99% factored (but still short of B-L-S N-1 proof; N+1 side doesn't help enough). Here's the 30.2% factored N-1-based KP proof. Runs really fast, because the number is fairly small and within reach of Primo anyway.

[CRGreathouse]

Originally Posted by Batalov

If this is correct by the book, I could proceed adding it to PFGW. Does it look proper to you, Charles? Thanks.

It looks fine. There are things I would do differently, to wit:

test()={
  print("D Broadhurst, KP N-1 test, 8 Jul 2001");
  \\ checked for correctness as a Konyagin-Pomerance N-1 test

  \\ Origin: http://tech.groups.yahoo.com/group/primeform/files/KP/KonPom.gp
  \\ KP Combined test seemed to be a work in progress
  \\ it would need more operations as were in the original script 
  \\ (See Theorem 4.2.9, p.186 in Crandall-Pomerance PN-ACP)

  my(N=(8*10^27811-71)/9,p7334,p557,p202,lsm,F,R,c1,c4);

  \\ Supply N-1 factors in the lsm[] array:

  p7334 = subst(polcyclo(13905),x,10)/83431/333721;
  p557  = subst(polcyclo(927),x,10)/46351/181693/1405333/43604227/75178674467317/664958944842271489;
  p202  = subst(polcyclo(618),x,10)/619;

  lsm=[2,2,2,2,3,3,3,5,7,11,13,19,31,37,41,211,241,271,619,1031,1237,2161,7211,9091,29611,46351,52579,181693,238681,333667,380071,\
497491,704521,1358983,1405333,2906161,3762091,7034077,\
43604227,44092859,102860539,569836171,984385009,2013681931,5905014721,326345481191,8595591045751,8985695684401,75178674467317,\
301917502615411,612053256358933,2702394989404991,41343860637475219,167798066344417387,664958944842271489,4185502830133110721,\
244261115576975427841,290934235666553153131,3462948230491376473027,182725114866521155647161,567664121498896654896380203841,\
896048585318577702680084550566846611,1471865453993855302660887614137521979,7073532436575439739956079653594194521,\
85517237338860102283865323577960208481,5294796903161592416528456780680376286484870226446771978908657527791,\
153211620887015423991278431667808361439217294295901387715486473457925534859044796980526236853,\
292815973319957863643663090772618219691707118842359217315010640397456279188940711108958040881306849,\
757,5563,6481,83431,333721,1577071,16357951,22053331,70541929,310362841,2016447481,\
14175966169,18245696041,258360989311,577603663291,31023833790241,440334654777631,\
483418418597220677238517353915231961831,8610583349234340055547908764091017276717091,\
p202,p557,p7334];

  F=prod(k=1,length(lsm),lsm[k]);
  R=(N-1)/F;

  if(denominator(R)==1
  && gcd(F,R)==1
  && F^(1/3)/6>1
  && 1>3*N/F^(10/3),
  ,
    error("Fail")
  );

  \\ Here, log(F)/log(N) =~ 0.3299 > 3/10, but less than 1/3
  print("Factored part is: ", log(F)/log(N));

  c1=divrem(R,F);
  c4=c1[1];
  c1=c1[2];

  print();

  \\ Unsure why the square test was not performed six times in the original script
  \\ This is part 1 from Theorem 4.1.6 (or 4.2.9):
  for(t=0,5,
    if(issquare((c1+t*F)^2+4*t-4*c4),error("Fail @ "t))
  );

  my(b=contfrac(c1/F), v=0, vn=1, u=1, un=b[1], n=1);
  while(F^(1/3)>vn,
    bn=b[n++];
    uo=u;
    u=un;
    vo=v;
    v=vn;
    un=bn*un+uo;
    vn=bn*vn+vo
  );

  \\ Here, u/v is as per Part 2 of Theorem 4.1.6 

  my(d=round(c4*v/F), zs=[v,u*F-c1*v,c4*v-d*F+u,-d], q=0);
  for(k=1,4,
    print("zs["k"] = ", zs[k]);
    q += zs[k]*x^(4-k)
  );

  print();

  \\ required precision is dependent on the size of N
  default(realprecision, 10000);
  my(ps=polroots(q), r, q1, q2);

  for(k=1,3,
    r=floor(real(ps[k]));
    print(r);
    q1=substpol(q, x, r);
    q2=substpol(q, x, r+1);
    if(q1*q2>=0,error("Fail root ",k))
  );

  print("OK, no integral roots");
  \\ if there's a root a, check: if aF+1 is a nontrivial factor of N, then N is composite

  \\ Note to myself: for implementation in PFGW: one root is very close to 0 and can be found by Newton iter, with initial 0; then reduce to quadratic.
  print(ps[2]);
};

test()

but these are largely stylistic differences.

[3.14159]

Originally Posted by CRGreathouse

That's actually usual for Fermat's method, so this doesn't indicate that your version is bad.

Meh. I think I've read too many kook sites claiming to have some drastic improvement to Fermat's algorithm in one way or another, and wished to entertain the thought of programming those in, then trying to see if the algorithm, with these "improvements, actually competes with trial factoring, for a fairly large semiprime number.

Here's one of them.

Some sort of iterated algorithm to try, for some fairly large semiprime n (Ex: 5403615307171909)

.. Oh, dear. This tutorial of theirs is a nightmare.

[science_man_88]

Originally Posted by 3.14159

Meh. I think I've read too many kook sites claiming to have some drastic improvement to Fermat's algorithm in one way or another, and wished to entertain the thought of programming those in, then trying to see if the algorithm, with these "improvements, actually competes with trial factoring, for a fairly large semiprime number.

Here's one of them.

Some sort of iterated algorithm to try, for some fairly large semiprime n (Ex: 5403615307171909)

.. Oh, dear. This tutorial of theirs is a nightmare.

the main part of the site I don't get without the code is the table they list off. I should know vb.net a bit though it's only been under 2 years.

[science_man_88]

Figured it out I think. Want my take on it Pi ?

[3.14159]

Originally Posted by science_man_88

Figured it out I think. Want my take on it Pi ?

Sure.