ez:=proc(): print(`NR(f,z,z0,eps), pn(E), DP(A,B,M)  `): end:

NR1:=proc(f,z,z0):
evalf(z0-subs(z=z0,f)/subs(z=z0,diff(f,z))):
end:

NR:=proc(f,z,z0,eps) local z1,z2,z3:

z1:=z0:
z2:=NR1(f,z,z1):

while abs(z1-z2)>eps do
 z3:=NR1(f,z,z2):
 z1:=z2:
 z2:=z3:
od:
z3:
end:


with(linalg):
#pn(E)
pn:=proc(E) local i,j,n:

n:=nops(E):

det([seq([i*E[i],seq(E[i-j],j=1..i-1),1,0$(n-i-1)],i=1..n-1),[n*E[n],seq(E[n-j],j=1..n-1) ]]):

end:

with(numtheory):

#DP(A,B,M): the product of two Dirichlet polynomials given by their list of coefficients, up to the M-th term. Try:
#DP([a1,a1],[b1,b2],4);
DP:=proc(A,B,M) local N,n,p,lu,lu1,c:

N:=min(nops(A)*nops(B),M):

p:=[]:

for n from 1 to N do

 lu:=divisors(n):


c:=0:
 
for lu1 in lu do

   if lu1<=nops(A) and n/lu1<=nops(B) then
    c:=c+A[lu1]*B[n/lu1]:
   fi:
od:
    p:=[op(p),c]:   
od:
p:

end: