#############################################################################
##
#W  cfrac.gi               GAP4 Package `FactInt'                 Stefan Kohl
##
##  This file contains functions for factorization using the
##  Continued Fraction Algorithm (Brillhard-Morrison Algorithm).
##
##  Argument of FactorsCFRAC:
##
##  <n>        the integer to be factored
##
##  The result is returned as a list of the prime factors of <n>.
##
#############################################################################

BindGlobal("CFRACSplit", function (n)

  local Digits,
        Multiplier,FactorBase,FactorBaseSize,MaxFactorBaseEl,pi,sqrt,
        Ai,Bi,Ci,Pi,Aj,Bj,Cj,Pj,Ak,Bk,Ck,Pk,Step,CiBuf,CiFact,
        abort1,abort2,AbortingTime1,AbortingTime2,LargePrimeLimit,
        CollectingInterval,CollectingIntervals,NextCollectionAt,
        Factored,FactoredLarge,LargeFactors,UsableFactors,
        FactoredLargeUsable,RelationsLarge,RelPos,Factorizations,
        RelsFB,RelsLarge,RelsTotal,Progress,Remaining,Required,
        LeftMatrix,RightMatrix,CanonFact,Line,Col,Lines,Cols,
        X,Y,YQuadFactors,FactorsFound,p,DependencyNr,Ready,Result,
        Pair1,Pair2,Fact,q,pos,zero,one,i,j,StartingTime,UsedTime;

  StartingTime := Runtime();

  Digits := LogInt(n,10) + 1;

  Info(IntegerFactorizationInfo,2,"CFRAC for n = ",n);
  Info(IntegerFactorizationInfo,2,"Digits              : ",Digits);

  zero := Zero(GF(2)); one := One(GF(2));
  CollectingIntervals := [1000,2000,5000,10000,20000,50000,100000,200000,
                          500000,1000000];
  CollectingInterval := CollectingIntervals[Minimum(Int(Digits/5),
                                            Length(CollectingIntervals))];
  NextCollectionAt := CollectingInterval;

  # Generate the factor base

  Multiplier := n mod 8; n := n * Multiplier;
  FactorBaseSize := QuoInt(Digits^3,200);
  FactorBase := [-1]; pi := 1;
  for i in [2..FactorBaseSize] do
    repeat
      pi := NextPrimeInt(pi);
    until Legendre(n,pi) = 1;
    FactorBase[i] := pi;
  od;
  MaxFactorBaseEl := FactorBase[FactorBaseSize];
  LargePrimeLimit := MaxFactorBaseEl^2;
  Required        := FactorBaseSize + 20;

  Info(IntegerFactorizationInfo,2,"Size of factor base : ",
                                   FactorBaseSize);
  Info(IntegerFactorizationInfo,3,"\nFactor base : \n",FactorBase,"\n");

  # Some initializations

  Factored := [];
  LargeFactors := []; FactoredLarge := [];
  FactoredLargeUsable := []; RelationsLarge := [];

  sqrt := RootInt(n);

  Aj   := sqrt;
  Bk   := 0; Bj := sqrt;
  Ck   := 1; Cj := n - sqrt^2;
  Pk   := 1; Pj := sqrt;

  Step := 1;

  # Abort Trial Division after dividing Ci by the first <AbortingTime> 
  # primes in the factor base if the the unfactored part after that
  # is larger than <abort>

  abort1        := RootInt(n,12)^5;
  abort2        := RootInt(n,3); 
  AbortingTime1 := Maximum(10,QuoInt(FactorBaseSize,20));
  AbortingTime2 := QuoInt(FactorBaseSize,5);

  # Calculate the continued fraction expansion of the square root of n
  # until there are enough factored Ci's
 
  repeat
    
    # Generate the next set of values of the considered five sequences,
    # especially the next Ci

    Step := Step + 1;
    Ai := QuoInt(sqrt + Bj,Cj);   
    Bi := Ai * Cj - Bj;
    Ci := Ck + Ai * (Bj - Bi);
    Pi := (Pk + Ai * Pj) mod n;

    # Trial divide Ci

    if Step > 100 then
      CiBuf := Ci; 
      if Step mod 2 = 0 then CiFact := [];
                        else CiFact := [-1]; 
      fi;
      for i in [2..AbortingTime1] do
        while CiBuf mod FactorBase[i] = 0 do
          CiBuf := CiBuf/FactorBase[i];
          Add(CiFact,FactorBase[i]);
        od;
      od;
      if CiBuf <= abort1 then
        for i in [AbortingTime1 + 1..AbortingTime2] do
          while CiBuf mod FactorBase[i] = 0 do
            CiBuf := CiBuf/FactorBase[i];
            Add(CiFact,FactorBase[i]);
          od;
        od;
        if CiBuf <= abort2 then
          for i in [AbortingTime2 + 1..FactorBaseSize] do
            while CiBuf mod FactorBase[i] = 0 do
              CiBuf := CiBuf/FactorBase[i];
              Add(CiFact,FactorBase[i]);
            od;
          od;
          if CiBuf < LargePrimeLimit then
            if CiBuf = 1 then
              Add(Factored,[Pi,CiFact]);
            else
              Add(LargeFactors,CiBuf);
              Add(CiFact,CiBuf);
              Add(FactoredLarge,[Pi,CiFact]);
            fi;
          fi;   
        fi;
      fi;
    fi;
    
    # Look for usable factorizations with a large factor

    if Step >= NextCollectionAt
    then

      Info(IntegerFactorizationInfo,3,"");
      Info(IntegerFactorizationInfo,3,
           "Collecting relations with a large factor");

      Sort(LargeFactors);
      UsableFactors := Set(List(Filtered(Collected(LargeFactors),
                                Pair1->Pair1[2] > 1),Pair2->Pair2[1]));
      FactoredLargeUsable := 
      Filtered(FactoredLarge,
               Fact->ForAll(Fact[2],q ->   q <= MaxFactorBaseEl 
                                        or q in UsableFactors));
      Sort(FactoredLargeUsable,
           function(f1,f2)
             return   Intersection(f1[2],UsableFactors)[1]
                    < Intersection(f2[2],UsableFactors)[1];
           end);
      RelationsLarge := []; pos := 1; RelPos := 1;
      while pos < Length(FactoredLargeUsable) do
        if   Intersection(FactoredLargeUsable[pos    ][2],UsableFactors)
           = Intersection(FactoredLargeUsable[pos + 1][2],UsableFactors)
        then
          RelationsLarge[RelPos] := 
          [FactoredLargeUsable[pos][1] * FactoredLargeUsable[pos + 1][1],
           Concatenation(FactoredLargeUsable[pos    ][2],
                         FactoredLargeUsable[pos + 1][2])];
          RelPos := RelPos + 1;
        fi;
        pos := pos + 1;
      od;

      RelsFB    := Length(Factored);
      RelsLarge := Length(RelationsLarge);
      RelsTotal := RelsFB + RelsLarge;
      Remaining := Maximum(0,Required - RelsTotal);
      UsedTime  := Int((Runtime() - StartingTime)/1000);
      Progress  := Minimum(100,Int(100 * RelsTotal/Required));

      if InfoLevel(IntegerFactorizationInfo) < 3 then
        Info(IntegerFactorizationInfo,2,
             "Step ",String(Step,9)," : Factored/FB.: ",String(RelsFB,4),
              ", w. Large Fact.: ",String(RelsLarge,4),
              ", Progress :",String(Progress,3),"%");
      fi;  

      Info(IntegerFactorizationInfo,3,
           "Steps                                          : ",
           String(Step,10));
      Info(IntegerFactorizationInfo,3,
           "Complete factorizations over the factor base   : ",
           String(RelsFB,10));
      Info(IntegerFactorizationInfo,3,
           "Total factorizations with a large prime factor : ",
           String(Length(FactoredLarge),10));
      Info(IntegerFactorizationInfo,3,
           "Relations with a large prime factor            : ",
           String(RelsLarge,10));
      Info(IntegerFactorizationInfo,3,
           "Relations remaining to be found                : ",
           String(Remaining,10));
      Info(IntegerFactorizationInfo,3,
           "Elapsed runtime                                : ",
           String(UsedTime,10)," sec.");
      Info(IntegerFactorizationInfo,3,
           "Progress (relations)                           : ",
           String(Progress,10)," %");
      Info(IntegerFactorizationInfo,3,"");

      NextCollectionAt := NextCollectionAt + CollectingInterval;

    fi;

    Ak := Aj; Bk := Bj; Ck := Cj; Pk := Pj;
    Aj := Ai; Bj := Bi; Cj := Ci; Pj := Pi;   

  until Length(Factored) + Length(RelationsLarge) >= Required;

  # Create exponent matrix

  Info(IntegerFactorizationInfo,2,"Creating the exponent matrix");
  LeftMatrix := []; Cols := FactorBaseSize;
  Factorizations := Concatenation(Factored,RelationsLarge);
  Lines := Length(Factorizations);
  for Line in [1..Lines] do
    LeftMatrix [Line] := ListWithIdenticalEntries(Cols,zero);
    CanonFact := Collected(Factorizations[Line][2]);
    for i in [1..Length(CanonFact)] do
      if   CanonFact[i][2] mod 2 = 1 
      then LeftMatrix[Line]
                     [Position(FactorBase,CanonFact[i][1])] := one;
      fi;
    od;
  od;

  # Do Gaussian Elimination

  Info(IntegerFactorizationInfo,2,
       "Doing Gaussian Elimination, #rows = ",Lines,
       ", #columns = ",Cols);
  RightMatrix := NullspaceMat(LeftMatrix);

  # Calculate X and Y such that X^2 = Y^2 mod n 
  # and check if 1 < gcd(X - Y,n) < n

  Info(IntegerFactorizationInfo,2,"Processing the zero rows"); 
  p := 1; FactorsFound := []; DependencyNr := 1; Line := 1; Ready := false;

  while Line <= Length(RightMatrix) and not Ready do 

    X := 1; Y := 1;
    for Col in [1..Lines] do 
      if   RightMatrix[Line][Col] = one
      then X := X * Factorizations[Col][1] mod n; fi;
    od;
    YQuadFactors := 
    Collected(Concatenation(List(Filtered([1..Lines],
                                          i->RightMatrix[Line][i] = one),
                                 j->Factorizations[j][2])));
    for i in [1..Length(YQuadFactors)] do
      Y := Y * YQuadFactors[i][1]^(YQuadFactors[i][2]/2) mod n;
    od;
    if  (X^2 - Y^2) mod n <> 0 
    then Error("Internal Error : X^2 - Y^2 mod n <> 0"); fi;

    p := Gcd(X - Y,n/Multiplier);

    if not p in [1,n/Multiplier] 
    then 
      Info(IntegerFactorizationInfo,2,
           "Dependency no. ",DependencyNr," yielded factor ",p); 
      Add(FactorsFound,p); FactorsFound := Set(FactorsFound);
      if FactorsTD(n/Multiplier,FactorsFound)[2] = [] 
      then Ready := true; fi;
    else
      Info(IntegerFactorizationInfo,2,
           "Dependency no. ",DependencyNr," yielded no factor"); 
    fi;

    Line         := Line + 1; 
    DependencyNr := DependencyNr + 1;

  od;
  
  if   FactorsFound = []
  then Error("\nSorry, CFRAC has failed ...\n",
             "Perhaps the continued fraction expansion of the square\n",
             "root of the number you attempted to factor has properties\n",
             "disadvantageous for obtaining a factorization\n",
             "with this algorithm, try using the MPQS\n\n"); return [n]; fi;
      
  Result := Flat(FactorsTD(n/Multiplier,FactorsFound));

  Info(IntegerFactorizationInfo,1,"The factors are\n",Result);
  Info(IntegerFactorizationInfo,2,"Digit partition : ",
       List(Result,p -> LogInt(p,10) + 1));
  UsedTime := Runtime() - StartingTime;
  Info(IntegerFactorizationInfo,2,
       "CFRAC runtime : ",TimeToString(UsedTime),"\n");

  return Result;
end);

