Let $f$ be a weight $k$ cusp form with Fourier expansion

$$ f(z) = \sum_{n \geq 1} a(n) e(nz). $$Deligne has shown that $a(n) \ll n^{\frac{k-1}{2} + \epsilon}$. It is conjectured that

$$ S_f^1(n) := \sum_{m \leq X} a(m) \ll X^{\frac{k-1}{2} + \frac{1}{4} + \epsilon}. $$It is known that this holds on average, and we recently showed that this holds on average in short intervals. (See HKLDW1, HKLDW2, and HKLDW3 for details and an overview of work in this area). This is particularly notable, as the resulting exponent is only 1/4 higher than that of a single coefficient. This indicates extreme cancellation, far more than what is implied merely by the signs of $a(n)$ being random.

It seems that we also have

$$ \sum_{m \leq X} S_f^1(m) \ll X^{\frac{k-1}{2} + \frac{2}{4} + \epsilon}. $$That is, the sum of sums seems to add in only an additional 1/4 exponent. This is unexpected and a bit mysterious.

The purpose of this notebook is to explore this and higher conjectures. Define the $j$th iterated sum as

$$ S_f^j(X) := \sum_{m \leq X} S_f^{j-1} (m).$$Then we numerically estimate bounds on the exponent $\delta(j)$ such that

$$ S_f^j(X) \ll X^{\frac{k-1}{2} + \delta(j) + \epsilon}. $$In [1]:

```
# This was written in SageMath 7.3 through a Jupyter Notebook.
# Jupyter interfaces to sage by loading it as an extension
%load_ext sage
# sage plays strangely with ipython. This re-allows inline plotting
from IPython.display import display, Image
```

We first need a list of coefficients of one (or more) cusp forms. For initial investigation, we begin with a list of 50,000 coefficients of the weight $12$ cusp form on $\text{SL}(2, \mathbb{Z})$, $\Delta(z)$, i.e. Ramanujan's delta function. We will use the data associated to the 50,000 coefficients for pictoral investigation as well.

We will be performing some numerical investigation as well. For this, we will use the first 2.5 million coefficients of $\Delta(z)$

In [2]:

```
# Gather 10 coefficients for simple checking
check_10 = delta_qexp(11).coefficients()
print check_10
fiftyk_coeffs = delta_qexp(50000).coefficients()
print fiftyk_coeffs[:10] # these match expected
twomil_coeffs = delta_qexp(2500000).coefficients()
print twomil_coeffs[:10] # these also match expected
```

In [3]:

```
# Function which iterates partial sums from a list of coefficients
def partial_sum(baselist):
ret_list = [baselist[0]]
for b in baselist[1:]:
ret_list.append(ret_list[-1] + b)
return ret_list
print check_10
print partial_sum(check_10) # Should be the partial sums
```

In [4]:

```
# Calculate the first 10 iterated partial sums
# We store them in a single list list, `sums_list`
# the zeroth elelemnt of the list is the array of initial coefficients
# the first element is the array of first partial sums, S_f(n)
# the second element is the array of second iterated partial sums, S_f^2(n)
fiftyk_sums_list = []
fiftyk_sums_list.append(fiftyk_coeffs) # zeroth index contains coefficients
for j in range(10): # jth index contains jth iterate
fiftyk_sums_list.append(partial_sum(fiftyk_sums_list[-1]))
print partial_sum(check_10)
print fiftyk_sums_list[1][:10] # should match above
twomil_sums_list = []
twomil_sums_list.append(twomil_coeffs) # zeroth index contains coefficients
for j in range(10): # jth index contains jth iterate
twomil_sums_list.append(partial_sum(twomil_sums_list[-1]))
print twomil_sums_list[1][:10] # should match above
```

*sizes* of the coefficients.
But in HKLDW3, we investigated some of the sign changes of the partial sums.

In [5]:

```
norm_list = []
for n,e in enumerate(fiftyk_coeffs, 1):
normalized_element = 1.0 * e / (1.0 * n**(5.5))
norm_list.append(normalized_element)
print norm_list[:10]
```

In [6]:

```
# Make a quick display
normed_coeffs_plot = scatter_plot(zip(range(1,60000), norm_list), markersize=.02)
normed_coeffs_plot.save("normed_coeffs_plot.png")
display(Image("normed_coeffs_plot.png"))
```

In [7]:

```
# We will be using matplotlib itself to make some of our plots.
# This import is necessary for us to interface between sage notebook and matplotlib.
from matplotlib.backends.backend_agg import FigureCanvasAgg
# We will use matplotlib's excellent pyplot
import matplotlib.pyplot as plt
# And numpy
import numpy as np
from scipy.optimize import curve_fit
```

In [8]:

```
fig = plt.figure(facecolor='0.97') # I use just off-white backgrounds
ax = plt.subplot(111) # create a plot
ax.scatter(range(1, len(norm_list)+1), # x coords
norm_list, # y coords
edgecolor='', # remove black outline of points
color='black', # make inner points black
s=1) # set pointsize to 1
ax.set_ylim([-2.5, 2.5]) # from quick figure above, most data is here
ax.set_xlim([-1000, 50000]) # seperate left border from data
ax.set_frame_on(False) # remove bounding box in plot
ax.xaxis.tick_bottom() # only have x-ticks and labels on bottom
ax.yaxis.tick_left() # only have y-ticks and labels on left
# Now let's save it
fig.set_canvas(FigureCanvasAgg(fig))
fig.savefig('quebecmaine_normalized.png') # lossful but fast
fig.savefig('quebecmaine_normalized.pdf') # lossless but heavy
display(Image("quebecmaine_normalized.png"))
# Doesn't that look better?
```