Three Sons Problem

magic(n) = n^2*(n^2+1)*(n^2-2*n+2);

three_sons(n) = {
	local(nfactor, f, count_outer, count_inner, xx, zz, k, x, y, z, b, c, det, detrt, ret);

	if(moebius(n) != 0, return(0));

	nfactor = factor(n);

	f = core(n,1)[2];

	if(magic(f)<n, return(0));

	count_outer = 0;

	fordiv(f, xx,
		if(xx == 1, next);
		if(xx^3 > n, break);

		if(magic(xx)<n, next);

		zz = n \ xx^2;
		k = 2*xx + zz;

		count_inner = 0;

		fordiv(nfactor, x,
			if(x^3 >= n || 3*x >=k, break);

			b = -(k-x);
			c = n \ x;
			det = b^2 - 4*c;

			if(!issquare(det,&detrt), next);

			y=(-b-detrt)\2;
			z=(-b+detrt)\2;

			if(y<=x, next);

			\\print([x,y,z]);
			ret=[x,y,z,xx,f];
			count_inner++;

			if(count_inner>1,
				\\print("Two per k!");
				return(0));

			);

		if(count_inner==1, count_outer++);

		if(count_outer>1,
			\\print("Two per n!");
			return(0));

		);

	if(count_outer != 1, return(0) );

	\\print(n" "ret);
	return(1);
	}