This Maple file accompanies the paper "The Tu-Deng Conjecture holds almost surely" from Lukas Spiegelhofer and Michael Wallner.#############################################################################################################################Abstract:The Tu-Deng Conjecture is concerned with the sum of digits of n in base 2 (the Hamming weight w(n) of the binary expansion of n) and states the following:assume that k is a positive integer and 0 < t < 2^k-1. Then|{ (a,b) in {0,...,2^k-2}^2 : a+b = t mod 2^k-1, w(a)+w(b)<k }| < 2^{k-1}+1.We prove that the Tu-Deng conjecture holds almost surely in the following sense: The proportion of t in [1,2^k-2] is such that the above inequality holds approaches one as k tends to infinity.Moreover, we prove that the Tu-Deng conjecture implies a conjecture due to T. W. Cusick concerning the sum of digits of n and n+t.#############################################################################################################################In this worksheet we compute the asymptotic expansions of certain generating function associated with the expected value and the variance of a random variable. These results are then used to prove the asymptotic version of the Tu-Deng Conjecture by showing concentration of the random variable with Chebyshev's inequality.For the computuations the gfun package needs to be loaded (mostly for the gdev command). For the latest version see Salvy's homepage http://perso.ens-lyon.fr/bruno.salvy/software/the-gfun-package/. For this computations gfun verison 3.76 was used.
Additionally, we use the gdev command part of the gdev package and available e.g., in the algolib package: http://algo.inria.fr/libraries/.
(Note that if you load the algolib package, the gfun package must be loaded first, because algolib contains a possibly older version of gfun.)restart;
#libname := "<path>\134\134gfun.mla", "<path>\134\134algolib.mla", libname;
with(gfun): gfun:-version;Proof of Proposition 4.2 - The asymptotics of the first momentWe want to compute [x^n y^(n-1)] G(x,y)G:=(2-4*x*y-x-x*y^2)/((1-y)*(1-2*x*y)*(1-x*(1+y)^2));
map(simplify,convert(gdev(convert(gdev(G,x,6),polynom),y,6),polynom));
add(coeff(coeff(%,x,n),y,n-1)*x^n*y^(n-1),n=1..5);First, we shift the function G by y in order to be able to compute the diagonal [x^n y^n] G(x,y) and not [x^n y^(n-1)] G(x,y)Then, we substitute one variable with t and the other one with x/t. The constant term of this function is the desired diagonal.
For more information on this method see Stanleys book "Enumerative Combinatorics Volume 2", Chapter 6.3.
This is the only choice for the substitution such that roots of denom are power series in x.
The leading coefficient of the denominator with respect to t is (-1+x).GT:=simplify(subs({x=t,y=x/t},y*G));
den:=[solve(denom(%),t)];
map(gdev,den,x,2);
We perform a partial fraction decomposition with respect to t.The variables A,B, and C are constants in x but independent of t.Note that we omitted the factor 1/(2*x-1) here.Note that at the end we are interested in extracting the coefficient of the constant term wrt to t. This is the diagonal of the initial BGF.Furthermore, our solution is a power series in x. Thus, in the following Ansatz only the second term with a B gives a function with nonnegative terms in t. In the first one, in order to ensure a power series in x, we have only negative powers of t. The same holds for the last term. The reason for this is the lack of a constant term in the expansions for den[1] and den[3] above. Yet den[2] has a constant term 2. Does it has nonnegative coefficients in t.In other words, we are at the end only interested in the second term involving B.ansatz := A/(t-den[1]) + B/(t-den[2]) + C/(t-den[3]);We solve for the unknowns A,B, and C.simplify(GT)*(2*x-1) = ansatz;
GTeq:=numer(GT) = (A*(t-den[2])*(t-den[3]) + B*(t-den[1])*(t-den[3]) + C*(t-den[1])*(t-den[2])):
ABCsubs:= map(simplify,solve({seq(map(coeff,GTeq,t,i),i=0..2)},{A,B,C}));
map(gdev,subs(ABCsubs,[A,B,C]),x=0,2);Inserting t=0 into the part involving B (which is legitimate as there are only non-negative powers of t) gives the diagonal generating function.
Note that we need to add the earlier omitted factor 1/(2*x-1).
In order to simplify it we compute an algebraic equation for it which we then solve. We have to choose the solution with a zero constant term.subs(ABCsubs,1/(2*x-1)*B/(-den[2]));
gdev(%,x,8);
algfuntoalgeq(%%,y(x));
solve(%,y);This is the final solution of the diagonal of y*G.GGx := (1/2)*(1/(1-4*x) - 1/sqrt(1-4*x));
gdev(GGx,x,6);We can however directly extract the coefficients in a closed form due tosum(binomial(2*n,n)*z^n,n=0..infinity);
series(sum(binomial(2*n,n)*z^n,n=0..10),z,6);
gdev((1-4*z)^(-1/2),z,6);This givesGN := n-> 1/2*(4^n-binomial(2*n,n));
seq(GN(n),n=1..10);And the asymptotics is directly computed from the known asyptotics of the central binomial coefficient. Note that we can compute it to any chose order.4^n/2*(1-map(simplify,asympt(binomial(2*n,n)/4^n,n,4))) assuming n::posint;Proof of Proposition 4.4 - The asymptotics of the second momentyzsubs := {y = 1 + I*s/sqrt(n), z = 1 + I*t/sqrt(n) };B0 := 1-x*y^2*z^2/(1-4*x*y*z)*(1 + 2*x*y/(1-2*x*y*(1+y*z)) + 2*x*z/(1-2*x*z*(1+y*z)));
B1 := 1/(1-2*x*z*(1+y*z))*(1 - x*(1+y*z)^2 - x*y*z/(1-2*x*y*(1+y*z)) - x*y*z);
B2 := 1/(1-2*x*y*(1+y*z))*(1 - x*(1+y*z)^2 - x*y*z/(1-2*x*z*(1+y*z)) - x*y*z);
B3 := 1-x/(1-4*x*y*z)*(1 + 2*x*y^2*z/(1-2*x*y*(1+y*z)) + 2*x*y*z^2/(1-2*x*z*(1+y*z)));We need [x^n y^(n-1) z^(n-1)] F(x,y,z) which is equivalent to [x^n y^n z^n] yzF(x,y,z) being the diagonal of a function.F := y*z*(B0+B1+B2+B3);We use the notation H for the function D(x,y,z) in the paper as D is the operator for differentiation in Maple.H := 1 - x*(1+y*z)^2 - x*y*z/(1-2*x*z*(1+y*z)) - x*y*z/(1-2*x*y*(1+y*z));Expand at (x,y,z) = (1/8,1,1)solve(subs([y=1,z=1],H),x);#degree of expansion
N := 8:
#expansion point
x0 := 1/8:
y0 := 1:
z0 := 1:
A := 'A':
Y := 'Y': Z := 'Z':
n := 'n': k := 'k':
#expansion point of x
xB := x0:
#shift y and z such that Y=y-y0, Z=z-z0
#then expansion at Y=0,Z=0
Hs := subs([y=Y+y0,z=Z+z0],H):
#degree of expansion
for n from 1 to N do
#all possibilites for coefficients of Y=(y-1) and Z=(z-1)
for k from 0 to n do
#new ansatz
xB := xB + A*Y^k*Z^(n-k);
#substitute ansatz into implicit equation H
impSubs := subs(x=xB, Hs);
#distinguish cases (problems with diff if 0 times differentiated)
if k = 0 then
impDif := diff(impSubs,Z$(n-k));
else if n = k then
impDif := diff(impSubs,Y$k);
else
impDif := diff(diff(impSubs,Y$k),Z$(n-k));
end if; end if;
#solve for variable A
tmp := [solve(eval(impDif,[Y=0,Z=0])=0,A)];
if nops(tmp) > 1 then
printf("ERROR: At ansatz for (y-1)^%d*(z-1)^%d %d branches detected", k,n-k,nops(tmp));
break;
end if;
XB[n,k] := op(1,tmp);
xB := subs(A=XB[n,k],xB);
end do;
end do;
n := 'n': k := 'k':
printf("This is f(Y,Z) with Y=(y-1) and Z=(z-1):");
sort(xB,[Z,Y],ascending);logxB := sort(expand(convert(gdev(convert(gdev(log(xB),Y=0,N+1),polynom),Z=0,N+1),polynom)),[Z,Y],ascending);Substitute the values of Y and Zfyz := subs([Y=y-y0,Z=z-z0],xB):
logfyz := subs([Y=y-y0,Z=z-z0],logxB):logglambdayz := convert(gdev(convert(gdev(logfyz + log(y) + log(z), z=1,N+1),polynom),y=1,N+1),polynom);
subs(yzsubs,%):
exptmp := map(simplify@(x->-(n+1)*x),convert(gdev(%,n=infinity,N),polynom));
Extract the terms which are responsible for the exponential growth and the ones leading to a normal approximation.expgrowth := subs([s=0,t=0],convert(exptmp,polynom));
#domterms := expgrowth-(1/4)*t^2-(1/4)*s^2-I*c*t-I*c*s;
domterms := expgrowth-(1/4)*t^2-(1/4)*s^2;Compute the other factor of the integrand: -F/(y*z*H_x)
without y*z in denominator! (extract [x^(n-1) y^(n-1) z^(n-1)])mtaylor(subs([y=Y+1,z=Z+1],subs(x=fyz,-F/(diff(H,x)))),[Y=0,Z=0],N+1):
sort(%,[Z,Y],ascending);
expand(subs(yzsubs,subs([Y=y-1,Z=z-1],%))):
sort(%,n);
Pnst := convert(gdev(%,n=infinity,N-2),polynom);Expansion in n, careful: square-roots in nNORD := N:
XX := exptmp-domterms:
map(simplify,gdev(Pnst*(add(XX^k/factorial(k),k=0..NORD)),n=infinity,NORD-1)) assuming n::posint;
intfact := convert(%,polynom):terms := [op(intfact)]:The generic result of the needed integrals. The value just depends on the needed power (i.e., moment).
See Lemma 4.7.intval := proc(k::integer)
local ret;
if k mod 2 = 1 then
ret := 0;
else
ret := sqrt(Pi)*2*factorial(k)/factorial(k/2);
end if;
return ret;
end proc:There is a shift in k needed as we divide by k (or s, respectively)intvalshift := proc(k)
if k=0 then
return -Pi*I;
else if k >0 then
return intval(k-1);
else
printf("ERROR: negative k in intval!");
end if;
end if;
end proc:This is the term-wise computation of the integral.subexpasy := 0:
for tmp in terms do
tmp;
#TODO: not optimal, depends on structure of terms,
#better: extract coefficients with respect to n!!!
if not has(tmp,n) then
#case without n
prefac := 1;
ordn := 1;
stpart := tmp;
else
prefac := op(1,tmp);
ordn := op(3,tmp);
stpart := expand(op(2,tmp));
end if;
if not evalb(whattype(stpart)=`+`) then
liststpart := [stpart];
else
liststpart := [op(stpart)];
end if;
liststpart;
if nops(tmp) > 3 then
printf("ERROR: tmp has more then 3 terms!");
break
end if;
cur := 0;
for tt in liststpart do
sdeg := degree(tt,s);
tdeg := degree(tt,t);
cur := cur + intvalshift(sdeg)*intvalshift(tdeg)*tt/s^sdeg/t^tdeg;
end do:
sdeg,tdeg,cur,prefac,ordn;
cur := simplify(cur*prefac*ordn);
subexpasy := subexpasy + cur;
end do:
subexpasy := map(simplify@(x->x/(2*Pi*I)^2),subexpasy) assuming n::posint;The main result: The asymptotics of [x^(n) y^(n-1) z^(n-1)]F(x,y,z)asyn := simplify(exp(expgrowth))*(subexpasy) assuming n::posint;and a simplified expressionsimplify(exp(expgrowth)/8)*(map(x->8*x,subexpasy)) assuming n::posint;LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=