Sum of Moebius function of n over its divisors

This text shows the calculation of sum of möbius function of n over the divisors of n.

Maxima 5.20.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 19f (19F)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) block(load(("/Users/yasube/Programming/imaxima-imath/imaxima.lisp")), linenum:0)$

(%i1) dsum(exp, var, N) ::=
if integerp(N)
then buildq([exp:exp,var:var,N:N],block([res:0,divlist:listify(divisors(N))], for var in divlist do res:res+exp, return(res)))
else apply(nounify(dsum),[exp,var,N])$

(%i2) texput(nounify(dsum),
lambda([arglist], block([e,v,pm],[e,v,pm]:args(arglist),
concat("\\sum_{",tex1(v) ,"|",tex1(pm),"}",tex1(e)))))$

(%i3) texput(moebius, "\\mu");

We already learned the following equation, representing the Dirichlet series generating function for möbius function.

 

(%i4) moe0:sum(moebius(n)/n^s,n,1,inf)=1/zeta(s);


(%i5) moe1:moe0*zeta(s);


(%i6) subst(sum(1/n^s,n,1,inf),zeta(s),moe1);


 

We also learned that the multiplication of two Dirichlet series results also a Dirichlet series of convolution product.

(%i7) moe2:sum(a(n)/n^s,n,1,inf)*sum(b(n)/n^s,n,1,inf)=sum(dsum(a(d)*b(n/d),d,n)/n^s,n,1,inf);


 

Using this, we can calculate the left hand side of (%o6).

(%i8) moe3:moe2,a(n):=moebius(n),b(n):=1;


(%i9) rhs(moe3)=1;


 

 Hence, if n=1,

 (%i10) dsum(moebius(d),d,n)=1;


if n>1 then,

(%i11) dsum(moebius(d),d,n)=0;


(%i12)

 

dsum : Sum over divisors

This article demonstrates the defining of dsum, sum over divisors on given number N. As before, the tex output definition is also presented using texput(). Few of the trivial examples are given, too.

Maxima 5.20.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 19f (19F)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) block(load(("/Users/yasube/Programming/imaxima-imath/imaxima.lisp")), linenum:0)$

(%i1) kill(all);


(%i1) dsum(exp, var, N) ::=
if integerp(N)
then buildq([exp:exp,var:var,N:N],block([res:0,divlist:listify(divisors(N))], for var in divlist do res:res+exp, return(res)))
else apply(nounify(dsum),[exp,var,N])$

(%i2) texput(nounify(dsum),
lambda([arglist], block([e,v,pm],[e,v,pm]:args(arglist),
concat("\\sum_{",tex1(v) ,"|",tex1(pm),"}",tex1(e)))))$

(%i3) texput(moebius, "\\mu")$

(%i4) "A couple of examples of dsum()"$

(%i5) 'divisors(36)=divisors(36);


(%i6) 'dsum(d,d,36)=dsum(d,d,36);


(%i7) 'divisors(20)=divisors(20);


(%i8) 'dsum(f(d),d,20)=dsum(f(d),d,20);


(%i9) exp1:sum(f(n)/n^s,n,1,inf)*sum(g(n)/n^s,n,1,inf)=sum(dsum(f(d)*g(n/d),d,n)/n^s,n,1,inf);


(%i10)

Integration formula of zeta

In the very first page of the original Riemann's paper on zeta function, the integration formula of zeta function is presented. This formula is shown below.

Maxima 5.20.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 19f (19F)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) block(load(("/Users/yasube/Programming/imaxima-imath/imaxima.lisp")), linenum:0)$

(%i1) kill(all);


(%i1) "We start from the definition of Gamma function. Instead of using
the name Gamma, we use Pi(s) since Gamma is already defined as a
Maxima function."$

(%i2) f1:Pi(s)='integrate(exp(-x)*x^s,x,0,inf);


(%i3) f2:subst(s-1,s,f1);


(%i4) assume(n>0);


(%i5) f3:changevar(f2,x-n*t,t,x);


(%i6) f4:f3/n^s;


(%i7) "a[n] is an arbitrary number series."$

(%i8) f5:f4*a[n];


(%i9) f6:sum(f5,n,1,inf);


(%i10) "We want to evaluate the right hand side of the previous
equation. We start by exchaging the integration and summation."$

(%i11) f7:integrate(sum(first(rhs(f4))*a[n],n,1,inf),t,0,inf);


(%i12) "Let's define F(x) as follows:"$

(%i13) f71:F(x)=sum(a[n]*x^n,n,1,inf);


(%i14) "Then, we obtain the mellin transformation."$

(%i15) f72:integrate(t^(s-1)*F(exp(-t)),t,0,inf);


(%i16) f8:lhs(f6)=f7;


(%i17) "The simplest application if to assume that a[n] = 1 for all n.
We obtain the zeta function in the left side. Let's see how we
can evaluate the right hand side further."$

(%i18) f9:subst(1,a[n],f8);


(%i19) assume(exp(t)-1>0);


(%i20) f10:factor(f9),simpsum=true;


(%i21)

 

Mobius function and its generating function

The pproduct() can be used to show that the Mobius function's generating function is 1/zeta(s).

This can be easily shown by expanding first few terms. The file mobius.mac can be obtained from here.

Maxima 5.20.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 19f (19F)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) block(load(("/Users/yasube/Programming/imaxima-imath/imaxima.lisp")), linenum:0)$

(%i1) matchdeclare([a,b,c],true)$

(%i2) defrule(dsimp1,a/(b^s*c^s),a/(b*c)^s)$

(%i3) pproduct(e, var, pmax) ::=
if integerp(pmax)
then buildq([e:e,v:var,pm:pmax], block([res:1],for v:next_prime(1) while (v else apply(nounify(pproduct),[e,var,pmax])$

(%i4) texput(nounify(pproduct),
lambda([arglist], block([e,v,pm],[e,v,pm]:args(arglist),
if (pm=inf) then
concat("\\prod_{",tex1(v) ,":prime}",tex1(e))
else
concat("\\prod_{",tex1(v) ,":prime <",tex1(pm),"}",tex1(e)))))$

(%i5) texput(moebius, "\\mu");


(%i6) moe0:sum(moebius(n)/n^s,n,1,inf)=1/zeta(s);


(%i7) moe1:subst(15,inf,lhs(moe0));


(%i8) moe2:moe1,nouns;


(%i9) zeta(s)=pproduct(1/(1-1/p^s),p,inf);


(%i10) z1:1/zeta(s)=pproduct(1-1/p^s,p,inf);


(%i11) z2:subst(15,inf,rhs(z1));


(%i12) z3:z2,nouns;


(%i13) z4:expand(z3);


(%i14) z5:apply1(z4,dsimp1)$

(%i15) rest(z5,length(z5)-length(moe2));


(%i16)