/* ** 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(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] }