
# !pip install dependent_bterms tqdm


import dependent_bterms as dbt


%display typeset


N = 10000

var('t m j')
assume(m > 0, j+1 > 0, (j, 'integer'))


@parallel(ncpus=16)
def F(n):
    return sum(k*sigma(k)*(k**2 - 3*n + 2)*(2*k**2 - n)*binomial(2*n, n-k) for k in srange(1, n + 1))


%%time

from tqdm.notebook import tqdm

all(
    result < 0
    for (args, result) in tqdm(F(srange(5, N+1)), total=int(N-5))
)

  0%|          | 0/9995 [00:00<?, ?it/s]

CPU times: user 5.23 s, sys: 54.6 s, total: 59.8 s
Wall time: 11min 43s


CC = ComplexBallField(100)


# bound for sigma(n)

A = ceil(100 * (exp(euler_gamma) + 0.6483 / log(log(N))))/100
A, A.n()


all(sigma(kk) <= A*kk*log(log(N)) for kk in srange(1, N+1))


tail_error_end = A * log(log(m)) * m^7 * exp(-m/4)
tail_error_end


plot(tail_error_end, m, N/2, 2*N) + plot(m^2/8 - m/24, m, N/2, 2*N, color="red")


tail_error_end(m=CC(N))


alpha = 7/10
assert N^(2*alpha - 1) >= 3  # requirement in Section 4.2

AR, n, k = dbt.AsymptoticRingWithDependentVariable('n^QQ', 'k', 0, alpha, bterm_round_to=4, default_prec=7)


tail_error_mid = AR.B(4*A * (
        n^(6*alpha)
        + n/(2*n^alpha) * (6*n^3 + 6*n^(2 + 2*alpha)
                           + 3*n^(1 + 4*alpha) + n^(6*alpha))
    ),
    valid_from=N
)
tail_error_mid  # times log(log(n)) * exp(-n^(2*alpha - 1))

