Okay, there are one or two things that I didn't think of off the top of my head. Thanks.

This is a broader approach to the problem than I had in mind, I was more concerned about finding one (preferably the smallest) solution to a rather large (>10^12) given number.

Anyway the program I compiled during the weekend seems to work quite nicely already, though it still won't find a good solution for, say, n=Xpqr for large primes p,q,r and X a small prime or a product of small primes.

This is how far I've gotten on my own:

Code:

ai(n,u)=if(n<3||u<-1,break());t=0;if(u==0,t=1);if(u>0,t=u);r=0;if(t==1,q=2;while(!ispseudoprime(n-q-1),q=nextprime(q+1);if(q>n,break()));r=q*(n-q-1),if(u==-1,t=n%2);while(r<=0&&u!=0&&t<sqrt(n),if(u==-1,t+=2);s=sigma(t);m=factor((s-t)*(n-s)+s*s);v=[1];for(z=1,m[1,2],v=concat(v,m[1,1]^z));for(y=2,#m[,1],w=#v;for(z=1,m[y,2],for(x=1,w,v=concat(v,v[x]*m[y,1]^z))));for(x=1,#v,q=(v[x]-s)/(s-t);p=(n-s*(q+1))/((s-t)*q+s);if(q-floor(q)==0&&ispseudoprime(floor(q))&&p-floor(p)==0&&ispseudoprime(p),r=t*p*q;if(sigma(r)-r!=n,r=0);if(r>0,x=#v)));if(u>0,u=0)));if(u==-1&&r==0&&isprime(n-1),r=(n-1)*(n-1));if(r>0,print("solution: "r);print(factor(r)),print("no solution u*p*q for u="t))
\\ finds N such that sigma(N)-N=n via N=pq if n-q-1 is prime, optional u for additional factors N=upq; for u=-1, searches through u=2...inf until a solution is found

(Will modify it a little more later, for example trying a number of "special" small factors first.)