Internally, we have a function Zn_quad_roots that compute all the solution of x^2+b*x+c mod N
for composite N.
Maybe we could add it to GP if we find a GP interface to it.
Bill, I'd truly appreciate having such a function.
On a related note, below is my naive function computing all square roots of a given t_INTMOD. Somehow it badly loses to SageMath's built-in sqrt() function that provides an ability to get all the roots out of the box.
While there is no built-in alternative in PARI/GP, is there a way to speed up my function?
Thanks,
Max
? sum(r=0,10^5, #asqrtintmod(Mod(r,1000000000000001000000000000001)) )
%14 = 111953
? ##
*** last result: cpu time 16,544 ms, real time 16,580 ms.
sage: %time sum(len(sqrt(Mod(r,1000000000000001000000000000001),all=True,extend=False)) for r in (0..10^5))
CPU times: user 1.92 s, sys: 280 µs, total: 1.92 s
Wall time: 1.92 s
111953
===
{ asqrtintmod(q) = my(p, r, t, l);
if(!issquare(q),return([]));
p = factor(q.mod);
l = matsize(p)[1];
t = vector(l,i, sqrt(lift(q) + O(p[i,1]^p[i,2])) );
t = apply(x->Mod(x,x.p^padicprec(x,x.p)), t);
r = List();
forvec(s=vector(l,i,[0,1]),
listput(r, chinese(vector(l,i,t[i]*(-1)^s[i])) );
);
Set(r);
}
===