(* function that finds a non-trivial factor of n, or prove n prime; this is experimential code *)
factorTest[n_]:=
(
(* initialize variables *)
n0=Floor[Sqrt[n]];ncbrt=(17n)^(1./3);
x0=Min[Ceiling[50ncbrt],n0]; x=x0+2;HL=Floor[x/ncbrt];

(* trial division *)
k=2;While[1<k≤x0,{If[Mod[n/k,1]==0,{Print[k,"   ",n/k];Break[];}];k++;}];

(* loop through blocks of integers [x-H,x+H] to search for a factor of n, exit once a factor is found *)
counter=0;
While[x-HL≤n0,{
(* find a suitable continued fraction convergnet to approximate n/x^2 *)
convgnts=Convergents[n/x^2];r=Length[convgnts];
j=0;While[j+1≤ r &&Denominator[convgnts[[j+1]]]≤4HL ,j++];

(* construct the polynomial g_{n,x} *)
bq=convgnts[[j]];b=Numerator[bq];q=Denominator[bq];
qnx=q n/x;qnxx=qnx/x;a=Round[qnx];eps1=qnx-a;eps2q'=qnxx-b;
c0=x eps1;c1=eps1-x eps2q';c2=qnxx-eps2q';

(* solve the equation g_{n,x}(h)=0 and for each solution h test if x+h divides n *)
If[c2≠0,{sqrtdisc=Sqrt[c1^2-4c0 c2]/(2c2);hmx=-c1/(2c2);
h1 = hmx+sqrtdisc;h2=hmx-sqrtdisc;
If[Mod[h1,1]==0 && Mod[n/(x+h1),1]==0, {Print[x+h1,"   ",n/(x+h1)];Break[];}];
If[Mod[h2,1]==0 &&Mod[n/(x+h2),1]==0, {Print[x+h2,"   ",n/(x+h2)];Break[];}];}];
If[c2==0 && c1≠0,{h=-c0/c1;
If[Mod[h,1]==0 &&Mod[n/(x+h),1]==0,{Print[h,"   ",n/(x+h)];Break[];}]}];


(* increment x *)
HR1=Floor[0.4(1-Abs[eps1])/Abs[eps2q']];HR2=Floor[Sqrt[0.6(1-Abs[eps1])x/qnxx]];HR=Min[HR1,HR2];

x=x+HR+HL+1;HL=Floor[x/ncbrt];counter++;
}];

(* print out the trivial factors *)
Print[1,"   ",n,"   ",HL,"   ",HR,"    ",counter];
)

nn=Prime[20000000]Prime[20000001];
sqnn=Floor[Sqrt[nn]];
Timing[factorTest[nn]]