multfac(fm, fn)=
{
	if(fm==matrix(0,2), return(fn));
	matreduce(concat(fm~, fn~)~);
}
addhelp(multfac, "multfac(fm, fn): Multiplies two factored numbers together. If fm = factor(m) and fn = factor(n), then multfac(fm, fn) == factor(m*n).");

divisorsmax(n)=
{
	if(n<360,return(20));
	localprec(38);
	my(L=log(n),nats=L-log(log(L)),bits=nats\log(2));
	localprec(bits+16);
	n^(1.538*log(2)/log(log(n)))\1;
}
addhelp(divisorsmax, "divisorsmax(n): Gives an upper bound on the number of divisors of all m <= n.");

A348172(n)=
{
	if(n<4, return(2));
	my(d=numdiv(n),f=n/d,num=numerator(f),fnum=factor(num),denum=denominator(f),s,L);
	
	L=if(n<369,888,solve(x=n,n^2, x/divisorsmax(x)-f)\1);
	\\print("Searching up to "L);
	forfactored(k=1,L\num,
		my(t=k[1]*num,tfac=multfac(k[2],fnum),tdiv=numdiv(tfac));
		if(f==t/tdiv,s++)
	);
	s;
}