/*
inverse phi and sigma Max Alekseyev:
http://www.cse.sc.edu/~maxal/gpscripts/

ported to current pari/gp by joro
*/
{ invphi(N) = local(P, L, p, l, D, r, t); 
  P=Set(); 
  L=listcreate(10^5); 

  fordiv(N,n,
      p = n+1;
      if( !ispseudoprime(p), next );

      l=setsearch(P, p); 
      if(l==0, 
         l=setsearch(P, p, 1); 
         P=setunion(P, [p]); 
         listinsert(L, [], l);
      ); 
      L[l] = concat(L[l], vector(valuation(N,p)+1,i,[n*p^(i-1),p^i]) );
  );

  \\ dynamic programming
  D = Set(divisors(N));
  \\print("#D=",#D);
  r = vector(#D,i,[]);
  r[setsearch(D,1)] = [1];
  for(i=1,#L,
    t = r;       \\ stands for 1 in (1 + terms of L) 
    for(j=1,#(L[i]),
      fordiv(N/L[i][j][1],n,
        l = setsearch(D,n*L[i][j][1]);
        t[l] = vecsort(concat(t[l],r[setsearch(D,n)]*L[i][j][2]));
	\\print(t[l]);
      );
    );
    r = t;
    \\print(r);
  );
  listkill(L); 
  r[setsearch(D,N)]
}

{ invsigma(n,m=3) = local(ii,vv,v2,i,k,p,f,R=[]); 
  if(n==1, return([1]));
  fordiv(n,d,
    if(d<m,next);
    \\ checking that d is of the form d = 1 + p + p^2 + ... + p^k for some prime p and integer k
    f = factor(d-1)[,1];
    for(i=1,length(f),
      p = f[i];
      P = d*(p-1)+1;
      k = valuation(P,p) - 1;
      if( k<1 || P!=p^(k+1), next );
      vv=invsigma(n\d,d);
      v2=[];
      for(ii=1,length(vv),if(vv[ii]%p,v2=concat(v2,vv[ii])));
      R = concat(R, p^k*v2 );
    );
  );
  vecsort(R)
}