k := 41862295; m := 4; # These three parameters may be chosen arbitrarily c := evalf(0.9942920827819165/log(k)); T := evalf(0.8541193299758623/log(k)); # tau := 1.14/log(k); M0 := k*log(k)/(k-1); g := 1/(c+(k-1)*t); m2 := int(g*g, t=0..T ); mu := int(t*g*g, t=0..T )/m2; sigma2 := int(t*t*g*g, t=0..T)/m2 - mu*mu; tau := 1 - k*mu; denominator := (1+tau/2)*(1 - (k*sigma2)/((1+tau-k*mu)^2)); a := (r - k*mu) / T; S := r * (log((r-k*mu)/T) +(k*sigma2)/(4*a*a*T*T*log(a)) ) + r*r/(4*k*T); Sav := int(S, r=1..1+tau, numeric)/tau; Z3 := int( g*g*k*t*log(1+t/T), t=0..T, numeric) / m2; W := int( g*g*log(1+tau/(k*t)), t=0..T, numeric) / m2; X:= (log(k)/tau) * c^2; V := (c/m2) * int( g^2 / (2*c+(k-1)*t) , t=0..T, numeric); U := (log(k)/c) * int( (1+u*tau - (k-1)*mu - c)^2 + (k-1)*sigma2, u=0..1, numeric); numerator := Sav + Z3 + W*X + V*U; Delta := (k/(k-1)) * numerator/ denominator; # this needs to be positive evalf( 1 - k*mu - T ); evalf( 1 - k*mu - tau ); # this needs to be positive too evalf(denominator); # this is the lower bound for M_k evalf(M0-Delta); varpi := m/(M0-Delta) - 1/4; delta := T * ((1/4) + varpi); # if this is less than 1, we win evalf(1080*varpi/13+ 330*delta/13);