/home/behackl/git/papers/bona-dejonge-conjecture/code/dependent_bterms/dependent_bterms/structures.py:381: FutureWarning: This class/method/function is marked as experimental. It, its functionality or its interface might change without a formal deprecation.
See https://github.com/sagemath/sage/issues/31922 for details.
  super().__init__(


tail_error_mid_bterm = next(tail_error_mid.summands.elements())
tail_error_mid_sym = AR.create_summand(
    'exact',
    coefficient=tail_error_mid_bterm.coefficient,
    growth=tail_error_mid_bterm.growth
).subs(n=m) * log(log(m)) * exp(-m^(2*alpha - 1))
tail_error_mid_sym


plot(tail_error_mid_sym, m, N/2, 2*N) + plot(m^2/8 - m/24, m, N/2, 2*N, color="red")


tail_error_mid_sym(m=CC(N))


R = 9


sum_error = dbt.taylor_with_explicit_error(lambda t: 1/(1 - t), k^2/n^2, valid_from=N) * AR.B(2 / (R * n^R) * (k^R + k^(R+1)/(R+1)), valid_from=N)
sum_error


sum_expansion = -2*sum([sum(j^r, j, 1, k) / (r * n^r) for r in srange(1, R, 2)])
sum_expansion += k^2/n  # remove the asymptotic main term, it remains an exponential
sum_expansion += sum_error
sum_expansion


# check magnitude of B-term
dbt.simplify_expansion(sum_expansion, simplify_bterm_growth=True)


%%time
# need to find correct order of Taylor expansion. check collapsed B-term growth after calculations!

S_nk_and_error = dbt.taylor_with_explicit_error(exp, sum_expansion, order=6, valid_from=N)
S_nk_and_error *= dbt.taylor_with_explicit_error(lambda t: 1/(1-t), k/n, order=6, valid_from=N)
S_nk_and_error = dbt.simplify_expansion(S_nk_and_error)
S_nk_and_error

CPU times: user 4min 7s, sys: 1.88 s, total: 4min 9s
Wall time: 3min 52s


S_nk_and_error_collapsed = dbt.simplify_expansion(S_nk_and_error, simplify_bterm_growth=True)
S_nk_and_error_collapsed


S_nk = AR.zero()
for summand in S_nk_and_error.summands:
    if summand.is_exact():
        S_nk += AR(summand)
S_nk


expansion_error = dbt.simplify_expansion(S_nk_and_error.error_part() * (2*k^4 + 3*n^2) * k^2)
expansion_error


# we determine the integral bound symbolically:

f(t) = t^j * exp(-t^2/m)
t_0 = sqrt(j * m / 2)

integral_bound = f(t_0) + integrate(f(t), t, 0, oo)
integral_bound


error_after_summation = AR.zero()

for summand in expansion_error.summands:
    with assuming(k > 0):
        coef = summand.coefficient.simplify()
    for c, p in coef.coefficients(k):
        error_after_summation += c * AR.B(dbt.evaluate(integral_bound, m=n, j=p) * n^(summand.growth.exponent), valid_from=N)

error_after_summation *= A
error_after_summation # times log(log(n)).


error_after_summation_bterm = next(error_after_summation.summands.elements())
error_after_summation_sym = AR.create_summand(
    'exact',
    coefficient=error_after_summation_bterm.coefficient,
    growth=error_after_summation_bterm.growth
).subs(n=m) * log(log(m))
error_after_summation_sym


# not negligible, but still clearly dominated by main term.
plot(error_after_summation_sym, m, N/2, 2*N) + plot(m^2/8 - m/24, m, N/2, 2*N, color="red")


# the lazy way of plugging in k = n^(3/4): switch to
# the corresponding asymptotic ring.
AR_34, n_34, k = dbt.AsymptoticRingWithDependentVariable('n^QQ', 'k', 0, 3/4, bterm_round_to=3, default_prec=10)


AR_34(S_nk).O()  # good!


C_2 = dbt.expansion_upper_bound(AR_34(S_nk), numeric=True)  # <-- C2
C_2, C_2.n()


# symbolic integration, then back to produce a B-term.
tails_error_until_n34 = dbt.evaluate((
    m^(3/4) * (f(m^(alpha))(j=5) + integrate(abs(f(t)(j=5)), t, m^alpha, m^(3/4)))
    ).subs({
        exp(-sqrt(m)): -1,  # this is a bit delicate, but correct.
        exp(-m^(2/5)): 1,
    }),
    m=n
).B(valid_from=N)
tails_error_until_n34 *= 2*A
tails_error_until_n34 # times exp(-n^(2/5)) * log(log(n))


tails_error_until_n34_bterm = next(tails_error_until_n34.summands.elements())
tails_error_until_n34_sym = AR_34.create_summand(
    'exact',
    coefficient=tails_error_until_n34_bterm.coefficient,
    growth=tails_error_until_n34_bterm.growth,
).subs(n=m) * log(log(m)) * exp(-m^(2/5))
tails_error_until_n34_sym


# this error is negligible too!
(
    plot(tails_error_until_n34_sym, m, N/2, 2*N) 
    + plot(m^2/8 - m/24, m, N/2, 2*N, color="red")
)


AR(dbt.expansion_upper_bound(
    dbt.evaluate((S_nk.subs(n=m) * m^15 / k^20).expand()(k=m^3/4), m=n),
    numeric=True,
    valid_from=N
)).B(valid_from=N)


assert N^(3/4) > sqrt(27*N/2)


tails_error_after_n34 = (
    2 * m^(-15) / 10000 * (f(m^(3/4))(j=27) + integrate(f(t)(j=27), t, m^(3/4), oo))
).expand()
tails_error_after_n34 = dbt.evaluate((tails_error_after_n34 * exp(sqrt(m))).expand(), m=n).B(valid_from=N)
tails_error_after_n34  # times exp(-n^(1/2))


tails_error_after_n34_bterm = next(tails_error_after_n34.summands.elements())
tails_error_after_n34_sym = AR.create_summand(
    'exact',
    coefficient=tails_error_after_n34_bterm.coefficient,
    growth=tails_error_after_n34_bterm.growth,
).subs(n=m) * exp(-m^(1/2))
tails_error_after_n34_sym


# this error is negligible too!
(
    plot(tails_error_after_n34_sym, m, N/2, 2*N) 
    + plot(m^2/8 - m/24, m, N/2, 2*N, color="red")
)


mellin_summands = (k * S_nk * (k^2 - 3*n + 2) * (2*k^2 - n)).map_coefficients(lambda t: t.expand())
mellin_summands


var('t s')

def mellin_transform_from_summand(k_exp, n_exp):
    a = k_exp
    b = n_exp
    return zeta(2*s - 2*b - a - 1) * zeta(2*s - 2*b - a) * gamma(s - b) * t^(-s)


mellin_transform = 0

for summand in mellin_summands.summands:
    for coef, k_pow in summand.coefficient.coefficients(k):
        mellin_transform += coef * mellin_transform_from_summand(k_exp=k_pow, n_exp=summand.growth.exponent)

mellin_transform


# there are no residues at half-integers >= 5/2:
all(mellin_transform.residue(s==nm) == 0 for nm in srange(5/2, 30, 1/2))


main_asy = sum([
    mellin_transform.residue(s==res)
    for res in [2, 3/2, 1]
])
main_asy = dbt.evaluate(main_asy, t=1/n)
main_asy


zeta_fct = zeta(s).operator()
gamma_fct = gamma(s).operator()

def zeta_bound(sgm):
    if sgm >= 3/2:
        return ceil(zeta(3/2)*10^5)/10^5
    if sgm <= -1/2:
        return (
            ceil(
                2^(sgm + 1/2) * pi^(sgm - 1/2) * zeta(3/2)
                * exp(1/600) * (abs(1 - sgm - I*100)/100)^(1/2 - sgm)
                * 10^5
            ) / 10^5
            * t^(1/2 - sgm)
        )
    if sgm == 1/2:
        return 618/1000 * t^(1/2)
    
def gamma_bound(sgm):
    if sgm > 0:
        return (
            ceil(
                (2*pi)^(1/2) * exp(1/600)
                * (abs(sgm + 100*I)/100)^(max(sgm-1/2, 0))
                * 10^5
            ) / 10^5
            * t^(max(sgm-1/2, 0)) * exp(-pi/2 * t)
        )
    if sgm < 0:
        offset = 0
        while sgm + offset < 0:
            offset += 1
        return t^(-offset) * gamma_bound(sgm + offset)


%%time

integral_error = 0

for smd in mellin_transform.expand().operands():
    offset = 0
    sigma = 3/4
    while smd.residue(s==sigma - offset - 1/4) == 0:
        offset += 1/2  # shift as far to the left as long as there are no new residues
        if offset >= 100:  # some summands actually have no residues, trivial zeros of zeta cancel them.
            offset = 0
            break
    smd *= t^s
    smd = smd(s=sigma - offset + I*t)
    integral_contribution = CC.integral(
        lambda x, _: fast_callable(abs(smd), vars=[t])(x),
        -100, 100  # integral over Re(s) = sigma - offset, central part
    )
    
    # estimates for integral where |t| >= 100
    error_bound = 1
    for fct in smd.operands():
        if fct.operator() is None:
            error_bound *= abs(fct)
        if fct.operator() is not None:
            [arg] = fct.operands()
            if fct.operator() is zeta_fct:
                error_bound *= zeta_bound(arg(t=0))
            if fct.operator() is gamma_fct:
                error_bound *= gamma_bound(arg(t=0))
    [[err_t_coef, err_t_pow]] = error_bound.coefficients(t)
    integral_error_bound = 2 * CC(
        err_t_coef(t=0) 
        * CC((2/pi)^(err_t_pow + 1)) 
        * CC(err_t_pow + 1).gamma_inc(CC(50*pi))
    )
    integral_contribution += integral_error_bound
    integral_error += AR.B(abs(integral_contribution).upper() * n^(sigma - offset), valid_from=N)

integral_error /= 2*pi
integral_error

CPU times: user 1min 58s, sys: 21.4 ms, total: 1min 58s
Wall time: 1min 58s


integral_error_bterm = next(integral_error.summands.elements())
integral_error_sym = AR.create_summand(
    'exact',
    coefficient=integral_error_bterm.coefficient,
    growth=integral_error_bterm.growth,
).subs(n=m)
integral_error_sym


# again, the error is not negligible. it is still well dominated by the main term, however.
plot(integral_error_sym, m, N/2, 2*N) + plot(m^2/8 - m/24, m, N/2, 2*N, color='red')


# ratio of error term to main term at n=N
(integral_error_sym(m=N) / (N^2/8 - N/24)).n()


errors = [
    tail_error_end,
    tail_error_mid_sym,
    error_after_summation_sym,
    tails_error_until_n34_sym,
    tails_error_after_n34_sym,
    integral_error_sym,
]
total_error = sum(errors)
total_error


# crude estimates for n >= N:
assert log(log(N)) <= N^(1/10)
assert N > ((71/10) / (1/4)) # n^(71/10) * exp(-n/4) is decreasing
assert N > (97/20 / (2/5))^(5/2)  # n^(97/20) * exp(-n^(2/5)) is decreasing
assert N > (23/5 / (2/5))^(5/2)  # n^(23/5) * exp(-n^(2/5)) is decreasing
assert N > (11/2 / (1/2))^(2/1)  # n^(11/2) * exp(-n^(1/2)) is decreasing

total_error_estimate = 0
for err in errors:
    err = err.subs({log(log(m)): m^(1/10)})
    if exp in [op.operator() for op in err.operands()]:
        err = err(m=N)
        total_error_estimate += AR(err)
    else:
        total_error_estimate += dbt.evaluate(err, m=n)

total_error_estimate.B(valid_from=N)


bool(38755553/5000 * N^(3/4) < N^2/8 - N/24)


((38755553/5000 * N^(3/4)) / (N^2/8 - N/24)).n()


plot(total_error, m, N/2, 2*N) + plot(m^2/8 - m/24, m, N/2, 2*N, color='red')


plot(total_error / (m^2/8 - m/24), m, N/2, 2*N) + plot(1, m, N/2, 2*N, linestyle='dashed', color='black')

Example: On a conjecture of Bóna and DeJonge¶

Manual check for $n \leq N$¶

Rigorous asymptotic error bounds¶

Tail error, summation range $k \geq n/2$¶

Tail error, summation range $k\in [n^{\alpha}, n/2]$.¶

Summand approximations¶

Completion of the sum: error for completing the sum over $n^{\alpha} \leq k < n^{3/4}$¶

Completion of the sum: error for completing the sum over $k \geq n^{3/4}$¶

Mellin transform¶

Integral errors¶

Accumulated Error¶