/*
** invphi.gp ver. 1.3 (c) 2005,10 by Max Alekseyev
**
** invphi(n)	computes all solutions to the equation eulerphi(x)=n with respect to x.
** numinvphi(n)	computes the number of such solutions and is faster than #invphi(n).
** invsigma(n)	computes all solutions to the equation sigma(x)=n with respect to x.
** invphitau(n,m) computes all solutions to the equations eulerphi(x)=n, numdiv(x)=m with respect to x.
*/


\\ inverting multiplicative function   p^k -> (p^(k+1)-1)/(p-1)
{ invsigma(n,m=3) = local(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,#f,
      p = f[i];
      P = d*(p-1)+1;
      k = valuation(P,p) - 1;
      if( k<1 || P!=p^(k+1), next );
      R = concat(R, p^k*select(x->(x%p),invsigma(n\d,d)) );
    );
  );
  vecsort(R,,8)
}


{ numinvphi(N) = local(P, L, p, l, D, r, t); 
  P=Set(); 
  L=listcreate(); 

  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)) );
  );

  \\ dynamic programming
  D = Set(divisors(N));
  r = vector(#D);
  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],n,
        t[setsearch(D,n*L[i][j])] += r[setsearch(D,n)];
      );
    );
    r = t;
  );
  listkill(L); 
  r[setsearch(D,N)];
}


{ invphi(N) = local(P, L, p, l, D, r, t); 
  P=Set(); 
  L=listcreate(); 

  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));
  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]),,8);
      );
    );
    r = t;
  );
  listkill(L); 
  r[setsearch(D,N)]
}


\\ M is the number of divisors
{ invphitau(N,M) = local(P, L, p, l, D, r, t); 

  P=Set(); 
  L=listcreate(); 

  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], select(x->(M%x[3])==0, vector(valuation(N,p)+1,i,[n*p^(i-1),p^i,i+1]) ) );
  );

  \\ dynamic programming
  D = Set(divisors(N));
  r = matrix(#D,M,i,j,[]);
  r[setsearch(D,1),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]);
        fordiv(M/L[i][j][3],m,
          t[l,m*L[i][j][3]] = vecsort(concat(t[l,m*L[i][j][3]],r[setsearch(D,n),m]*L[i][j][2]),,8);
        );
      );
    );
    r = t;
  );
  listkill(L); 
  r[setsearch(D,N),M]
}