4413.sagews

f(x) = (3/x^3*(sin(x) - x*cos(x)))^2

def I(y):
    return numerical_integral( (1 - f(x*sqrt(y)))/(2*cosh(x/2)^2), 0, Infinity)[0] \
            + numerical_integral( (1 - f(x*sqrt(y)))/sinh(x), 0, Infinity)[0]
            
def C(Npp,y):
    return 4*log(1.*Npp)*sum([ 1/(1+Npp^j)*f(j*log(1.*Npp)*sqrt(1.*y)) for j in range(1,100)])

plot(I, (0.1,100))

def poitou_bound(n, y, S):
    if y == 0:
        return 0
    else:
        return exp(1 + RR(euler_gamma) + log(4*RR(pi))\
                     - 12*RR(pi)/(5*n*sqrt(1.*y)) - I(y) + 1/n*sum([C(q,y) for q in S]))

def best_poitou_bound(n, S, y_start = 0.5):
    y = y_start
    y_step = y_start/2
    D = poitou_bound(n, y, S)
    while y_step > 10^(-6):
        D_lower = poitou_bound(n, y - y_step, S)
        D_upper = poitou_bound(n, y + y_step, S)
        if D_lower > D:
            D = D_lower
            y = y - y_step
        else:
            D = D_upper
            y = y + y_step
        y_step /= 2
    return D, y

best_poitou_bound(8, [2])

(12.6276716475697, 0.875001907348633)

12.6276716475697^8

6.46526372054113e8

best_poitou_bound(8, [2,2,2,2])

(22.1250933743379, 0.632814407348633)

18.23468^8

1.22232071052188e10