"Forget I exist"
Jul 2009
Dumbassville
2^{6}×131 Posts

Quote:
Originally Posted by mart_r
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<3u<1,break());t=0;if(u==0,t=1);if(u>0,t=u);r=0;if(t==1,q=2;while(!ispseudoprime(nq1),q=nextprime(q+1);if(q>n,break()));r=q*(nq1),if(u==1,t=n%2);while(r<=0&&u!=0&&t<sqrt(n),if(u==1,t+=2);s=sigma(t);m=factor((st)*(ns)+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)/(st);p=(ns*(q+1))/((st)*q+s);if(qfloor(q)==0&&ispseudoprime(floor(q))&&pfloor(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(n1),r=(n1)*(n1));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 nq1 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.)

I should be thanking you, it was trying to come up with a "better" ali after looking and responding to this thread that I realized a property of partitions() I hadn't before allowing a partial solution up to a given number ( like vecmax would) and it allowed me to practice with select() as well.
here's what I came up with:
Code:
ali2000(y)=
a=select(v>isprime(v[1])==1 &&
#vecsort(v,,8)==#v &&
#setminus(Vec(divisors(v[#v])),Vec(v))==1,
partitions(y1));
for(x=1,#a,
if(sigma(a[x][1]*a[x][#a[x]])(a[x][1]*a[x][#a[x]])==y,
a[x]=a[x][1]*a[x][#a[x]],
a[x]=0)
);
a=vecsort(a,,8);
a=vector(#a1,n,a[n+1])
Last fiddled with by science_man_88 on 20130723 at 20:54
