|
#############################################################################
##
#W pplus1.gi GAP4 Package `FactInt' Stefan Kohl
##
## This file contains functions for factorization using a variant of
## Williams' p+1.
##
## Arguments of FactorsPplus1:
##
## <n> the integer to be factored
## <Residues> the number of residues that should be examined
## (the probability of hitting a usable one is
## approx. 1 - 1/2^<Residues>)
## <Limit1> the limit for the first stage
## <Limit2> the limit for the second stage
##
## The result is returned as a list of two lists. The first list
## contains the prime factors found, and the second list contains
## remaining unfactored parts of <n>, if there are any.
##
#############################################################################
BindGlobal("Pplus1Product", function (a,b,n)
local c,c12,a12b12;
a12b12 := a[1][2] * b[1][2];
c12 := (a[1][1] * b[1][2] + a[1][2] * b[2][2]) mod n;
c := [[(a[1][1] * b[1][1] + a12b12) mod n, c12],
[c12, (a12b12 + a[2][2] * b[2][2]) mod n]];
return c;
end);
BindGlobal("Pplus1Square", function (a,n)
local c,c12,a12a12;
a12a12 := a[1][2] * a[1][2];
c12 := a[1][2] * (a[1][1] + a[2][2]) mod n;
c := [[(a[1][1]^2 + a12a12) mod n, c12],
[c12, (a12a12 + a[2][2]^2) mod n]];
return c;
end);
BindGlobal("Pplus1Power", function (Base,exp,n)
local Power,BinExp,i;
BinExp := CoefficientsQadic(exp,2);
Power := Base;
for i in [Length(BinExp),Length(BinExp)-1..2] do
Power := Pplus1Square(Power,n);
if BinExp[i-1] = 1 then
Power := Pplus1Product(Power,Base,n);
fi;
od;
return Power;
end);
BindGlobal("Pplus1Split", function (n,Residues,Limit1,Limit2)
local Residue,ResNo,a,
PowerOfa,PowerAfterFirstStage,p,pExponent,
DiffPowers,DiffPos,DiffSum,NextDiff,DiffsLng,BufProd,
FactorFoundAndReady,FactorsFound;
FactorFoundAndReady := function (PowerProd)
local Result,i,j,k;
Result := Gcd(PowerProd,n);
if not Result in [1,n] then
Info(IntegerFactorizationInfo,1,LogInt(Result,10) + 1,
"-digit factor ",Result," was found");
Add(FactorsFound,Result); n := n/Result;
if IsProbablyPrimeInt(n) then Add(FactorsFound,n); return true; fi;
if IsBound(DiffPowers) then
for i in [1..Length(DiffPowers)] do
if IsBound(DiffPowers[i]) then
for j in [1,2] do for k in [1,2] do
DiffPowers[i][j][k] := DiffPowers[i][j][k] mod n;
od; od;
fi;
od;
fi;
fi;
return false;
end;
if IsProbablyPrimeInt(n) then return [n]; fi;
FactorsFound := [];
Info(IntegerFactorizationInfo,2,
"p+1 for n = ",n,"\nResidues : ",Residues,
", Limit1 : ",Limit1,", Limit2 : ",Limit2);
Residue := 1;
for ResNo in [1..Residues] do
Info(IntegerFactorizationInfo,2,"Residue no. ",ResNo);
a := [[Residue,1],[1,0]];
Info(IntegerFactorizationInfo,3,"First stage");
p := 2; PowerOfa := a;
while p <= Limit1 do
pExponent := LogInt(Limit1,p);
PowerOfa := Pplus1Power(PowerOfa,p^pExponent,n);
if FactorFoundAndReady(PowerOfa[1][2])
then return FactorsFound; fi;
p := NextPrimeInt(p);
od;
Info(IntegerFactorizationInfo,3,"Second stage");
PowerAfterFirstStage := PowerOfa; PowerOfa := [[1,0],[0,1]];
DiffPowers := []; DiffPos := 1; DiffSum := 0;
DiffsLng := Length(PrimeDiffs); BufProd := 1;
while (DiffSum <= Limit2) and (DiffPos <= DiffsLng) do
NextDiff := PrimeDiffs[DiffPos];
DiffSum := DiffSum + NextDiff;
if not IsBound(DiffPowers[NextDiff])
then DiffPowers[NextDiff] :=
Pplus1Power(PowerAfterFirstStage,NextDiff,n); fi;
PowerOfa := Pplus1Product(PowerOfa,DiffPowers[NextDiff],n);
BufProd := BufProd * PowerOfa[1][2] mod n;
if DiffPos mod 50 = 0 then
if FactorFoundAndReady(BufProd)
then return FactorsFound; fi;
fi;
DiffPos := DiffPos + 1;
od;
Residue := NextPrimeInt(Residue);
od;
Add(FactorsFound,n); return FactorsFound;
end);
#############################################################################
##
#F FactorsPplus1( <n>, [ [ <Residues> ], <Limit1>, [ <Limit2> ] ] )
##
InstallGlobalFunction( FactorsPplus1,
function ( arg )
local n, Residues, Limit1, Limit2, GetArg, ArgCorrect,
FactorsList, m, Split, q;
GetArg := function (ArgPos,ArgName,ArgDefault)
if IsBound(arg[ArgPos]) then
if arg[ArgPos] <> fail then return arg[ArgPos];
else return ArgDefault; fi;
fi;
if ValueOption(ArgName) <> fail then return ValueOption(ArgName); fi;
return ArgDefault;
end;
ArgCorrect := Length(arg) in [1..4];
if ArgCorrect then
n := arg[1];
if Length(arg) = 4
then
Residues := arg[2];
Limit1 := GetArg(3,"Pplus1Limit1",Int(PrimeDiffLimit/40));
Limit2 := GetArg(4,"Pplus1Limit2",40 * Limit1);
else
Residues := GetArg(5,"Pplus1Residues",2);
Limit1 := GetArg(2,"Pplus1Limit1",Int(PrimeDiffLimit/40));
Limit2 := GetArg(3,"Pplus1Limit2",40 * Limit1);
fi;
if not (IsInt(n) and n >= 1 and IsInt(Residues) and Residues >= 1
and IsPosInt(Limit1) and IsPosInt(Limit2))
then ArgCorrect := false; fi;
fi;
if not ArgCorrect
then Error("Usage : FactorsPplus1( <n>, [ [ <Residues> ], <Limit1>, ",
"[ <Limit2> ] ] ), where <n>, <Residues>, ",
"<Limit1> and <Limit2> have to be positive integers");
fi;
if IsProbablyPrimeInt(n) then return [[n],[]]; fi;
InitPrimeDiffs(Limit2);
FactorsList := FactorsTD(n);
if FactorsList[2] <> [] then
m := FactorsList[2][1];
if SmallestRootInt(m) < m
then ApplyFactoringMethod(FactorsPowerCheck,
[FactorsPplus1,"FactorsPplus1"],
FactorsList,infinity);
else
Split := Pplus1Split(m,Residues,Limit1,Limit2);
FactorsList[2] := [];
for q in Split do if IsProbablyPrimeInt(q)
then Add(FactorsList[1],q);
else Add(FactorsList[2],q); fi; od;
fi;
fi;
Sort(FactorsList[1]); Sort(FactorsList[2]);
FactorizationCheck(n,FactorsList);
return FactorsList;
end);
#############################################################################
##
#E pplus1.gi . . . . . . . . . . . . . . . . . . . . . . . . . . ends here
[ Dauer der Verarbeitung: 0.36 Sekunden
(vorverarbeitet)
]
|
