Product over primes

It is possible in Maxima to define pproduct which is the product over primes, similar to ordinary product(). Because I want the evaluation rule of the arguments of pproduct to be the same as product, I used buildq.

Also, the texput is used to define the output of TeX so that it can be rendered as expected in Imaxima. You can obtain the file pproduct.mac 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) 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])$

(%i2) 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)))))$

The function pproduct can be used to write down the famous Euler product of the Riemann zeta function.

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


(%i4) subst(10,inf,rhs(%));


(%i5) ev(%,nouns);


(%i6)

Value of zeta(2)

Here shows that

using the sine factorization formula. The session file 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/test-imaxima/cmucl/5.20.0/share/maxima/5.20.0/emacs/imaxima.lisp")), linenum:0)$

(%i1) powerdisp:true;


(%i2) sin(%pi*x)=%pi*x*product(1-x^2/n^2,n,1,inf);


(%i3) sin(%pi*x)=taylor(sin(%pi*x),x,0,10);


(%i4) coeff(rhs(%),x,3);

We see how the first n terms of the infinite product is expanded with n=2,3,4,5,6, and 10 and seet the coincidence of coeffcients of x^3 to the sum of first n terms of the infinite sum definition of zeta function.

(%i5) %pi*x*'product(1-x^2/n^2,n,1,k)=expand(%pi*x*product(1-x^2/n^2,n,1,k)),k=2;


(%i6) 'sum(1/n^2,n,1,k)=sum(1/n^2,n,1,k),k=2;


(%i7) %pi*x*'product(1-x^2/n^2,n,1,k)=expand(%pi*x*product(1-x^2/n^2,n,1,k)),k=3;


(%i8) 'sum(1/n^2,n,1,k)=sum(1/n^2,n,1,k),k=3;


(%i9) %pi*x*'product(1-x^2/n^2,n,1,k)=expand(%pi*x*product(1-x^2/n^2,n,1,k)),k=4;


(%i10) 'sum(1/n^2,n,1,k)=sum(1/n^2,n,1,k),k=4;


(%i11) %pi*x*'product(1-x^2/n^2,n,1,k)=expand(%pi*x*product(1-x^2/n^2,n,1,k)),k=5;


(%i12) 'sum(1/n^2,n,1,k)=sum(1/n^2,n,1,k),k=5;


(%i13) %pi*x*'product(1-x^2/n^2,n,1,k)=expand(%pi*x*product(1-x^2/n^2,n,1,k)),k=6;


(%i14) 'sum(1/n^2,n,1,k)=sum(1/n^2,n,1,k),k=6;


(%i15) %pi*x*'product(1-x^2/n^2,n,1,k)=expand(%pi*x*product(1-x^2/n^2,n,1,k)),k=10;


(%i16) 'sum(1/n^2,n,1,k)=sum(1/n^2,n,1,k),k=10;

Thus the following equation holds:

(%i17) -sum(1/n^2,n,1,inf)*%pi=-%pi^3/6;


(%i18) %/(-%pi);


(%i19)

Euler's Cotangent and Sine factorization formulae

The following session illustrates how to obtain Euler's cotangent formula and sine factorization formula:

They play important role relating to zeta function, which will be treated in other posts in the future.

The file EulerCot.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/test-imaxima/cmucl/5.20.0/share/maxima/5.20.0/emacs/imaxima.lisp")), linenum:0)$

(%i1) kill(all);


(%i1) load(fourie)$

(%i2) "In this file, the following two equations are proved.
They are originated by Euler."$

(%i3) %pi*cot(%pi*t)=2*t*sum(1/(t^2-n^2),n,1,inf)+1/t;


(%i4) sinfact:sin(%pi*t)=%pi*t*product(1-t^2/n^2,n,1,inf);


(%i5) assume(t>0);


(%i6) assume(t-n>0);


(%i7) declare(n,integer);


(%i8) an:fourier(cos(t*x),x,%pi);





(%i11) an0:rhs(ev(first(an)));


(%i12) an2:factor(trigexpand(rhs(ev(second(an)))));

As the Fourier coefficients are obtained for cos(t x), we obtain the following formula.

(%i13) fo:cos(t*x)=an0+sum(an2*cos(n*x),n,1,inf);


(%i14) fo2:fo,x=%pi;


(%i15) fo3:fo2*%pi/sin(%pi*t);


(%i16) eucot:trigreduce(expand(fo3));



(%i18) "With this, we can show the factor of sine function:"$

(%i19) sinfact;


(%i20) "Taking the differential of log of left hand side:"$

(%i21) trigreduce(diff(log(lhs(sinfact)),t));


(%i22) "Taking the differential of log of right hand side:"$

(%i23) expand(factor(diff(log(rhs(sinfact)),t)));


(%i24) "According to the Euler's cotangent formula, the previous
two formula are eaual. Hence,"$

(%i25) sf1:sin(%pi*t)=C*t*product(1-t^2/n^2,n,1,inf);


(%i26) sf2:limit(sf1/t,t,0);

Maxima cannot compute this limit. We replace inf with 100 and how the limit is computed.

(%i27) sf3:subst(100,inf,sf2);


(%i28) sf3,nouns;


(%i29) "Hence we obtain:"$

(%i30) sf1, C=%pi;

Let's take a look at the graph of left and right hand side of the equation to make sure this holds.

(%i31) wxdraw2d(explicit(lhs(sinfact),t,-5,5), explicit(ev(subst(100,inf,rhs(sinfact)),simpproduct),t,-5,5));



(%i32)