Zeta function with Maxima/Macsyma in imaxima.el: 2010

Dec 12, 2010

Totient function and its identity

Maxima 5.22.1 http://maxima.sourceforge.net
using Lisp CMU Common Lisp Snapshot 2010-05 (20A Unicode)
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) n=dsum(totient(d),d,n);

(%i2) listify(divisors(24));

(%i3) map(totient, %);

(%i4) apply("+", %);

(%i5) pair2frac(x):=first(x)/second(x);

(%i6) orderbydenom(x,y):=orderlessp(denom(x),denom(y));

(%i7) create_list([k,n],k,1,n),n=10;

(%i8) map(pair2frac,%);

(%i9) sort(%,orderbydenom);

(%i10) create_list([k,n],k,1,n),n=24;

(%i11) map(pair2frac,%);

(%i12) sort(%,orderbydenom);

(%i13)

May 23, 2010

Generating function of Mangoldt function

From the Mangoldt function definition with dsum(), it is easy to obtain the generating function of Mangoldt function.

The following file contains the definition of dsum().
(%i1) load("elementaryNumberTheory.mac");

(%i2) mangoldt(N):=if integerp(N)
then apply(dsum,[moebius(N/d)*log(d),d,N])
else dsum(moebius(N/d)*log(d),d,N)$

(%i3) texput(nounify(mangoldt), "\\Lambda")$

The definition of Mangoldt function is given above.
(%i4) 'mangoldt(n)=mangoldt(n);

We can see that the definition works fine.

(%i5) 'mangoldt(27)=mangoldt(27);

(%i6) 'mangoldt(27)=radcan(mangoldt(27));

(%i7) 'mangoldt(28)=mangoldt(28);

(%i8) 'mangoldt(28)=radcan(mangoldt(28));

The following holds for any number-theoretic functions:

(%i9) dprod;

Having a(n) being Moebius function and b(n) being logarithmic function, the following holds:
(%i10) dprod,a(n):=moebius(n),b(n):=log(n);

Since (%o11) and (%o12) holds, and right hand side of (%o10) is already the generating function of Mangoldt function, we obtain (%o13).

(%i11) 'diff(zeta(s),s)=sum(diff(1/n^s,s),n,1,inf);

(%i12) 1/zeta(s)=sum(moebius(n)/n^s,n,1,inf);

(%i13) -diff(zeta(s),s)/zeta(s)=sum('mangoldt(n)/n^s,n,1,inf);

(%i14)

Mantoldt function

This article defines the Mangoldt function. The file of this session can be found here.

The following file defines the function dsum(). You can find the definition of dsum() in the previous article. 
(%i1) load("elementaryNumberTheory.mac");

(%i2) matchdeclare([a,b],true)$

(%i3) defrule(logexpon, a*log(b), log(b))$
The above rule logexpon defines a pattern matching and transformation function that appears in the standard definition of the Mangoldt function, which is as follows:

(%i4) mangoldt1(N):= if integerp(N)
then block([exp], if (exp:logexpon(radcan(log(N)))) = false then 0 else exp)
else apply(nounify(mangoldt1), [N])$

The above defines that if N is power of a prime, then the log of the prime is returned, else 0 is returned.

(%i5) mangoldt(N):=if integerp(N)
then apply(dsum,[moebius(N/d)*log(d),d,N])
else dsum(moebius(N/d)*log(d),d,N)$
The above defines that mangoldt function of N is defined to be sum of convolution multiplication of moebius() and log() over the divisors of N.
(%i6) texput(nounify(mangoldt), "\\Lambda")$

(%i7) a:radcan(log(11^3));

(%i8) logexpon(a);

(%i9) 'mangoldt1(11^3)=mangoldt1(11^3);

The above shows that the mangoldt1() works fine if the number is power of a prime.

(%i10) a:radcan(log(2^2*3^3));

(%i11) logexpon(a);

(%i12) 'mangoldt1(2^2*3^3)=mangoldt1(2^2*3^3);

The above shows that mangoldt1() works fine for numbers other than the power of prime.

The following session shows that the definition based on dsum() function works equally.

(%i13) a:11^3;

(%i14) 'mangoldt(a)=mangoldt(a);

(%i15) radcan(%);

(%i16) a:2^2*3^3;

(%i17) 'mangoldt(a)=mangoldt(a);

(%i18) radcan(%);

(%i19)

Mar 27, 2010

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)

Mar 14, 2010

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.

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