alias(K=EllipticK):alias(Kc=EllipticCK):alias(E=EllipticE):alias(Ec=EllipticCE):alias(Li=polylog): alias(AG=GaussAGM):alias(p=pochhammer): jac:=numtheory[jacobi]; Phi2:=(x,y)->log((cos(Pi*y)+cosh(Pi*x))/(cosh(Pi*x)-cos(Pi*y)))/4/Pi-2/Pi*sum(cos(Pi*(2*n+1)*y)*cosh(Pi*(2*n+1)*x)/(1+exp(Pi*(2*n+1)))/(2*n+1),n=0..infinity): tphi:=proc(a,b,d,D) local SS;Digits:=D+10:SS:=(x,y)->log((cos(Pi*y)+cosh(Pi*x))/(cosh(Pi*x)-cos(Pi*y)))/4/Pi-2/Pi*sum(cos(Pi*(2*n+1)*y)*cosh(Pi*(2*n+1)*x)/(1+exp(Pi*(2*n+1)))/(2*n+1),n=0..infinity);sort(PolynomialTools[MinimalPolynomial](evalf(exp(Pi*SS(a,b)),D),d));end: cphi:=proc(x,y,a,b) local M;M:=log(a)/Pi/b; print(evalf[50](Phi2(x,y)=M));M;end: tphi8:=proc(a,b,d,D) local SS;Digits:=D+10:SS:=(x,y)->log((cos(Pi*y)+cosh(Pi*x))/(cosh(Pi*x)-cos(Pi*y)))/4/Pi-2/Pi*sum(cos(Pi*(2*n+1)*y)*cosh(Pi*(2*n+1)*x)/(1+exp(Pi*(2*n+1)))/(2*n+1),n=0..infinity);sort(PolynomialTools[MinimalPolynomial](evalf(exp(8*Pi*SS(a,b)),D),d));subs(_X=x,%);end: UN:=(theta,N)->Li(N,exp(I*theta)): dL1:=z->sum(Zeta(1,1-n)*log(z)^n/n!,n=1..infinity)-gamma(1)-1/2*gamma^2-Zeta(2)/2- gamma*(log(-log(z)))-1/2*log(-log(z))^2: Ld:=proc (m,d, s) local r, L; L := add(numtheory[jacobi](d, k)*Zeta(m, r, k/abs(d)), k = 1 .. abs(d)-1)/abs(d)^r; L:=diff(L,r\$m);subs(r = s, L); if type(s,numeric) then L:=evalf(%) fi;L;end proc; Lsd:=proc(m,d,s) local r,L;L:=add(numtheory[jacobi](d,k)*Zeta(0,r,k/abs(d)),k=1..abs(d)-1)/abs(d)^r; if m >0 then L:=diff(L,r\$m):fi;evalf[50](subs(r=s,L));end: LSD:=proc(m,d,s) local k,j;add(jac(d,k)*add(binomial(m,j)*Zeta(m-j,s,k/abs(d))*(-log(abs(d)))^j,j=0..m),k=1..abs(d)-1)/abs(d)^s; end; EQL:=(m,d,s)->zeta[m](d,s)=ie(Lsd(m,d,s)): Ls := proc (z, s) options operator, arrow; #s not integer Li(s, z) = sum(Zeta(s-n)*log(z)^n/factorial(n), n = 0 .. infinity)- GAMMA(1-s)*(-log(z))^(s-1) end proc: LL:=proc(d,s) add(numtheory[jacobi](d,k)*Zeta(0,s,k/abs(d)),k=1..abs(d)-1)/abs(d)^s;end: Lin := proc (z,N) options operator, arrow; Sum(Zeta(N-n)*log(z)^n/factorial(n),n = N .. infinity)+ Sum(Zeta(N-n)*log(z)^n/factorial(n), n = 0 .. N-2) +log(z)^(N-1)/(N-1)!*(harmonic(N-1)-log(-log(z)));end proc: Lse2 := proc (z) options operator, arrow; -aze(2)-log(2)*log(z)-sum(aze(2-n)*log(z)^n/factorial(n), n = 2 .. infinity) end proc: aze:=proc(s) local a: if s=1 then a:=log(2) else a:=(1-2^(1-s))*Zeta(s) fi;a;end: CC:=proc(theta,M) local m;Cl(1,theta+Pi)=add((-1)^(m+1)*aza(1-2*m)*theta^(2*m+1)/(2*m+1)!,m=1..M)-log(2)*theta;end: alias(sinc=unapply(sin(x)/x,x)): Ch:=proc(n,m) if m=0 then 1 else binomial(n,m) fi;end: b2:=n->binomial(2*n,n)/4^n: z11:=(s,t)->(-1)^(s-1)/GAMMA(s)*Int((log(x))^(s-1)*Li(t,x)/(1-x),x=0..1): f2w:=proc(a0,N,col) local P,W, i,r, x,y,xx,yy,w,k,f,zz: with(plots):Digits:=N+10:f:=proc(a) if a=0 then [1,0] elif a=1 then [0,1] elif a=2 then [-1,0] else [0,-1] fi; end;x[1]:=0;y[1]:=0; r:=ListTools[Reverse](map(f,convert(trunc(10^N*evalf[N](a0)),base,4)));for k from 2 to N do x[k]:=x[k-1]+r[k][1];y[k]:=y[k-1]+r[k][2] od;xx:=[seq([x[k],y[k]],k=1..N)]: P:=pointplot([0,1],style = point, symbol=cross, color = black,symbolsize = 100): W:=pointplot(xx, axes=none,style = point, symbol=circle, color =col,symbolsize = 1, labels = ["x(t)", "y(t)"], caption = typeset("A uniform random walk on the first ",N," decimal digits of ",a0, "."),font = [Times, Roman, 12]) ;display(W,P) end: anim:=proc(num,a,b,s) local P,k,COL; COL:=proc(n) if irem(n,3)=0 then blue elif irem(n,3)=1 then green else yellow fi;end; for k from a to b do P[k]:=f2w(num,trunc(s^k),COL(k)): od; display([seq(P[k],k=a..b)],insequence=true);end: ip:=N->parse(cat(2,seq(ithprime(n),n=2..N))): ic:=N->parse(cat(0,seq(n,n=1..N))): Lia:=(s,z)->z/GAMMA(s)*Int((-log(z))^(s-1)/(1-z*t),t=0..1): L1:=(s,z)->(1/2)*z*Int(exp(-t)*t^(s-1)*(log(t)-Psi(s))*coth((1/2)*t-(1/2)*ln(z)), t = 0 .. infinity)/GAMMA(s): LI := proc (s,z) local a;options operator, arrow;a:= Re(GAMMA(1-s)*(I^(1-s)*Zeta(0, 1-s, 1/2-((1/2)*I)*ln(-z)/Pi)+I^(s-1)*Zeta(0, 1-s, 1/2+((1/2)*I)*ln(-z)/Pi))/(2*Pi)^(1-s)); if type(s,float) then a:=evalf(%) fi;a; end proc: dist:=proc(f,b,N) local s,d,p,t;Digits:=N+5;t:=time();s:=trunc(b^N*evalf[round(N*log[10](b))](f)); s:=convert(s,base,b):d:=seq(numboccur(k,s),k=0..b-1);p:=convert([d],`+`);print(p,`digits in`,time()-t,`seconds`);evalf(d/p,6);end: Lambda:=proc(x,k) local j,ans; ans:=(-1)^k*log(abs(x))^k/k; if k>1 then ans:=ans+(k-2)!*sum((-1)^j/j!*log(abs(x))^j*Li(k-j,x),j=0..k-2);fi;ans;end: LN:=(n,x)->Li(n,-x)+(-1)^n*Li(n,-1/x)=2*sum(log(x)^(n-2*r)*Li(2*r,-1)/(n-2*r)!,r=1.. trunc(n/2))-log(x)^n/n!: L21n:=(x,n)->Zeta(n+1)-sum((-1)^(n-m)/(n-m)!*log(1-x)^(n-m)*Li(m+1,1-x),m=0..n): CF:=(x,n)->numtheory[cfrac](x,n): cl:=(n,theta)->Im(I^(n mod 2)*Li(n,exp(I*theta))): alias(Li=polylog);Li2:=(r,t)->Re(Li(2,r*exp(I*t))); Li3:=(r,t)->Re(Li(3,r*exp(I*t))): Cl2:=(r,t)->Im(Li(2,r*exp(I*t))): Cl4:=(r,t)->Im(Li(4,r*exp(I*t))): J:=BesselJ: H:=(x,y,z)->hypergeom(x,y,z): HH:=(a,b,c,z)->GAMMA(c)/GAMMA(b)/GAMMA(c-b)*Int(t^(b-1)*(1-t)^(c-b-1)/(1-z*t)^a,t=0..1): t3:=q->JacobiTheta3(0,q):t2:=q->JacobiTheta2(0,q): t4:=q->JacobiTheta4(0,q): kay:=N->(t2/t3)(exp(-Pi*sqrt(N)))^2: kc:=N->(t4/t3)(exp(-Pi*sqrt(N)))^2: a:=proc(q) t3(q)*t3(q^3)+t2(q)*t2(q^3);end: b:=proc(q) 3/2*a(q^3)-a(q)/2;end: c:=proc(q)1/2*a(q^(1/3))-1/2*a(q);end: ie:=identify@evalf: timed:=proc(X) local Y:global t; t:=time();print(t);Y:=X:t:=time();print(`cpu`,t);Y;end: Kay:=proc() local d,x,x1,a,b;d:=Digits; if nargs>1 then d:=args[2];fi; x:=evalf[d](args[1]);a:=evalf[d](Pi/2); while evalf(x-10^(1-d))>0 do x1:=sqrt(1-x^2); x:=(1-x1)/(1+x1);a:=a*(1+x);od;evalf[d](a);end: Ti2:=x->Int(arctan(y)/y,y=0..x); CL:=proc(r,t) local w; w:=arctan(r*sin(t)/(1-r*cos(t)));Cl2(r,t)=(Cl2(1,2*w)+Cl2(1,2*t)-Cl2(1,2*w+2*t))/2+w*log(r);end: CLa:=proc (r, t) local w; w := -Pi+arctan(r*sin(t)/(1-r*cos(t))); Cl2(r, t) = (1/2)*Cl2(1,2*w)+(1/2)*Cl2(1,2*t)-(1/2)*Cl2(1,2*w+2*t)+w*log(r) end: E1:=x->Li2(x,Pi/4)=1/4*Li2(sqrt(2)*x-x^2,0)-1/2*Li2(x/(x-sqrt(2)),0)+1/8*Li2(-x^2/(x-sqrt(2))^2,0):#x<1 E2:=(r,t)->Li2(r,t)=-Li2(1/r,t)+(Pi-t)^2/2-log(r)^2/2-Pi^2/6: E3:=(r,t)->Li3(r,t)=Li3(1/r,t)+(3*(Pi-t)^2-Pi^2)*log(r)/6-log(r)^3/6: E4:=t->Ti2(tan(t))=t*ln(abs(tan(t)))+1/2*Cl2(1,2*t)+1/2*Cl2(1,Pi-2*t): E5:=t->Li2(2*cos(t),t)=(Pi/2-t)^2: E6:=t->Li2(cos(t),t)=1/4*Li(2,cos(t)^2)+1/2*(Pi/2-t)^2: E7:=t->Li2(sec(t),t)=5*Pi^2/24-1/4*Li(2,cos(t)^2)-1/2*ln(cos(t))^2-1/2*Pi*t: E8:=t->Li2(sec(t)/2,t)=Pi^2/12-1/2*ln(2*cos(t))^2-1/2*t^2: E9 := t->Re(Li(3, 1/2-(1/2)*exp((2*I)*t))) = (7/16)*Zeta(3)+(1/8) *Li(3, sin(t)^2)+(1/2)*t^2*log(sin(t))-(1/4)*Re(Li(3, exp((2*I)*t)))-(1/8)*Re(Li(3, exp((4*I)*t))): E10:=t->Re(Li(3,1-(1)*exp(I*t)))=1/2*Zeta(3)-1/2*Re(Li(3,exp(I*t)))+1/4*t^2*log(2*sin(t/2)): E11:=(x,t)->Li2(x,t)+Li2(2*cos(t)-x,t)=(Pi/2-t)^2+1/2*Li(2,2*x*cos(t)-x^2): RED:=proc(Q,n,d) local k;collect(reduce(add(op(Q)[k],k=1..n),d)+add(op(Q)[k],k=n+1..nops(Q)),polylog):end: reduce:=proc(a,d) local b; b:=a;simplify(b,{test(expand(b),d)});end: test:=proc (eq) Pslq(lhs(eq), [op(rhs(eq))]) end proc: Test:=proc (eq,L) Pslq(eq, L) end proc: shorten:=proc(A,n,d,N) local ts,e,k,R,Q;Q:=A;e:=0; for k to N do print(k,"length",length(Q),nops(Q));ts:=time(); R:=RED(Q,n,d);e:=evalf[d+10](R-A); if e<10^(-d) then print(k,"error",evalf[10](e),"length",length(Q),nops(Q),"time",time()-ts);Q:=R; save(Q,temp2);else print("failure");return; fi;od;print("termination");Q;end: G:=proc(m,n,a,b,z) local p,q,k;p:=nops(a);q:=nops(b); MeijerG([[seq(a[k],k=1..n)],[seq(a[k],k=n+1..p)]],[[seq(b[k],k=1..m)],[seq(b[j],j=m+1 ..q)]],z);end: KK:=proc() local n,N, d,x,x1,a,b;d:=Digits;N:=6; if nargs>1 then N:=args[2];fi;x:=evalf[d](args[1]); a:=evalf[d](Pi/2); for n to N do x1:=sqrt(1-x^2); x:=evalf[d]((1-x1)/(1+x1));a:=a*(1+x);od;evalf[d](a);end: cs:=f->convert(f,StandardFunctions): ch:=f->convert(f,hypergeom): cv:=IntegrationTools[Change]: cg:=x->convert(x,GAMMA): L3:=sum(1/(3*n+1)^2-1/(3*n+2)^2,n=0..infinity): er:=proc(f) expand(rationalize(f));end: ab:=evalf@abs: REAL:=Q->convert(map(Re,[op(Q)]),`+`): UN:=proc(A) local t;subs(Re=(t->t),A):%;end: Cl:=proc(r,t); Im(polylog(2,r*exp(I*t)));end: CLa:=proc (r, t) local w; w := -Pi+arctan(r*sin(t)/(1-r*cos(t))); Cl2(r, t) = (1/2)*Cl2(1, 2*w)+(1/2)*Cl2(1, 2*t)-(1/2)*Cl2(1, 2*w+2*t)+w*log(r) end proc: Ti:=z->Sum((-1)^n*z^(2*n+1)/(2*n+1)^2,n=0..infinity): pa:=A->pslq(A,[Pi, arctan(2), arctan(sqrt(2))]): cr:=f->convert(evalf(f),rational): cc:=f->convert(evalf(f),confrac): fval:=proc(f) f=evalf[6](f) end: val:=proc(f) f=value(f) end: dp:=proc(a,b) local k:=nops(a);add(a[k]*b[k],k=1..nops(a));end: bisect:=proc (f, a0, b0, N) local l, m,mp, r, a, b, c, p,eps;Digits:=2*N;eps:=.1^N;a := a0; b := b0; l := f(a); r := f(b); while eps < b-a do c := (1/2)*a+(1/2)*b; m := f(c);mp:=f(c+eps); if m >mp then a := c; l:= m else b := c;r := m; end if end do; print('function-values',[l, r]); [a, b] end proc: reshape:=proc(L,n,m) local k;linalg[matrix](n,m,[seq(L[k],k=1..m*n)]);end: r2p:=proc(A) global x; sqrfree(evala(Norm(convert(x-A,RootOf))),x)[2][1][1];end: f2p:=proc() local A,d; A:=evalf(args[1]);if nargs>1 then d:=args[2];else d:=trunc(sqrt(A));fi; PolynomialTools[MinimalPolynomial](A,d);end: IND:=proc() local P,Q,eps,k,D; if nargs>1 then Digits:=args[2]; fi;eps:=10.^(Digits-3);P:=args[1]:Q:=IntegerRelations[PSLQ](P);print(nops(Q)); D:=add(Q[k]*P[k],k=1..nops(P));print(evalf(D)); if evalf(D)>eps then print('independent');P; elif Q[nops(P)]=0 then ListTools[Rotate](P,1) else [seq(P[n],n=1..nops(P)-1)];fi;end: BtoA := proc(n,B,A,U) local bk,j,k; description "n is a positive integer", "B is an input Vector(n) of numbers", "A is an output Vector(n)", "U is a storage Vector(n)", "1+sum(A[i]*q^i,i=1..n)=product((1-q^i)^(-B[i]),i=1..n)", "as a power series in q, up to the q^n term"; bk := B[1]; for j to n do U[j] := bk; od; for k from 2 to n do bk := k*B[k]; for j from k by k to n do U[j] := U[j]+bk; od; od; for k to n do A[k] := 1/k*(add(U[j]*A[k-j],j=1..k-1)+U[k]); od; n; end: s2p := proc(S,q) local n,A,B,U,i; description "S is a polynomial in q", "q is the name of the variable", "output is a q-product"; if not type(S,polynom(anything,q)) then error "S must be polynomial in q"; fi; n := degree(S,q); if tcoeff(S,q)-1<>0 then error "the trailing coefficient must be 1"; fi; A := Vector(n); for i to n do A[i] := coeff(S,q,i); od; B := Vector(n); U := Vector(n); AtoB(n,A,B,U); mul((1-q^i)^(-B[i]),i=1..n); end: AtoB := proc(n,A,B,U) local k,bk,j; description "n is a positive integer", "A is an input Vector(n) of numbers", "B is an output Vector(n)", "U is a storage Vector(n)", "1+sum(A[i]*q^i,i=1..n)=product((1-q^i)^(-B[i]),i=1..n)", "as a power series in q, up to the q^n term"; for k to n do U[k]:= k*A[k]-add(A[j]*U[k-j],j=1..k-1); od; bk := A[1]; B[1] := bk; for j to n do U[j] := U[j]-bk; od; for k from 2 to n do bk := U[k]; B[k] := bk/k; for j from k by k to n do U[j] := U[j]-bk; od; od; n; end: p2s := proc(P,q) local m,n,PP,A,B,U,i,T,T1,T2,Q,Q1,Q2; description "P is a q-product", "q is the name of the variable", "output is a q-series"; if type(P,`*`) then m := nops(P); PP := [op(P)]; else m := 1; PP := [P]; fi; n := 0; for i to m do T := PP[i]; if type(T,`^`) then T1 := op(1,T); T2 := op(2,T); else T1 := T; T2 := 1; fi; if has(T2,q) then error "q must not be in the exponent"; fi; Q :=1-T1; if type(Q,`^`) then Q1 := op(1,Q); Q2 := op(2,Q); if Q1<>q then error "must have base=q"; fi; if not type(Q2,posint) then error "q^z must have z integer"; fi; n := max(n,Q2); elif Q=q then n :=max(n,1); else error "expression must be a finite q-product"; fi; od; A := Vector(n); B := Vector(n); U := Vector(n); for i to m do T := PP[i]; if type(T,`^`) then T1 := op(1,T); T2 := op(2,T); else T1 := T; T2 := 1; fi; if has(T2,q) then error "q must not be in the exponent"; fi; Q := 1-T1; if type(Q,`^`) then Q1 := op(1,Q); Q2 := op(2,Q); if Q<>q then error "must have base=q"; fi; if not type(Q2,posint) then error "q^z must have z integer"; fi; if B[Q2]<>0 then error "2 of the factors have the same base"; fi; B[Q2] := -T2; elif Q=q then if B[1]<>0 then error "2 of the factors have the same base"; fi; B[1] := -T2; else error "expression must be a finite q-product"; fi; od; BtoA(n,B,A,U); 1+add(A[i]*q^i,i=1..n); end: Ps:= proc() local r,a,D,d,b,c; r:=args[1];a:=args[2]; c:=[-r,op(a)];D:=Digits;if nargs>2 then D:=args[3];fi; d := IntegerRelations[PSLQ](c); print(d,"Error is", evalf[D](linalg[dotprod](d,c)), "checking to",D,`places`); if d[1]=0 then print("Dependency detected");dp(d,c)=0 else r= linalg[dotprod](d[2 .. nops(d)]/d[1],c[2..nops(d)]) fi; end: Pslq:= proc() local r,a,D,d,b,c; r:=args[1];a:=args[2]; c:=[-r,op(a)];D:=Digits;if nargs>2 then D:=args[3];fi; d := IntegerRelations[PSLQ](c); print(d,"Error is", evalf[D+10](linalg[dotprod](d,c)), "checking to",D+10,`places`); if d[1]=0 then print("Dependency detected");dp(d,c)=0 else r= linalg[dotprod](d[2 .. nops(d)]/d[1],c[2..nops(d)]) fi; end: pslq:= proc() local r,a,D,d,b,c; r:=args[1];a:=args[2]; c:=map(Re,[-r,op(a)]);D:=Digits;if nargs>2 then D:=args[3];fi; d := IntegerRelations[PSLQ](c); print(d,"Error is", evalf[D+10](linalg[dotprod](d,c)), "checking to",D+10,`places`); r= linalg[dotprod](d[2 .. nops(d)]/d[1],c[2..nops(d)]) end: IND:=proc() local P,Q,eps,k,D; if nargs>1 then Digits:=args[2]; fi;eps:=10.^(Digits-3);P:=args[1]:Q:=IntegerRelations[PSLQ](P);print(nops(Q)); D:=add(Q[k]*P[k],k=1..nops(P));print(evalf(D)); if evalf(D)>eps then print('independent');P; elif Q[nops(P)]=0 then ListTools[Rotate](P,1) else [seq(P[n],n=1..nops(P)-1)];fi;end: test:=proc(x,d) Pslq(x,[op(x)],d);end: with(IntegerRelations): with(PolynomialTools): with(gfun):with(IntegrationTools); with(MmaTranslator): m2m:=proc(s); FromMma(s);end: _EnvExplicit:=true; with(combinat): with(IntegerRelations): with(PolynomialTools): with(student): with(gfun): print("1. combinat,gfun,student,IntegerRelations,PolynomialTools loaded"); print("2. CF,p2s,s2p,r2p,f2p,pslq,find loaded"); # bzeta(p,a) # Input: p numeric # a list of numerics # Output: # sum(1/n^p/binomial(2*n,n), product(sum(1/j^a[k],j=1..n-1),k=1..L), # n=2..infinity) `evalf/bzeta` := proc(p::numeric,a::list) local dig,os,b,j,s,i,n,S,c,t; dig := Digits; Digits := Digits + length(Digits+9)+4;b:=map(signum,a); c := 2; n := nops(a); if nops(b)<>n then ERROR(`lists must have the same number of elements`); fi; os := -1.; if b=[] or {op(b)}={0} then S := .5 else S := 0; fi; for j to n do s[j] := 0 od; for i while os<>S do os := S; t := 1.; for j to n do s[j] := s[j] + 1./i^a[j]; t := t*s[j]^b[j]; od; c := iquo((4*i+2)*c,i+1); S := S+t/c/(i+1)^p; od; Digits := dig; evalf(S); end: # binzeta(p,a,b) # Input: p numeric # a list of numerics # b list of numerics (of same length as a) # Output: # sum(1/n^p/binomial(2*n,n), product(sum(1/j^a[k],j=1..n-1)^b[k],k=1..L), # n=2..infinity) `evalf/binzeta` := proc(p::numeric,a::list,b::list) local dig,os,j,s,i,n,S,c,t; dig := Digits; Digits := Digits + length(Digits+9)+4; c := 2; n := nops(a); if nops(b)<>n then ERROR(`lists must have the same number of elements`); fi; os := -1.; if b=[] or {op(b)}={0} then S := .5 else S := 0; fi; for j to n do s[j] := 0 od; for i while os<>S do os := S; t := 1.; for j to n do s[j] := s[j] + 1./i^a[j]; t := t*s[j]^b[j]; od; c := iquo((4*i+2)*c,i+1); S := S+t/c/(i+1)^p; od; Digits := dig; evalf(S); end: # azeta(p,a) # Input: p numeric # a list of numerics # Output: # sum((-1)^n/n^p/binomial(2*n,n), product(sum(1/j^a[k],j=1..n-1),k=1..L), # n=2..infinity) `evalf/azeta` := proc(p::numeric,a::list) local dig,os,j,s,i,n,S,c,t,z,b; dig := Digits; Digits := Digits + length(Digits+9)+4; c := 2; n := nops(a);b:=map(signum,a); if nops(b)<>n then ERROR(`lists must have the same number of elements`); fi; z := 1; os := -1.; if b=[] or {op(b)}={0} then S := .5 else S := 0; fi; for j to n do s[j] := 0 od; for i while os<>S do os := S; t := 1.; for j to n do s[j] := s[j] + 1./i^a[j]; t := t*s[j]^b[j]; od; c := iquo((4*i+2)*c,i+1); z := -z; S := S+z*t/c/(i+1)^p; od; Digits := dig; evalf(S); end: # abinzeta(p,a,b) # Input: p numeric # a list of numerics # b list of numerics (of same length as a) # Output: # sum((-1)^(n+1)/n^p/binomial(2*n,n), # product(sum(1/j^a[k],j=1..n-1)^b[k],k=1..L), n=2..infinity) `evalf/abinzeta` := proc(p::numeric,a::list,b::list) local dig,os,j,s,i,n,S,c,t,z; dig := Digits; Digits := Digits + length(Digits+9)+4; c := 2; n := nops(a); if nops(b)<>n then ERROR(`lists must have the same number of elements`); fi; z := 1; os := -1.; if b=[] or {op(b)}={0} then S := .5 else S := 0; fi; for j to n do s[j] := 0 od; for i while os<>S do os := S; t := 1.; for j to n do s[j] := s[j] + 1./i^a[j]; t := t*s[j]^b[j]; od; c := iquo((4*i+2)*c,i+1); z := -z; S := S+z*t/c/(i+1)^p; od; Digits := dig; evalf(S); end: # bbinzeta(p,a,b) # Input: p numeric # a list of numerics # b list of numerics (of same length as a) # Output: # sum(1/(2*n+1)^p/binomial(2*n,n), #product(sum(1/(2*j+1)^a[k],j=0..n-1)^b[k],k=1..L), # n=1..infinity) `evalf/bbinzeta` := proc(p::numeric,a::list,b::list) local dig,os,j,s,i,n,S,c,t; dig := Digits; Digits := Digits + length(Digits+9)+4; c := 1; n := nops(a); if nops(b)<>n then ERROR(`lists must have the same number of elements`); fi; os := -1.; S:=0; #if b=[] or {op(b)}={0} then S := 1. else S := 0; fi; for j to n do s[j] := 0 od; for i from 0 while os<>S do os := S; t := 1.; for j to n do s[j] := s[j] + 1./(2*i+1)^a[j]; t := t*s[j]^b[j]; od; S := S+t/c/(2*i+1)^p; c := iquo((4*i+2)*c,i+1); od; Digits := dig; evalf(S); end: ind := proc (a, b, n,d) local t, s,S, k, j; global x,y; Digits:=d+10: s := evalf[d+10]([seq(seq(b^j*a^k, k = 0 .. n), j = 0 .. n)]); S:= [seq(seq(y^j*x^k, k = 0 .. n), j = 0 .. n)]; Digits:=d:t := IntegerRelations:-PSLQ(s); print(`error`,evalf[d+10]((add(t[k]*s[k], k = 1 .. nops(s))) = 0 ) ,`polynomial:`); add(t[k]*S[k], k = 1 .. nops(s)) = 0 end proc: alias(Az=abinzeta):alias(Sz=binzeta):alias(kappa=bzeta):alias(delta=azeta): ############################################################################# Col3:=proc(n) if irem(n,3)=0 then blue elif irem(n,3)=1 then red else black fi;end: anim:=proc(num,a,b,c,d) local P,k;for k from c to d do P[k]:=pointplot({seq([n,cos(k/10*n^(num))],n=a..b)},symbol=circle,symbolsize=1,colour=Col3(k),caption = typeset("Values of ",cos(k/10*n^num),"."),font = [Times, Roman, 12]):od;display([seq(P[k],k=c..d)],insequence=true);end: print("3. anim and binzetas loaded"); ####################################################