ultimate_perfect_mod(p,m,N,aa,bb):= for a:aa thru bb 
do(mm:mod(a,N),if mm>0 then ( mxp:maxpp(a),w:(p-1)*divsum(a)-p*a-(p-2)*mxp+m*(p-1),
if w=0 then (print(a),display(factor(a))) else 1=1));

ultimate_perfect(p,m,aa,bb):= for a:aa thru bb 
do(mxp:maxpp(a),w:(p-1)*divsum(a)-p*a-(p-2)*mxp+m*(p-1),
if w=0 then (print(a),display(factor(a))) else 1=1);

maxpp(a):= (c:ifactors(a),cc:last(c), c1:cc[1]);

perfect_pms(p,m,aa,bb):= for e:aa thru bb do(w:p-1,x:(p^(e+1)-1)/w+m,
if factor(x)=x then (y:p^e*x,display(factor(y)),print(y)) else 1=1);

perfect_pms(2,0,2,80);

perfect_pms(p,m,aa,bb):= for e:aa thru bb do(w:p-1,x:(p^(e+1)-1)/w+m,
if factor(x)=x then (y:p^e*x,display(factor(y)),print(y)) else 1=1);

perfect_pms(2,2,1,30);

perfect_pms(2,8,2,50);

perfect_pms(2,0,2,50);

6