Division free ECM factoring

[Bill Hart]

Hi all,

I just got done writing a Pari script for factoring integers with the
ECM method without using modular inverses.

It's not really something I can publish (I'm sure it's a well-known
trick). So I decided to post it here for fun.

Clearly it can be converted to a Python script (which I leave as a
trivial exercise), if you don't want to use Pari/GP mode (sage -gp).

It doesn't have a stage 2, so it can only find up to about 15 digit
factors easily.

Note I also start it by giving it 300MB of memory. Also note it may
have bugs!

allocatemem(300000000)

mul8(L,k,m) = {

if(k==1, return(L));

Z=mul8(L, floor(k/2), m);

if (Z[1] == "done", return(Z));

s=Z[1]*Z[2]*Z[3]*Z[4]*Z[5]*Z[6]*Z[7]*Z[8];

if(k%2 == 0,
if(s != 0,
z = (Z[2]*Z[4]*Z[5]^2-Z[4]*Z[6]*Z[3]^2)%m;
s = gcd(z,m);
if (s != 1, return(["done",s,0,0,0,0,0,0])));

return([(Z[4]*Z[2]^3-Z[3]^3*Z[1])%m,
(Z[3]*(Z[5]*Z[2]^2-Z[1]*Z[4]^2))%m,
(Z[5]*Z[3]^3-Z[4]^3*Z[2])%m,
(Z[4]*(Z[6]*Z[3]^2-Z[2]*Z[5]^2))%m,
(Z[6]*Z[4]^3-Z[5]^3*Z[3])%m,
(Z[5]*(Z[7]*Z[4]^2-Z[3]*Z[6]^2))%m,
(Z[7]*Z[5]^3-Z[6]^3*Z[4])%m,
(Z[6]*(Z[8]*Z[5]^2-Z[4]*Z[7]^2))%m]),

if(s != 0,
z = (Z[2]*Z[4]*Z[6]^2-Z[5]*Z[7]*Z[3]^2)%m;
s = gcd(z,m);
if(s != 1, return(["done",s,0,0,0,0,0,0])));

return([(Z[3]*(Z[5]*Z[2]^2-Z[1]*Z[4]^2))%m,
(Z[5]*Z[3]^3-Z[4]^3*Z[2])%m,
(Z[4]*(Z[6]*Z[3]^2-Z[2]*Z[5]^2))%m,
(Z[6]*Z[4]^3-Z[5]^3*Z[3])%m,
(Z[5]*(Z[7]*Z[4]^2-Z[3]*Z[6]^2))%m,
(Z[7]*Z[5]^3-Z[6]^3*Z[4])%m,
(Z[6]*(Z[8]*Z[5]^2-Z[4]*Z[7]^2))%m,
(Z[8]*Z[6]^3-Z[7]^3*Z[5])%m])
)
}

mul(L,n,m) = {
if(n < 2^1024,
return(mul8(L, n, m)),

M = mul(L, floor(n/(2^1024)), m);
return(mul8(M, 2^1024+n%(2^1024), m))
)
}

fac(m) = {
num=1000;
j=0; i=1; a=2^8*3^6*5^4*7^2; p1=2;
while (i < num, i = i+1; a = a*p1; p1 = nextprime(p1+1));

while(1,k=random(m);
L=[m-1,m-1,0,1,1,1,k,(k-1)%m];
L=mul(L,a,m);

if (L[1] == "done", return(L[2]));

j=j+1;
if (j%50 == 0,
for (k = num, 3*num,
a = a*p1;
p1 = nextprime(p1)
);
num = num*3;
j=0
)
)
}

As an example of its use:

fac(nextprime(random(10^15))*nextprime(random(10^16)))

Enjoy!

Bill.