A faster way to factor n
This is an alternative to the solution posted by the_one_and_only_az.
We can write $q = (1+I) p + r$ with $r$ a rather small gaussian integer. We can use a version of the classical Fermat factoring algorithm as follows, inspired by a paper of de Weger [AAECC, 13, 2002, p. 17-29].
Write $$4 (1+I) n = x^2 - y^2$$ or $$x^2 - 4 (1+I)n = y^2.$$ We iterate on $x$ until the difference is the square of a gaussian number. The value of $y$ will be small, since $|q - (1+I)p|$ is (this is the analogous of $\Delta$ in the above article).
We start from a rounded version of a complex approximation to $2 \sqrt{(1+I) n}$ and we increment $x$ in any direction until we hit a square in $Z[I]$.
Checking that $z$ is square in $Z[I]$ requires that the norm of $z$ is a square in $Z$ and there are a lot of classical tricks to speed up the process like testing quadratic residuocity modulo small primes.
In our case, the search finished rather quickly. I add a Maple code for this.
The first part finds a best approximation to the squareroot and then
we look for ds and dt until we hit a square. After that, we recover
$p$ and $q$ using a rational approach to compute the squareroot in $Z[I]$
using some ancilary functions.
Code
Inorm := proc(x) evalc(x * subs(I = -I, x)); end:
Gdiv2 := proc(x)
if coeff(x, I, 0) mod 2 = 1 or coeff(x, I, 1) = 1 then FAIL;
else x/2;
fi;
end:
# a = u+I*v = (x+I*y)^2 that is u = x^2-y^2, v = 2*x*y
# or x^4-u*x^2-v^2/4 = 0; with delta = u^2+v^2, this leads to
# x^2 = (u+/-sqrt(delta))/2
Gsqrt := proc(a)
local u, v, delta, rd, tmp, x;
u:=coeff(a, I, 0);
v:=coeff(a, I, 1);
delta:=(u*u+v*v);
rd:=isqrt(delta);
printf("sqrt(%a)=%a\n", delta, rd);
tmp:=u+rd;
# printf("u+sq=%a\n", tmp);
if tmp mod 2 = 1 then
ERROR("odd");
RETURN(FAIL);
fi;
tmp:=iquo(tmp, 2);
x:=isqrt(tmp);
if x*x <> tmp then ERROR("sq"); fi;
# if debug:
# print("x="+str(x));
if v mod (2*x) <> 0 then ERROR("y"); fi;
y:=iquo(v, (2*x));
# print("square root in G: "+str(G(x, y)));
# return G(x, y);
x+I*y;
end:
# INPUT: LHS = st^2 -4*n = no = r^2, where LHS = 4*n, n = (1+I)*n0
# OUTPUT: the factors p and q of n0 (and not n).
# TODO: use known squareroot of norm
find_p_q := proc(st, LHS, no)
local r, p, q, u;
r:=Gsqrt(LHS);
printf("r:=%a;\n", r);
printf("chk: %a\n", evalc(r*r-LHS));
# st^2-4*n = r^2 or 4*n = r^2-st^2
# up to +/- r?
# up to exchange of s and t?
# r-st=2*q
# r+st=2*(I+1)*p
p:=r+st;
printf("p0=%a\n", p);
p:=Gdiv2(p);
q:=r-st;
printf("q0=%a\n", q);
q:=Gdiv2(q);
if p = FAIL or q = FAIL then
ERROR("p and q are odd");
fi;
printf("p1:=%a;\n", p);
printf("q1:=%a;\n", q);
# p:=p / (1+I):=p*(1-I)/2
u:=Gdiv2(p*(1-I));
if u = FAIL then
printf("swap p and q\n");
u:=p; p:=q; q:=u;
u:=Gdiv2(p*(1-I));
if u = FAIL then
ERROR("Gasp");
fi;
fi;
p:=u;
printf("p2=%a\n", p);
printf("q2=%a\n", q);
# we might want p >> 0?
if coeff(p, I, 0) < 0 and coeff(p, I, 1) < 0 then
p:=-p; q:=-q;
fi;
[p, q];
end:
IFermat := proc(n0, S, T)
local n, x0, s0, t0, s, t, no, LHS, dig, tmp, tmp0, ds, dt;
dig:=400;
n:=evalc(n0 * (1+I));
x0:=evalf(sqrt(abs(4*n))*exp(I*argument(4*n)/2), dig);
print(evalf(4*n-x0^2, dig));
s:=floor(coeff(x0, I, 0));
t:=floor(coeff(x0, I, 1));
tmp:=evalc(4*n-(s+I*t)^2);
printf("s=%a, t=%a => %a\n", s, t, tmp);
s0:=s; t0:=t; tmp0:=tmp;
s:=floor(coeff(x0, I, 0));
t:=ceil(coeff(x0, I, 1));
tmp:=evalc((s+I*t)^2-4*n);
printf("s=%a, t=%a => %a\n", s, t, tmp);
if abs(tmp)^2 < abs(tmp0)^2 then
printf("new best fc\n");
s0:=s; t0:=t; tmp0:=tmp;
fi;
s:=ceil(coeff(x0, I, 0));
t:=floor(coeff(x0, I, 1));
tmp:=evalc((s+I*t)^2-4*n);
printf("s=%a, t=%a => %a\n", s, t, tmp);
if abs(tmp)^2 < abs(tmp0)^2 then
printf("new best cf\n");
s0:=s; t0:=t; tmp0:=tmp;
fi;
s:=ceil(coeff(x0, I, 0));
t:=ceil(coeff(x0, I, 1));
tmp:=evalc((s+I*t)^2-4*n);
printf("s=%a, t=%a => %a\n", s, t, tmp);
if abs(tmp)^2 < abs(tmp0)^2 then
printf("new best ff\n");
s0:=s; t0:=t; tmp0:=tmp;
fi;
for ds from 0 to S
do
if ds mod 10^2 = 0 then
printf("ds=%a at %a\n", ds, time());
fi;
s:=s0+ds;
for dt from 0 to T
do
t:=t0+dt;
LHS:=(s+I*t)^2-4*n;
no:=Inorm(LHS);
if issqr(no) then
printf("ds:=%a; dt:=%a; s:=%a; t:=%a;\n", ds, dt, s, t);
RETURN(find_p_q(s+I*t, LHS, no));
fi;
od;
od;
FAIL;
end:
chal_n:=-1407117267259556515600776848775735215276485947638218330895429449228218\
6438059768663873855254160723159027292054036744986124692686209556946009925974804\
689837 - 6610317000442669421334350882238890617453420384375090003163635532360574\
3323313950321796270758826382118980344498322657078158896922616675836987252304039\
22888*I;
IFermat(chal_n, 10^6, 10^6);