#############################################################################
##
#F  FactorsCFRAC( <n>, [ <ContinueFromFile>, <PagingDir> ] )
##
InstallGlobalFunction( FactorsCFRAC,

function ( n )

  local  FactorsList, StandardFactorsList, m,
         Ready, Pos, Passno;

  if   not (IsInt(n) and n > 1)
  then Error("Usage : FactorsCFRAC( <n> ), ",
             "where <n> has to be an integer > 1"); fi;

  if   ValueOption("NoPreprocessing") = true 
  then FactorsList := Flat(FactorsTD(n));
  else Info(IntegerFactorizationInfo,2,
            "Doing some preprocessing using Pollard's Rho");
       FactorsList := Flat(FactorsRho(n,1,16,8192)); 
  fi;

  Passno := 0;
  repeat
    Passno := Passno + 1;
    Info(IntegerFactorizationInfo,3,"Pass no. ",Passno);

    StandardFactorsList := [[],[]];
    for m in FactorsList do 
      if IsProbablyPrimeInt(m) then Add(StandardFactorsList[1],m);
                       else Add(StandardFactorsList[2],m); fi;
    od;
    ApplyFactoringMethod(FactorsPowerCheck,[FactorsCFRAC,"FactorsCFRAC"],
                         StandardFactorsList,infinity);
    FactorsList := Flat(StandardFactorsList);

    for Pos in [1..Length(FactorsList)] do 
      if   not IsProbablyPrimeInt(FactorsList[Pos]) 
      then FactorsList[Pos] := CFRACSplit(FactorsList[Pos]); fi;
    od;
    FactorsList := Flat(FactorsList);
    Ready := ForAll(FactorsList,IsProbablyPrimeInt);
  until Ready;

  Sort(FactorsList);
  FactorizationCheck(n,FactorsList);
  return FactorsList;
end);

#############################################################################
##
#E  cfrac.gi . . . . . . . . . . . . . . . . . . . . . . . . . . .  ends here