# Machine Learning amplitudes for faster event generation

Fady Bishara<EMAIL_ADDRESS>DESY, Notkestraße 85, 22607 Hamburg,
Germany  Marc Montull<EMAIL_ADDRESS>DESY, Notkestraße 85, 22607
Hamburg, Germany  Paul Scherrer Institut, Forschungsstraße 111, 5232 Villigen
PSI, Swizerland  Physik-Institut, Universität Zürich, Winterthurerstrasse 190,
CH-8057 Zürich, Switzerland

###### Abstract

We propose to replace the exact amplitudes used in MC event generators for
trained Machine Learning regressors, with the aim of speeding up the
evaluation of slow amplitudes. As a proof of concept, we study the process
$gg\to ZZ$ whose LO amplitude is loop induced. We show that gradient boosting
machines like XGBoost can predict the fully differential distributions with
errors below $0.1\%$, and with prediction times $\mathcal{O}(10^{3})$ faster
than the evaluation of the exact function. This is achieved with training
times $\sim 7$ minutes and regressors of size $\lesssim 30$ Mb. These results
suggest a possible new avenue to speed up MC event generators.

††preprint: DESY 19-232

##  I Introduction

The success of the LHC in discovering the Higgs boson is a testament to the
impressive advancements made by the HEP community in understanding
accelerators, detectors, and to make accurate Standard Model (SM) predictions.
As a result, the LHC is rapidly evolving from an energy frontier machine,
capable of discovering new resonances, to a precision machine, capable of
measuring small deviations over precise SM predictions.

Due to the key role of higher order corrections in precision physics, there
has been an Herculean effort in recent years to compute, store, and automate
higher loop calculations for SM and Beyond the SM (BSM) predictions Alwall
_et al._ (2014); Hirschi and Mattelaer (2015); Boughezal  _et al._ (2017);
Campbell and Neumann (2019); Nagy (2003, 2002); Gavin  _et al._ (2011);
Camarda  _et al._ (2019); Catani  _et al._ (2009, 2015); de Florian  _et al._
(2011); Actis  _et al._ (2017); Hirschi  _et al._ (2011); Li  _et al._ (2016);
Actis  _et al._ (2017); Naterop  _et al._ (2019); Peraro (2019); Camarda  _et
al._ (2019); Bern  _et al._ (2014); Bertone  _et al._ (2017); Maitre  _et al._
(2016); Spira (2017). As impressive as this has been, the use of N(N)LO
results by the broader HEP community has been relatively low, in part due to
the long times required to evaluate amplitudes beyond tree level. This
evaluation time increases dramatically with the loop order, and makes certain
Monte Carlo (MC) event simulations at one loop already unfeasible.
Nonetheless, higher loop effects will become more important as the precision
from the experimental and theoretical sides keeps improving. This calls for
innovations to reduce evaluation times for slow amplitudes. One possible
avenue to do just that is to improve the traditional tools and techniques – an
effort that is well under way. In this work, however, we take a new and
different approach to address these issues.

The main goal of this work is to show that thanks to the advances in Machine
Learning (ML) algorithms and tools, it is now possible to train ML regressors
with pre-computed slow amplitudes, and use them to predict the same amplitudes
accurately and in a fraction of the time.

As a proof of concept we study the $gg\to ZZ$ process which is loop induced at
LO (see Fig. 1). We find that ML regressors can achieve prediction times
$\mathcal{O}(10^{3})$ faster than traditional tools while the predicted values
for single and double differential distributions have errors below $0.1\%$.
This was achieved with training times $\lesssim$ 7 minutes on a single CPU
core, and with a disk size for the trained regressors of a few to tens of
megabytes (Mb).

Machine Learning algorithms are constantly finding new applications in HEP
research (see Bellagente  _et al._ (2019); Andreassen and Nachman (2019);
Nachman (2019); Bendavid (2017); Otten  _et al._ (2019); Klimek and Perelstein
(2018); Bothmann and Debbio (2019); Ball  _et al._ (2017, 2015); Santos  _et
al._ (2017); Farina  _et al._ (2018); D’Agnolo and Wulzer (2019); Blance  _et
al._ (2019); Cerri  _et al._ (2019); Brochet  _et al._ (2019); Banerjee  _et
al._ (2019a); Komiske  _et al._ (2017); Coccaro  _et al._ (2019); Andreassen
_et al._ (2019a, b); de Oliveira  _et al._ (2016); Metodiev  _et al._ (2017);
Collins  _et al._ (2018, 2019); Andreassen  _et al._ (2019c); Paganini  _et
al._ (2018a, b); Chakraborty  _et al._ (2020); Lim and Nojiri (2018) for
concrete applications and Luisoni  _et al._ (2016); Buckley (2019); Abdughani
_et al._ (2019); Albertsson  _et al._ (2018); Albrecht  _et al._ (2019);
Carrazza (2018); Larkoski  _et al._ (2017); Bourilkov (2019); Brehmer  _et
al._ (2019a, b); Carleo  _et al._ (2019) for recent reviews). Nonetheless, we
are not aware of any work where these tools have been used to speed up time
consuming amplitudes 111The only similar works we found is Ref. Otten  _et
al._ (2018) where DNN are used to predict NLO cross sections for the pMSSM-19
and Ref. Chawdhry  _et al._ (2019) ML regressors were used as interpolators to
store the NNLO QCD amplitude for $pp\to 3\gamma$..

There are many processes where speeding up MC event generation can be
immediately useful. Therefore, a next step after this work would be to test
the ML regressors on other processes and implement them into a MC generator.
We comment on further applications at the end of this letter.

Figure 1: SM LO diagrams for $gg\to ZZ$, up to fermion momentum flow and
crossings.

##  II A proof of concept with $gg\to ZZ$

We chose to test the performance of ML regressors in approximating the $gg\to
ZZ$ squared amplitude for several reasons. First of all, the LO contribution
to this process arises at one loop (Fig. 1) and therefore is slow. Secondly,
it was shown in Ref. Cascioli  _et al._ (2014) that this process contributes
the bulk ($\leavevmode\nobreak\ \sim 60\%$) of the full NNLO correction of
hadronic Z-boson pair production, making its computation imperative when
performing phenomenlogical studies to test the SM or to search for New Physics
(NP). In addition, it is relevant for NP searches where it constitutes a
background to $pp\to ZH$ with $H$ decaying to $\bar{b}b$ or to invisible new
particles Khachatryan  _et al._ (2017); Aaboud  _et al._ (2018); Sirunyan  _et
al._ (2018); Banerjee  _et al._ (2018, 2019b); Arcadi  _et al._ (2019). At the
same time, this process is simple enough to avoid unnecessary complications:
the squared amplitude only depends on two variables, the center of mass energy
and the polar angle $\theta$, i.e.
$|\mathcal{M}(\sqrt{\hat{s}},\cos\theta)|^{2}$. Furthermore, when the the pair
of $ZZ$ bosons are on-shell, there are no resonant peaks. We leave the study
of processes with $s$-channel resonances for future work. Furthermore, since
the $\alpha_{s}$ dependence amounts to an overall rescaling of the amplitude
squared, we can approximate the function using a fixed value of $\alpha_{s}$
and restore the scale dependence afterwards.

##  III ML Algorithm and training

###  III.1 Choosing a Machine Learning algorithm

The problem we are trying to address here requires, above all, two features
from an ML algorithm: first, it must be able to approximate the true function
over the entire domain as accurately as possible; second, it must be be able
to do so faster than existing dedicated programs $\sim 5\cdot 10^{-3}\;[s]$
per phase space (PS) point 222We determined this evaluation time using
OpenLoops Buccioni  _et al._ (2019) to generate the exact $gg\to ZZ$ squared
amplitudes for the training and prediction sets. The CPU core used for this
timing is the same one used for all timings in this letter, see main text.. An
additional bonus feature is for the model to be lightweight, i.e. to have a
small disk size, $\lesssim\mathcal{O}(100)$ Mb, so that it is easy to
distribute quickly.

With this in mind, we evaluated several algorithms suited for regression in
the early stages of this work. In particular, we tested deep neural networks
(DNN) with TensorFlow Abadi  _et al._ (2015), random forests Wright and
Ziegler (2017); Ho (1995); Breiman (2001), and gradient boosting machines
(GBM) Breiman (1997); Friedman (2001). From the outset, GBMs as implemented in
XGBoost Chen and Guestrin (2016) outperformed the others by far in terms of
speed, accuracy, and robustness against overfitting with very little tuning
333XGBoost rose to prominence by winning the Kaggle Higgs Callenge. Since then
it has been consistently on the top of the ladder in a large number of the
Kaggle competitions https://www.kaggle.com/sudalairajkumar/winning-solutions-
of-kaggle-competitions ; https://www.kaggle.com/rajiao/winning-solutions-of-
kaggle-competitions outperforming other ML architectures, like DNN’s or
SVM’s.. Therefore, all the results presented in this letter were obtained with
XGBoost via the scikit-learn API.

As discussed above, we use the default or close to the default values for the
hyper-parameters, except for the number of estimators ($n$), maximum depth of
the trees ($m_{d}$), and the learning rate ($l_{r}$), for which we performed a
small scan $n\in[10,1000]$, $m_{d}\in[10,800]$, and $l_{r}\in[0.01,0.3]$.
Based on this bare bones optimization, the final set of parameters used in
this work are given in Table 1.

XGBoost parameter |  Value
---|---
n_estimators |  200
max_depth |  50
learning rate |  0.1
min_child_weight |  1
$\gamma$ |  0
colsample_bytree |  1
subsample |  0.75
booster |  gbtree
objective |  reg:squarederror
Table 1: Hyper-parameter settings used for all XGBoost regressors in this
work. The parameters we attempted to optimize are shown above the split while
the values for the parameters below the split are the XGBoost default ones
with the exception of ‘subsample’, see text for details.

###  III.2 Datasets for Training and Prediction

To train and test the XGBoost regressor, we generated $18$ million (18M) pairs
of phase space points, $(\sqrt{\hat{s}},\cos\theta)$, uniformly distributed in
the region defined by,

$\begin{split}\sqrt{\hat{s}}&\in\left[2\sqrt{m_{Z}^{2}+p_{T,\text{cut}}^{2}},3\,\mathrm{TeV}\right]\,,\\\
\cos\theta&\in[-1,1]\times\sqrt{1-\frac{4p^{2}_{T,\text{cut}}}{\hat{s}-4m_{Z}^{2}}}\,,\end{split}$
(1)

with $p_{T}^{\text{cut}}=1$ GeV to regulate the singularity in
$\langle\left|\mathcal{M}\right|^{2}\rangle$ in the limit $p_{T}\rightarrow 0$
(similarly to what is done in MCFM Campbell and Ellis (2010) and
Madgraph_aMC@NLO Hirschi and Mattelaer (2015)). We chose
$(\sqrt{\hat{s}}\,)_{max}$ = 3 TeV as an arbitrary cutoff relevant for LHC
physics. Nevertheless, it is straightforward, and inconsequential, to extend
the cutoff to the collider center of mass energy; we checked this explicitly
up to 14 TeV.

We then computed the corresponding squared amplitudes required for training
and testing XGBoost , using the OpenLoops software 444We also tested MadLoops
Hirschi and Mattelaer (2015) and found it to be a factor of a few slower than
OpenLoops Buccioni  _et al._ (2019) for the $gg\to ZZ$ process. For this
reason, we used the faster evaluation time of OpenLoops as a benchmark.. The
full sample of $18$M points was split into a training and prediction dataset
with 3M and 15M points, respectively. To ensure that the data sets are
statistically independent, we generated the phase space points using the
Python implementation of the Mersenne Twister algorithm Matsumoto and
Nishimura (1998) which, when initialized properly, has a period of
$2^{19937}-1$.

###  III.3 Phase Space partitioning and multiple regressors

The function that we are trying to approximate,
$\langle\left|\mathcal{M}\right|^{2}\rangle$, is peaked at $\cos\theta\to\pm
1$. This motivates an ansatz to break up the full phase space into smaller
sub-regions with roughly equal integrated
$d\langle\left|\mathcal{M}\right|^{2}\rangle/d\mbox{PS}$ with the purpose of
training one regressor per sub-region.

For example, choosing $\cos\theta$ and $\sqrt{\hat{s}}$ regions defined by

$\displaystyle\cos\theta\in\\{-1,-0.94,-0.7,0.7,0.94,1\\}\times\sqrt{1-\frac{4p_{T,\text{cut}}^{2}}{\hat{s}-4m_{Z}^{2}}}\,,$
$\displaystyle\sqrt{\hat{s}}\in\\{2\sqrt{m_{Z}^{2}+p_{T,\text{cut}}^{2}},1.3\;\mbox{[TeV]},3\mbox{[TeV]}\\}\,,$
(2)

with the idea of decreasing as much as possible the variation of the squared
amplitude in each sub-region. This partitions the full phase space into 10
sub-regions each with its own dedicated XGBoost regressor that is trained on,
and predicts in, only one sub-region. These partitions are delineated by
dashed gray lines in the right half of Fig. 3.

For the remainder of this letter, we will refer to the ansatz with 10
regressors as the ‘10 regions’ regressor, and to the one trained on the full
domain defined by Eq. (1) as ‘1 region’ .

To compare the performance of the ‘1 region’ and ‘10 regions’ regressors, we
first train the ‘1 region’ regressor on a given dataset. Then, to train each
of the ten regressors that make up the ‘10 regions’ regressor, we split the
same dataset according to the regions defined in Eq. (2). In the end, each of
these ten regressors making the ‘10 regions’ regressor is only trained on a
fraction between 2.4% and 24% of the total dataset, corresponding to the
fraction of its phase space area (since the PS is uniformly sampled).

###  III.4 Training time

Figure 2: Prediction time per point as a function of the size of training set.
The crosses (open circles) correspond to the results for the ‘1 region’ (‘10
regions’ ) regressors. The solid curves are simple power law fits and are
shown in the legend.

We benchmark the time it took to train the ‘1 region’ and each of the ‘10
regions’ regressors on a single CPU core of an Intel® Xeon® CPU model
E5-2640V4 @ 2.40 GHz on x86_64 architecture. Since XGBoost can train and
predict on multiple cores by default, the times reported here are quite
conservative. In practice, modern desktop machines with at least four cores
are increasingly common and so training and prediction times can be easily be
improved by a factor of a few to ten.

The results of the timing tests for the training phase are shown in Fig. 2
(top panel) for both the ‘1 region’ and ‘10 regions’ regressors; in addition,
a simple power law fit to the points is given. For a training dataset size of
3M PS points, the ‘1 region’ (‘10 regions’ ) regressors took $\sim 16\,(7)$
minutes to train. In the case of the ‘10 regions’ regressors – there are 10 of
them – we added up the times it took to train each one of them.

##  IV Results

Figure 3: Absolute value of the percentage relative error per bin of the
double differential distribution,
$d^{2}\langle\left|\mathcal{M}\right|^{2}\rangle/d\cos\theta d\sqrt{\hat{s}}$.
Each bin has size $140\times 0.2$ (GeV, $\cos\theta$). The total number of
training (prediction) points is 3M (16M). Left: ‘1 region’ regressor. Right:
‘10 regions’ regressor with dashed gray lines showing the sub-regions defined
in Eq. (2).

In order to benchmark the trained ‘1 region’ and ‘10 regions’ regressors
defined above, we study the relative error of their predictions, and measure
their evaluation times. The relative error is defined as,

$\varepsilon=\frac{\langle\left|\mathcal{M}\right|^{2}\rangle_{\mathtt{OpenLoops}}-\langle\left|\mathcal{M}\right|^{2}\rangle_{\mathtt{XGBoost}}}{\langle\left|\mathcal{M}\right|^{2}\rangle_{\mathtt{OpenLoops}}}\,.$
(3)

As for the training times, we measure the prediction times on a single core of
the same CPU described above.

###  IV.1 Accuracy of predictions

Figure 3 shows the relative error on the sum of the predicted
$d^{2}\langle\left|\mathcal{M}\right|^{2}\rangle/d\mbox{PS}$ values in each
bin for the ‘1 region’ and ‘10 regions’ regressors. Each of the regressors was
trained on 3M points and the relative errors were computed from predictions of
15M points. Each bin has a size of 140 GeV $\times\,0.2$, which is appropriate
for phenomenological studies at the LHC. The left and right sides of the plot
correspond to the ‘1 region’ and ‘10 regions’ regressors, respectively, put
together for easier comparison since the amplitude is symmetric under
$\cos\theta\to-\cos\theta$. In addition, the right panel is overlaid with the
boundaries of the sub-regions defined in Eq. (2). We find that the ‘1 region’
regressor has a maximum relative error per bin of $0.3\%$ while the ‘10
regions’ regressor has a maximum error of $0.03\%$.

Figure 4: Relative error of
$d\langle\left|\mathcal{M}\right|^{2}\rangle/d\sqrt{\hat{s}}$. The red and the
blue curves correspond to the ‘1 region’ and ‘10 regions’ regressors
respectively (see Sec. III.3).

For phenomenological studies, another important differential distribution is
the singly differential one with respect to $\sqrt{\hat{s}}$. The relative
error for this distribution is shown in Fig. 4 and is of
$\mathcal{O}$(percent) and $\mathcal{O}$(permille) for the ‘1 region’ and ‘10
regions’ regions respectively.

In order to assess the effect of the size of the training set on the
performance of the machines, we show in Fig. 5 the fraction of points with
relative error greater than 1%, 5%, and 10% using the full 15M phase space
point prediction dataset. The dashed (solid) curves correspond to the full
(subdivided) phase space. Again, it is clear that subdividing the phase space
is very effective in reducing the errors. Figure 5 also shows that, for the
chosen hyperparameters (Table 1), there is little benefit from using training
datasets larger than 1M points. Furthermore, we find no over-training even up
to training datasets of 3M points.

Figure 5: Percentage of predicted points with an error greater than $1\%$
(blue), $5\%$ (red), and $10\%$ (black) as a function of the number of trained
points (the number of predicted points is 15 million). The solid (dashed)
curves correspond to the ‘10 regions’ (‘1 region’ ) regressors.

From the results presented in this section, the benefit of subdividing the
phase space and training separate machines on the subregions is clear: the
error between the ‘1 region’ and ‘10 regions’ is reduced by an order of
magnitude, see Fig. 3 and Fig. 4, while the training and prediction times are
reduced by a factor of two, see Fig. 2.

###  IV.2 Prediction speed

The time required to predict one phase space point is a crucial performance
metric for the trained machine. Clearly it must be much faster than the time
required to evaluate the true function (we use OpenLoops as our benchmark).
The results of the timing tests for the training phase are shown in Fig. 2 for
both the ‘1 region’ and ‘10 regions’ regressors. In addition, a simple power
law fit to the points is shown for each one on the plot. For the ‘10 regions’
regressor trained on 1M points, the prediction time is $\sim 1\times 10^{-5}$
seconds in comparison to $8.7\times 10^{-3}$ seconds for Fortran interface of
OpenLoops – i.e., a factor of $\sim 1000$ speedup.

Note that the trained regressors can be packaged as a single, standalone, C
library. We checked that calling this library during an event generation run
has negligible overhead.

###  IV.3 Disk size

Another desirable feature for the standalone packaged library is to be
lightweight in terms of disk size. We find that for 1M points, the ‘1 region’
(‘10 regions’ ) regressor has a size is 2.6 (19) Megabytes. This makes these
regressors ultra portable and could be downloaded on the fly during MC event
generation – we envision machines of this type to be an option given to the
user when generating events with popular MC event generators such as
MG5_aMCNLO Alwall  _et al._ (2014), Sherpa Bothmann  _et al._ (2019), and
Whizard Kilian  _et al._ (2018).

##  V Summary and conclusions

The idea of using ML regressors to approximate squared amplitudes proposed in
this work is a new application of Machine Learning techniques in HEP. Our goal
is to accurately predict the trained squared-amplitudes in a fraction of the
time it takes to evaluate the exact ones.

As a proof of concept, we studied the accuracy and speed of the XGBoost
regressor to predict the squared amplitudes for the $gg\to ZZ$ process which
at LO is generated at one loop. Our results show that the XGBoost regressors
deliver a 1000-fold speedup in evaluation time with respect to OpenLoops with
no more than $\mathbf{0.03\%}$ relative error with respect to the true double
differential distribution binned as in Fig. 3.

Another convenient feature of the XGBoost regressor studied in this letter, is
its reduced training speed. Using the hyper-parameters given Table 1, training
on 1M uniformly generated PS points takes about 2 minutes on one CPU core.
Moreover, since XGBoost is by default able to train and predict on multi-core
CPUs, actual training and prediction times will be in practice faster by a
factor proportional to the number of available CPU cores. XGBoost can also run
on GPUs with some minor modifications, otherwise, LightGBM Ke  _et al._ (2017)
works on GPUs by default and could even be a better performing regressor.

In addition, the disk size of the trained XGBoost regressor for this process,
is at most $\mathbf{30\;\mbox{Mb}}$, making it easy to distribute on the fly
during process generation in MC event generators.

Another important result of this work is to demonstrate that the errors on the
predictions of the XGBoost regressor can be reduced by an order of magnitude
by training independent regressors on separate sub-regions of the full phase
space. A bonus feature of training more regressors on sub-regions is that
their aggregate training and prediction times for a given dataset are reduced
with respect to training a single regressor on the full phase space

In Table 2 we summarize the aforementioned performance benchmarks for one
XGBoost regressor (‘1 region’ ) trained on the full phase space region, and
for ten regressors (‘10 regions’ ) each trained on a sub-region.

|  ‘1 region’ |  ‘10 regions’ |
---|---|---|---
$\left|\varepsilon^{\text{bins}}_{min}\right|\leavevmode\nobreak\ [\%]$ |  $7\cdot 10^{-5}$ |  $3\cdot 10^{-5}$ |  Fig. 3
$\left|\varepsilon^{\text{bins}}_{max}\right|\leavevmode\nobreak\ [\%]$ |  $0.3$ |  $0.03$ |  Fig. 3
$t_{\textrm{predict}}^{\text{(1 core)}}\leavevmode\nobreak\ [s/\textrm{point}]$ |  $2\cdot 10^{-5}$ |  $10^{-5}$ |  Fig. 2
$t_{\mathrm{train}}^{\text{(1 core)}}\leavevmode\nobreak\ [s]$ |  977 |  390 |  Fig. 2
$\mathrm{Size}\;[\mathrm{Mb}]$ |  $4.8$ |  28 |
Table 2: Main characteristics of the two ML regressors trained on 3M points
and predicting on 15M. The relative errors $\varepsilon^{\text{bins}}_{max}$
and $\varepsilon^{\text{bins}}_{min}$ stand for the relative errors in the
bins of size $140\leavevmode\nobreak\ \times\leavevmode\nobreak\ 0.2$
($\sqrt{\hat{s}}$ [GeV] $\times\,\cos\theta$).

The success of the proof of concept studied in this work suggests many
applications and further ideas to explore:

  * •

To test the performance of ML regressors on qualitatively different channels
with slow amplitudes. For instance: i) Amplitudes with resonant s-channels ii)
N(N)LO amplitudes iii) $2\to n$ processes.

  * •

Test and benchmark other ML algorithms.

  * •

Implement the trained ML regressors into an MC event generator.

  * •

On a side note, it would be interesting to test the performance of GBM on
interpolating PDFs and NNLO grids.

###### Acknowledgements.

#### Acknowledgements.

We are especially thankful to P. Englert for his contributions during the
initial stages of this project. We are happy to thank M. F. Zoller and J.
Schlenk for useful discussions. We would like to thank J. Zupan for
encouraging us to think about new applications of Machine Learning in HEP.
F.B. was supported by the ERC Starting Grant NewAve (638528). M.M. was
supported by the Swiss National Science Foundation, under Project Nos.
PP00P2176884. This research was partially supported by the Munich Institute
for Astro- and Particle Physics (MIAPP) which is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC-2094 – 390783311, the Aspen Center for Physics,
which is supported by National Science Foundation grant PHY-1607611, and the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306.

## References

  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Hirschi and Mattelaer (2015) V. Hirschi and O. Mattelaer, JHEP 10, 146 (2015), arXiv:1507.00020 [hep-ph] .
  * Boughezal  _et al._ (2017) R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu, F. Petriello, and C. Williams, Eur. Phys. J. C77, 7 (2017), arXiv:1605.08011 [hep-ph] .
  * Campbell and Neumann (2019) J. Campbell and T. Neumann, (2019), arXiv:1909.09117 [hep-ph] .
  * Nagy (2003) Z. Nagy, Phys. Rev. D68, 094002 (2003), arXiv:hep-ph/0307268 [hep-ph] .
  * Nagy (2002) Z. Nagy, Phys. Rev. Lett. 88, 122003 (2002), arXiv:hep-ph/0110315 [hep-ph] .
  * Gavin  _et al._ (2011) R. Gavin, Y. Li, F. Petriello, and S. Quackenbush, Comput. Phys. Commun. 182, 2388 (2011), arXiv:1011.3540 [hep-ph] .
  * Camarda  _et al._ (2019) S. Camarda _et al._ , (2019), arXiv:1910.07049 [hep-ph] .
  * Catani  _et al._ (2009) S. Catani, L. Cieri, G. Ferrera, D. de Florian, and M. Grazzini, Phys. Rev. Lett. 103, 082001 (2009), arXiv:0903.2120 [hep-ph] .
  * Catani  _et al._ (2015) S. Catani, D. de Florian, G. Ferrera, and M. Grazzini, JHEP 12, 047 (2015), arXiv:1507.06937 [hep-ph] .
  * de Florian  _et al._ (2011) D. de Florian, G. Ferrera, M. Grazzini, and D. Tommasini, JHEP 11, 064 (2011), arXiv:1109.2109 [hep-ph] .
  * Actis  _et al._ (2017) S. Actis, A. Denner, L. Hofer, J.-N. Lang, A. Scharf, and S. Uccirati, Comput. Phys. Commun. 214, 140 (2017), arXiv:1605.01090 [hep-ph] .
  * Hirschi  _et al._ (2011) V. Hirschi, R. Frederix, S. Frixione, M. V. Garzelli, F. Maltoni, and R. Pittau, JHEP 05, 044 (2011), arXiv:1103.0621 [hep-ph] .
  * Li  _et al._ (2016) Z. Li, J. Wang, Q.-S. Yan, and X. Zhao, Chin. Phys. C40, 033103 (2016), arXiv:1508.02512 [hep-ph] .
  * Naterop  _et al._ (2019) L. Naterop, A. Signer, and Y. Ulrich, (2019), arXiv:1909.01656 [hep-ph] .
  * Peraro (2019) T. Peraro, JHEP 07, 031 (2019), arXiv:1905.08019 [hep-ph] .
  * Bern  _et al._ (2014) Z. Bern, L. J. Dixon, F. Febres Cordero, S. Höche, H. Ita, D. A. Kosower, and D. Maitre, Comput. Phys. Commun. 185, 1443 (2014), arXiv:1310.7439 [hep-ph] .
  * Bertone  _et al._ (2017) V. Bertone, S. Carrazza, and N. P. Hartland, Comput. Phys. Commun. 212, 205 (2017), arXiv:1605.02070 [hep-ph] .
  * Maitre  _et al._ (2016) D. Maitre, G. Heinrich, and M. Johnson,  _Proceedings, 13th DESY Workshop on Elementary Particle Physics: Loops and Legs in Quantum Field Theory (LL2016): Leipzig, Germany, April 24-29, 2016_ , PoS LL2016, 016 (2016), arXiv:1607.06259 [hep-ph] .
  * Spira (2017) M. Spira, Prog. Part. Nucl. Phys. 95, 98 (2017), arXiv:1612.07651 [hep-ph] .
  * Bellagente  _et al._ (2019) M. Bellagente, A. Butter, G. Kasieczka, T. Plehn, and R. Winterhalder, (2019), arXiv:1912.00477 [hep-ph] .
  * Andreassen and Nachman (2019) A. Andreassen and B. Nachman, (2019), arXiv:1907.08209 [hep-ph] .
  * Nachman (2019) B. Nachman, (2019), arXiv:1909.03081 [hep-ph] .
  * Bendavid (2017) J. Bendavid, (2017), arXiv:1707.00028 [hep-ph] .
  * Otten  _et al._ (2019) S. Otten, S. Caron, W. de Swart, M. van Beekveld, L. Hendriks, C. van Leeuwen, D. Podareanu, R. Ruiz de Austri, and R. Verheyen, (2019), arXiv:1901.00875 [hep-ph] .
  * Klimek and Perelstein (2018) M. D. Klimek and M. Perelstein, (2018), arXiv:1810.11509 [hep-ph] .
  * Bothmann and Debbio (2019) E. Bothmann and L. Debbio, JHEP 01, 033 (2019), arXiv:1808.07802 [hep-ph] .
  * Ball  _et al._ (2017) R. D. Ball _et al._ (NNPDF), Eur. Phys. J. C77, 663 (2017), arXiv:1706.00428 [hep-ph] .
  * Ball  _et al._ (2015) R. D. Ball _et al._ (NNPDF), JHEP 04, 040 (2015), arXiv:1410.8849 [hep-ph] .
  * Santos  _et al._ (2017) R. Santos, M. Nguyen, J. Webster, S. Ryu, J. Adelman, S. Chekanov, and J. Zhou, JINST 12, P04014 (2017), arXiv:1610.03088 [hep-ex] .
  * Farina  _et al._ (2018) M. Farina, Y. Nakai, and D. Shih, (2018), arXiv:1808.08992 [hep-ph] .
  * D’Agnolo and Wulzer (2019) R. T. D’Agnolo and A. Wulzer, Phys. Rev. D99, 015014 (2019), arXiv:1806.02350 [hep-ph] .
  * Blance  _et al._ (2019) A. Blance, M. Spannowsky, and P. Waite, JHEP 10, 047 (2019), arXiv:1905.10384 [hep-ph] .
  * Cerri  _et al._ (2019) O. Cerri, T. Q. Nguyen, M. Pierini, M. Spiropulu, and J.-R. Vlimant, JHEP 05, 036 (2019), arXiv:1811.10276 [hep-ex] .
  * Brochet  _et al._ (2019) S. Brochet, C. Delaere, B. François, V. Lemaître, A. Mertens, A. Saggio, M. Vidal Marono, and S. Wertz, Eur. Phys. J. C79, 126 (2019), arXiv:1805.08555 [hep-ph] .
  * Banerjee  _et al._ (2019a) S. Banerjee, R. S. Gupta, J. Y. Reiness, S. Seth, and M. Spannowsky, (2019a), arXiv:1912.07628 [hep-ph] .
  * Komiske  _et al._ (2017) P. T. Komiske, E. M. Metodiev, and M. D. Schwartz, JHEP 01, 110 (2017), arXiv:1612.01551 [hep-ph] .
  * Coccaro  _et al._ (2019) A. Coccaro, M. Pierini, L. Silvestrini, and R. Torre, (2019), arXiv:1911.03305 [hep-ph] .
  * Andreassen  _et al._ (2019a) A. Andreassen, I. Feige, C. Frye, and M. D. Schwartz, Eur. Phys. J. C79, 102 (2019a), arXiv:1804.09720 [hep-ph] .
  * Andreassen  _et al._ (2019b) A. Andreassen, I. Feige, C. Frye, and M. D. Schwartz, Phys. Rev. Lett. 123, 182001 (2019b), arXiv:1906.10137 [hep-ph] .
  * de Oliveira  _et al._ (2016) L. de Oliveira, M. Kagan, L. Mackey, B. Nachman, and A. Schwartzman, JHEP 07, 069 (2016), arXiv:1511.05190 [hep-ph] .
  * Metodiev  _et al._ (2017) E. M. Metodiev, B. Nachman, and J. Thaler, JHEP 10, 174 (2017), arXiv:1708.02949 [hep-ph] .
  * Collins  _et al._ (2018) J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. Lett. 121, 241803 (2018), arXiv:1805.02664 [hep-ph] .
  * Collins  _et al._ (2019) J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. D99, 014038 (2019), arXiv:1902.02634 [hep-ph] .
  * Andreassen  _et al._ (2019c) A. Andreassen, P. T. Komiske, E. M. Metodiev, B. Nachman, and J. Thaler, (2019c), arXiv:1911.09107 [hep-ph] .
  * Paganini  _et al._ (2018a) M. Paganini, L. de Oliveira, and B. Nachman, Phys. Rev. Lett. 120, 042003 (2018a), arXiv:1705.02355 [hep-ex] .
  * Paganini  _et al._ (2018b) M. Paganini, L. de Oliveira, and B. Nachman, Phys. Rev. D97, 014021 (2018b), arXiv:1712.10321 [hep-ex] .
  * Chakraborty  _et al._ (2020) A. Chakraborty, S. H. Lim, and M. M. Nojiri, JHEP 19, 135 (2020), arXiv:1904.02092 [hep-ph] .
  * Lim and Nojiri (2018) S. H. Lim and M. M. Nojiri, JHEP 10, 181 (2018), arXiv:1807.03312 [hep-ph] .
  * Luisoni  _et al._ (2016) G. Luisoni, S. Poslavsky, and Y. Schroder,  _Proceedings, 17th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2016): Valparaiso, Chile, January 18-22, 2016_ , J. Phys. Conf. Ser. 762, 012077 (2016), arXiv:1604.03370 [hep-ph] .
  * Buckley (2019) A. Buckley, in  _19th International Workshop on Advanced Computing and Analysis Techniques in Physics Research: Empowering the revolution: Bringing Machine Learning to High Performance Computing (ACAT 2019) Saas-Fee, Switzerland, March 11-15, 2019_ (2019) arXiv:1908.00167 [hep-ph] .
  * Abdughani  _et al._ (2019) M. Abdughani, J. Ren, L. Wu, J. M. Yang, and J. Zhao, Commun. Theor. Phys. 71, 955 (2019), arXiv:1905.06047 [hep-ph] .
  * Albertsson  _et al._ (2018) K. Albertsson _et al._ ,  _Proceedings, 18th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2017): Seattle, WA, USA, August 21-25, 2017_ , J. Phys. Conf. Ser. 1085, 022008 (2018), arXiv:1807.02876 [physics.comp-ph] .
  * Albrecht  _et al._ (2019) J. Albrecht _et al._ (HEP Software Foundation), Comput. Softw. Big Sci. 3, 7 (2019), arXiv:1712.06982 [physics.comp-ph] .
  * Carrazza (2018) S. Carrazza,  _Proceedings, 18th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2017): Seattle, WA, USA, August 21-25, 2017_ , J. Phys. Conf. Ser. 1085, 022003 (2018), arXiv:1711.10840 [hep-ph] .
  * Larkoski  _et al._ (2017) A. J. Larkoski, I. Moult, and B. Nachman, (2017), arXiv:1709.04464 [hep-ph] .
  * Bourilkov (2019) D. Bourilkov, (2019), arXiv:1912.08245 [physics.data-an] .
  * Brehmer  _et al._ (2019a) J. Brehmer, F. Kling, I. Espejo, and K. Cranmer, (2019a), arXiv:1907.10621 [hep-ph] .
  * Brehmer  _et al._ (2019b) J. Brehmer, S. Dawson, S. Homiller, F. Kling, and T. Plehn, JHEP 11, 034 (2019b), arXiv:1908.06980 [hep-ph] .
  * Carleo  _et al._ (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019), arXiv:1903.10563 [physics.comp-ph] .
  * Note (1) The only similar works we found is Ref. Otten  _et al._ (2018) where DNN are used to predict NLO cross sections for the pMSSM-19 and Ref. Chawdhry  _et al._ (2019) ML regressors were used as interpolators to store the NNLO QCD amplitude for $pp\to 3\gamma$.
  * Cascioli  _et al._ (2014) F. Cascioli, T. Gehrmann, M. Grazzini, S. Kallweit, P. Maierhöfer, A. von Manteuffel, S. Pozzorini, D. Rathlev, L. Tancredi, and E. Weihs, Phys. Lett. B735, 311 (2014), arXiv:1405.2219 [hep-ph] .
  * Khachatryan  _et al._ (2017) V. Khachatryan _et al._ (CMS), JHEP 02, 135 (2017), arXiv:1610.09218 [hep-ex] .
  * Aaboud  _et al._ (2018) M. Aaboud _et al._ (ATLAS), Phys. Lett. B786, 59 (2018), arXiv:1808.08238 [hep-ex] .
  * Sirunyan  _et al._ (2018) A. M. Sirunyan _et al._ (CMS), Phys. Rev. Lett. 121, 121801 (2018), arXiv:1808.08242 [hep-ex] .
  * Banerjee  _et al._ (2018) S. Banerjee, C. Englert, R. S. Gupta, and M. Spannowsky, Phys. Rev. D98, 095012 (2018), arXiv:1807.01796 [hep-ph] .
  * Banerjee  _et al._ (2019b) S. Banerjee, R. S. Gupta, J. Y. Reiness, and M. Spannowsky, Phys. Rev. D100, 115004 (2019b), arXiv:1905.02728 [hep-ph] .
  * Arcadi  _et al._ (2019) G. Arcadi, A. Djouadi, and M. Raidal, (2019), arXiv:1903.03616 [hep-ph] .
  * Note (2) We determined this evaluation time using OpenLoops Buccioni  _et al._ (2019) to generate the exact $gg\to ZZ$ squared amplitudes for the training and prediction sets. The CPU core used for this timing is the same one used for all timings in this letter, see main text.
  * Abadi  _et al._ (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * Wright and Ziegler (2017) M. N. Wright and A. Ziegler, Journal of Statistical Software 77, 1 (2017).
  * Ho (1995) T. K. Ho, in  _Proceedings of the Third International Conference on Document Analysis and Recognition (Volume 1) - Volume 1_, ICDAR ’95 (IEEE Computer Society, Washington, DC, USA, 1995) pp. 278–.
  * Breiman (2001) L. Breiman, Machine Learning 45, 5 (2001).
  * Breiman (1997) L. Breiman,  _Arcing the edge_ , Tech. Rep. (Technical Report 486, Statistics Department, University of California at …, 1997).
  * Friedman (2001) J. H. Friedman, Ann. Statist. 29, 1189 (2001).
  * Chen and Guestrin (2016) T. Chen and C. Guestrin, Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining - KDD ’16 (2016), 10.1145/2939672.2939785.
  * Note (3) XGBoost rose to prominence by winning the Kaggle Higgs Callenge. Since then it has been consistently on the top of the ladder in a large number of the Kaggle competitions https://www.kaggle.com/sudalairajkumar/winning-solutions-of-kaggle-competitions ; https://www.kaggle.com/rajiao/winning-solutions-of-kaggle-competitions outperforming other ML architectures, like DNN’s or SVM’s.
  * Campbell and Ellis (2010) J. M. Campbell and R. K. Ellis,  _Proceedings, 10th DESY Workshop on Elementary Particle Theory: Loops and Legs in Quantum Field Theory: Woerlitz, Germany, April 25-30, 2010_ , Nucl. Phys. Proc. Suppl. 205-206, 10 (2010), arXiv:1007.3492 [hep-ph] .
  * Note (4) We also tested MadLoops Hirschi and Mattelaer (2015) and found it to be a factor of a few slower than OpenLoops Buccioni  _et al._ (2019) for the $gg\to ZZ$ process. For this reason, we used the faster evaluation time of OpenLoops as a benchmark.
  * Matsumoto and Nishimura (1998) M. Matsumoto and T. Nishimura, ACM Trans. Model. Comput. Simul. 8, 3 (1998).
  * Bothmann  _et al._ (2019) E. Bothmann _et al._ , SciPost Phys. 7, 034 (2019), arXiv:1905.09127 [hep-ph] .
  * Kilian  _et al._ (2018) W. Kilian, S. Brass, T. Ohl, J. Reuter, V. Rothe, P. Stienemeier, and M. Utsch, in  _International Workshop on Future Linear Collider (LCWS2017) Strasbourg, France, October 23-27, 2017_ (2018) arXiv:1801.08034 [hep-ph] .
  * Ke  _et al._ (2017) G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, in  _Advances in Neural Information Processing Systems 30_, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017) pp. 3146–3154.
  * Otten  _et al._ (2018) S. Otten, K. Rolbiecki, S. Caron, J.-S. Kim, R. Ruiz De Austri, and J. Tattersall, (2018), arXiv:1810.08312 [hep-ph] .
  * Chawdhry  _et al._ (2019) H. A. Chawdhry, M. L. Czakon, A. Mitov, and R. Poncelet, (2019), arXiv:1911.00479 [hep-ph] .
  * Buccioni  _et al._ (2019) F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang, and M. F. Zoller, Eur. Phys. J. C79, 866 (2019), arXiv:1907.13071 [hep-ph] .
  * (87) https://www.kaggle.com/sudalairajkumar/winning-solutions-of-kaggle-competitions, .
  * (88) https://www.kaggle.com/rajiao/winning-solutions-of-kaggle-competitions, .


it also holds, more generally, if $\ga_j\ge0$ for all $j$ and at
least two $\ga_j$ are non-zero.

Assume (<ref>). Then there exists a unique continuous function
$\gf(t)$ on $\oi$ such that $\gf(0)=0$, $\gf(1)=1$, and for

\begin{align}\label{q4}
\gf(t)
= \frac{\ga_{q-1-j}}{A}\gf(qt-j)
+ \frac{\sum_{q-j\le i<q}\ga_i}{A},
\qquad \text{if}\quad
\frac{j}{q}\le t\le \frac{j+1}{q}.
\end{align}

Moreover, if  $\ga_j\ge0$ for all $j$, then $\gf$ is strictly

This follows with only notational changes as in our third proof of
Lemma2 based on the recursive construction (<ref>); we
obtain by an analogue of (<ref>) a sequence of continuous
functions $\gf_k$ that, using (<ref>), converge uniformly to a
function $\gf(t)$ satisfying (<ref>).

We then extend $f(n)$ and $g(n)$ to functions of a real variable
$x\ge1$ by (<ref>) as before, and it is easily verified that
(<ref>) generalises to

\begin{align}\label{q5}
\qquad x\ge q.
\end{align}

We may now argue as in S:recurrence and prove extensions of
Theorem1 and its corollaries for the recursion (<ref>). We
now define

\begin{align}\label{q3}
\rho:=\log_q A
=\log\Bigpar{\sum_{0\le j<q}\ga_j}.
\end{align}

The simplest situation is when $g(n)=O(n^{\rho-\eps})$ for some
$\eps>0$; then

\begin{align}\label{q6}
f(x)=x^\rho P(\log_q x)-Q(x),
\qquad x\ge1,
\end{align}

where $P(t)$ is a continuous $1$-periodic function, and

\begin{align}\label{q7}
Q(x):=\sum_{m\ge1} q^{-\rho m} g(q^mx) = o(x^\rho).
\end{align}

We leave further details to the reader and content ourselves with the
discussion of two classes of examples.

§.§ Binomial coefficients not divisible by a prime $q$

Let $f(n)$ denote the number of binomial coefficients $\binom{m}{k}$,
$0\le k\le m<n$, that are not divisible by a given prime $q$. This
sequence has a long history, at least dating back to Fine's
[17] observation that almost all binomial coefficients are
even; see, e.g., [8] and the references therein. It
equals A006046 (see Example <ref>) when $q=2$. The case
$q=3$ corresponds to A006048, while $q=5$ gives A194458. We then
deduce the recurrence

\begin{equation}\label{qa1}
= \sum_{0\le j<q}(q-j)f\lpa{\tr{\tfrac{n+j}q}},
\end{equation}

with $f(0)=0$ and $f(1)=1$. This is (<ref>) with $\ga_j=q-j$ and
$g(n)=0$; furthermore, in this example (<ref>) holds for all
$n\ge0$. We have $A=\binom{q+1}2$ and $\rho := \log_q A$. Stein
[37] proved that

\begin{equation}\label{qa2}
\frac1A\le \frac{f(n)}{n^\rho}\le 1,
\end{equation}

and extended $f(n)$ to a continuous function $f(x)$; see
[21, 40, 42] for finer lower bounds. Our
general approach yields the same continuous extension $f(x)$ as in
[37]; we obtain (by (<ref>))

\begin{equation}\label{qa3}
\frac{f(n)}{n^{\rho}}=P(\log_qn)\qquad(n\ge1),
\end{equation}

where $P(t):=A^{-\frax t}f(q^{\frax t})$ is a continuous $1$-periodic
function. Moreover, since (<ref>) holds for all $n\ge1$, it is
easily verified (using (<ref>)) that $f(x)=A\gf(x/q)$ for
$x\in[1,q]$, and thus

\begin{equation}
P(t) = A^{1-\{t\}} \varphi(q^{\{t\}-1}).
\end{equation}

(Cf. Rb=1 for a similar simplification, for related reasons.)
Here $\varphi(t)$ satisfies $\varphi(0)=0$, $\varphi(1)=1$, and for
$0\le j<q$,

\begin{equation}
\varphi(t)
= \frac{j+1}{A}\varphi(\{qt\})
+ \frac{\binom{j+1}{2}}{A},
\qquad \text{if}\quad
\frac{j}{q}\le t\le \frac{j+1}{q}.
\end{equation}

By (<ref>) or (<ref>),
almost all binomial coefficients are divisible by any
given prime $q$ because

\begin{equation}
\rho = \log_q(q+1)-\log_q2+1 <2.
\end{equation}

$\varphi_3$     $\varphi_5$     $\varphi_7$     $\varphi_{11}$

$P_3$     $P_5$     $P_7$     $P_{11}$

The function $\varphi$ (SS:binomial-q) when $q=3, 5,
7, 11$ (upper half) and the periodic function $P$ for the same set of
values of $q$ (lower half).

Alternatively, it is known [8] that

\begin{equation}
f(n) = \frac12\sum_{0\le j\le s}
{\binom{p+1}2}^j b_j\prod_{j\le i\le s}(b_i+1),
\end{equation}

when $n = b_0 + b_1q+\cdots b_sq^s$ with $0\le b_j<q$ and
$s=\tr{\log_q n}$, from which we can obtain an alternative
representation for the periodic function $P$. Also the generating
function [37]

\begin{equation}
\sum_{n\ge0}f(n)z^n
= \frac{z}{1-z}\prod_{k\ge0}\sum_{0\le j<q}(j+1)z^{j\cdot q^k},
\end{equation}

is helpful in applying the Mellin transform approach; see
[20, 24, 25].

The more general problem of the number of multinomial coefficients
$\binom{m}{j_1,\dots,j_d}$ not divisible by a prime $q$, for $0\le
m<n$ and $j_1+\cdots+j_d=m$, $j_1,\dots,j_d\ge0$, where $d\ge1$ is
given (see [7, 8, 40]) can be similarly
dealt with. This number satisfies the recurrence

\begin{equation}
= \sum_{0\le j<q}\binom{q-j+d-2}{d-1}
\end{equation}

we then deduce, by (<ref>), the identity $f(n)= n^{\rho}
P(\log_qn)$, $n\ge1$, for some continuous periodic function $P$ with
$\rho = \log_q\binom{q+d-1}{d}$.

§.§ Generating polynomial of Gray codes

Gray codes of integers are strings of binary words in which
neighboring code words differ by one bit only; we already discussed
some properties of the binary reflected Gray codes in
Examples <ref> and <ref>. Here, we consider a simple
extended version of binary Gray codes to $q$-ary ones
(non-reflected); see [9, 36]. The construction is
as follows. If

\begin{equation}
n = \sum_{0\le j\le s}\kappa_j q^j,
\qquad(0\le \kappa_j<q),
\end{equation}

then the Gray code of $n$ is given by $(\kappa_s',\dots,\kappa_0')$,
where $\kappa_s' = \kappa_s$ and

\begin{equation}\label{E:kappa-p}
\kappa_j' := (\kappa_{j}-\kappa_{j+1})\bmod q
\qquad (0\le j<s).
\end{equation}

For simplicity, we consider the number of nonzero digits $\gamma(n):=
\sum_{0\le j\le s}\mathbf{1}_{\kappa_j'>0}$ in this Gray code
representation of $n$; other quantities such as the sum of digits
$\sum_{0\le j\le s}\kappa_j'$ can be considered similarly (a sketch
given below). Then by the recurrence

\begin{equation}
\gamma(qk+j)
= \gamma(k)+\begin{cases}
0, & \text{if } k \bmod q \equiv j;\\
1, & \text{otherwise},
\end{cases}
\end{equation}

for $0\le j<q$, we deduce that the generating polynomial $f(n) :=
\sum_{0\le k<n} \alpha^{\gamma(k)}$ of $\gamma(n)$ satisfies

\begin{equation}
= f\lpa{\ltr{\tfrac{n}{q}}}+
\alpha \sum_{1\le j<q} f\lpa{\ltr{\tfrac{n+j}{q}}}
\end{equation}

where $g$ can be expressed as

\begin{align}
= (1-\alpha)\alpha^{\gamma(qk+r)+1}
\qquad(k\ge0),
\end{align}

for $0\le r\le q-2$ and $r+1\le j\le q-1$, and $g(n)=0$ for all other
values of $n$. Alternatively, in terms of $q$-ary expansion, the
nonzero $g(n)$, $0\le n<q^2$, occurs when $n=(\kappa_1,\kappa_0)_q$
is of the form:

$(0,1)_q,(0,2)_q,(0,3)_q,(0,4)_q,\dots, (0,q-1)_q$

$(1,2)_q,(1,3)_q,(1,3)_q,\dots, (1,q-1)_q$

We then deduce from (<ref>) the exact and asymptotic expansion
(since $Q(n)=0$ by (<ref>))

\begin{equation}
f(n) = P(\log_qn)n^{\log_q(1+(q-1)\alpha)}
\qquad(n\ge1),
\end{equation}

whenever $\alpha>0$, for some continuous periodic function $P=P_\ga$;
note that then $g(n)=O\bigpar{(\ga\vee1)^{\log_q n}} =
O\bigpar{n^{0\vee\log_q\ga}}$ with $0\vee\log_q\ga<\rho

Periodic fluctuations of $f(n)n^{-\log_q(1+(q-1)\alpha)}$ in
the cases of nonzero digits of $q$-ary Gray codes for $q=3,4$ and
different $\alpha$ (SS:Gray).

Similarly, the generating polynomial of the sum-of-digits function in
such $q$-ary Gray codes $f(n) = \sum_{0\le k<n}\alpha^{\sigma(k)}$,
$\sigma(k) := \sum_{0\le j\le s}\kappa_j'$ when $k=\sum_{0\le j\le
s}\kappa_jq^j$, satisfies the recurrence

\begin{equation}
f(n) = \sum_{0\le j<q}\alpha^j f\lpa{\ltr{\tfrac{n+j}{q}}}
+ g(n),
\end{equation}

where $g$ can be expressed as

\begin{align}
\ga^{\gs(qk+r)}\frac{(1-\ga^j)(\ga^{q-r}-\ga^{q-j})}{1-\ga}
& (0\le j<r),\\
\ga^{\gs(qk+r)}\frac{(1-\ga^{j-r})(1-\ga^{q-j})}{1-\ga}
& (r\le j<q),
\end{cases}
\end{align}

for $k\ge0$, $0\le r\le q-1$, and $0\le j\le q-1$. This is derived by
the relation

\begin{align}
\sigma(qk+j) = \sigma(k) + (j-k) \bmod q,
\end{align}

which in turn follows from (<ref>). We then deduce that

\begin{equation}
f(n) = P(\log_qn)n^{\log_q(1+\alpha+\cdots + \alpha^{q-1})},
\end{equation}

for $n\ge1$.

§ MELLIN TRANSFORMS

Mellin transforms are another useful techniques in analysing
divide-and-conquer recurrences; see [18, 19, 20, 25, 28] and the
references therein for more information. Up to now most of the tools
we adopt to solve (<ref>) are direct and elementary in nature; it
is however possible to apply Mellin transforms for a more effective
characterisation of the underlying periodic oscillations, notably
calculations of the Fourier coefficients, as already observed before
in the literature (although an analytic approach often requires
stronger conditions).

Let again $\ga,\gb>0$, and $\rho:=\log_2(\ga+\gb)$, and assume that
$f(n)$ and $g(n)$ satisfy the recursion (<ref>). Extend again $f$
and $g$ to $[1,\infty)$ by (<ref>), with $g(1):=0$, and define
$f(x):=g(x):=0$ for $x\in[0,1)$. Denote the Mellin transform of
$f(x)$ by $f^*(s)$:

\begin{align}\label{mellin}
f^*(s):=\int_0^\infty f(x) x^{s-1}\dd x
= \int_1^\infty f(x) x^{s-1}\dd x,
\end{align}

for all complex $s$ such that the integral is absolutely convergent,
and similarly for $g^*(s)$. If $f(n)=O(n^c)$ for large $n$ and some
real $c$, then $f^*(s)$ exists at least in the half-plane $\Re s <
-c$, and is analytic there. If $f^*(s)$ or $g^*(s)$ extends
meromorphically to a larger domain, we use the same notation there.

Assume that $f^*(s)$ and $g^*(s)$ exist. Then (<ref>) implies that

\begin{align}\label{mel2}
&=\int_1^2 f(x)x^{s-1}\dd x
+ \int_2^\infty \Bigpar{2^\rho
f\Bigpar{\frac{x}{2}}+g(x)}x^{s-1}\dd x\notag\\
&=\int_1^2 \bigpar{f(x)-g(x)}x^{s-1}\dd x
+ 2^{\rho+s} f^*(s)+g^*(s).
\end{align}

Furthermore, (<ref>) or (<ref>) shows that if $1\le x<2$,
then $f(x)-g(x)=f(1)P_0(\log_2x)x^\rho$. Hence, recalling the
definition (<ref>),

\begin{align}\label{mel3}
\int_1^2 \bigpar{f(x)-g(x)}x^{s-1}\dd x &
=f(1)\int_1^2 P_0(\log_2 x) x^{\rho+s-1}\dd x
\notag\\&
=f(1)\int_1^2 \lrpar{1+(2^\rho-1)\gf\lrpar{x-1}} x^{s-1}\dd x.
\end{align}

Substituting this into (<ref>) yields

\begin{align}\label{mel4}
\lrpar{1-2^{\rho+s}}f^*(s)
=g^*(s) + f(1)\int_0^1\lrpar{1+(2^\rho-1)\gf(u)}
(1+u)^{s-1}\dd u.
\end{align}

Note that the right-hand side of (<ref>) equals, apart from a
factor $\frac1{\log 2}$, the right-hand side of (<ref>) at
$s=-\rho-\chi_k$. On the other hand, for such $s$, the left-hand side
factor $1-2^{\rho+s}$ equals zero. Hence, combining (<ref>) and
Theorem4<ref> yields the following (see
[19] for more information).

If $g(n)=O(n^{\rho-\eps})$ for some $\eps>0$, then $g^*(s)$ is
analytic in (at least) the half-plane $\Re s < -\rho+\eps$, and
$f^*(s)$ is meromorphic in the same half-plane with only simple
poles, which are at $-\rho-\chi_k=-\rho-\frac{2k\pi i}{\log 2}$ for
some $k\in\bbZ$. Furthermore, Theorem4<ref> applies, and

\begin{align}\label{mel5}
\hP(k)=-\Res\lrsqpar{f^*(s);s=-\rho-\chi_k}.
\end{align}

We do not claim that every $-\rho-\chi_k$ actually is a pole. In
fact, by (<ref>), $-\rho-\chi_k$ is a pole of $f^*(s)$ if
and only if $\hP(k)\neq0$.

For a better demonstration of the approach, we study the case when
$g(n)\sim n^\rho$, and thus Theorem1 does not directly apply
(although readily amenable) due to an extra logarithmic leading term
in the asymptotic approximation of $f(n)$. For a meromorphic function
$F(z)$, let $F(s)\FV$ denote finite value (or the constant
term) in the Laurent series expansion at $z=s$; thus $F(s)\FV=F(s)$
when the latter is finite.

Assume that $g(n)=n^\rho+g_0(n)$, where $g_0(n)=O(n^{\rho-\eps})$ for
some $\eps>0$. Then

\begin{equation}\label{tb1}
n^\rho\log_2 n +n^{\rho }P(\log_{2}n)-Q(n), \qquad n\ge1,
\end{equation}

where (as in Theorem1) $P(t)$ is a continuous $1$-periodic
function and $Q(n)=o\lrpar{n^\rho}$ as . The Fourier
coefficients of $P(t)$ are given by

\begin{equation}\label{c4B}
\widehat{P}(k)
=\frac{f(1)}{\log 2}\int_{0}^{1}
\frac{1+(\alpha +\beta -1)\varphi(u)}
{(1+u)^{\rho +\chi_{k}+1}}\dd u
\frac{1}{\log 2}g^*(-\rho-\chi_k), & k\neq0,
\\
\frac{1}{\log 2}g^*(-\rho)\FV+\frac12, & k=0,
\end{cases}
\end{equation}

where $g^*(s)$ is meromorphic in $\Re s<-\rho+(\eps\wedge1)$ with a
sole simple pole at $s=-\rho$. In particular, (<ref>) holds for

If, moreover, $g(n)=n^\rho$ for even $n\ge2$ (i.e., $g_0(2m)=0$),
then $Q(n)=0$ for $n\ge1$.

Assume, without loss of generality, $\eps<1$. The assumption and
(<ref>) then imply that we have, for $n\ge1$ and $t\in[0,1]$,

\begin{align}\label{tb2g}
\le |g(n+1)-g(n)|
\notag\\&
\end{align}

which together with a Taylor expansion of $(n+t)^{s-1}$ yields, for
$\Re s<-\rho$,

\begin{align}\label{mel12g}
=\sum_{n\ge1} \int_{0}^{1} g(n+t)(n+t)^{s-1}\dd t
=\sum_{n\ge1}\Bigpar{ g(n)n^{s-1}
+O\bigpar{n^{\rho+\Re s-1-\eps}}} \notag\\
&=\sum_{n\ge1} n^{\rho+s-1} +
\sum_{\ge 1} O\bigpar{n^{\rho+\Re s-1-\eps}}\notag\\
\sum_{n\ge1} O\bigpar{n^{\rho+\Re s-1-\eps}}.
\end{align}

Moreover, each term in the final sum is an entire function in $s$,
and the $O$ is uniform for $s$ in any compact set;
thus the sum converges to an analytic function in $H_\eps$.
Hence, $g^*(s)$ extends to a meromorphic function
in $H_\eps$, with a single simple pole at $s=-\rho$,
as asserted.

Let $f_1(n):=n^\rho\log_2 n$, and $f_2(n):=f(n)-f_1(n)$, and let
$g_j:=\gL_{\ga,\gb}[f_j]$, $j=1,2$. We have

\begin{align}\label{mel6}
\gL_{\ga,\gb}[f_1](2n)
\end{align}

and, with $\psi(x):=x^\rho\log_2(x)$, using the mean-value theorem.

\begin{align}\label{mel7}
\gL_{\ga,\gb}[f_1](2n+1)
&=(2n+1)^\rho+O\bigpar{n^{\rho-1}\log (n+1)}.
\end{align}

In other words, for all $n\ge2$,

\begin{align}\label{mel8}
g_1(n):=  \gL_{\ga,\gb}[f_1](n)
=n^\rho+O\bigpar{n^{\rho-1}\log n}
\end{align}

and thus

\begin{align}\label{mel9}
\end{align}

Consequently, $g_2^*(s)$ is analytic in the half-plane
$H_\eps:=\para{s:\Re s < -\rho+\eps}$, and, by C2,

\begin{equation}\label{mel10}
=x^{\rho }P_2(\log_{2}x)-Q_2(x), \qquad x\ge1,
\end{equation}

where $Q_2(x)=\sum_{m\ge1} 2^{-\rho m} g_2\bigpar{2^mx}$, and
$P_2(t)$ is a periodic continuous function. Since
$f(n)=f_1(n)+f_2(n)$, this shows (<ref>) with $P(t):=P_2(t)$ and

Furthermore, by LB1,

\begin{align}\label{mel11}
\hP(k)=
\hP_2(k)=-\Res\lrsqpar{f_2^*(s);s=-\rho-\chi_k}.
\end{align}

By (<ref>), $f_2^*(s)$ is meromorphic in $H_\eps$, with poles
only at $-\rho-\chi_k$. Moreover, similarly to (<ref>), it
follows from (<ref>) that
$|f_1(n+t)-f_1(n)|=O\lrpar{n^{\rho-\eps}}$ for $n\ge1$ and
$t\in[0,1]$, which in turn, similarly to (<ref>), yields

\begin{align}\label{mel12f}
&=\sum_{n\ge1} \int_{0}^{1} f_1(n+t)(n+t)^{s-1}\dd t
=\sum_{n\ge1} n^{\rho+s-1}\log_2n +
\sum_{n\ge1} O\bigpar{n^{\rho+\Re s-1-\eps}}
\notag\\
&=-\frac{1}{\log 2}\zeta'(1-\rho-s)+
\sum_{n\ge1} O\bigpar{n^{\rho+\Re s-1-\eps}}.
\end{align}

Again, each term in the final sum is an entire function in $s$, and
the $O$ is uniform for $s$ in any compact set; thus the sum converges
to an analytic function in $H_\eps$. Hence, $f_1^*(s)$ extends to a
meromorphic function in $H_\eps$, with a single (double) pole at
$s=-\rho$. Furthermore, the residue $\Res\lrsqpar{f_1^*(s);s=-\rho} =
\frac{-1}{\log 2} \Res\lrsqpar{\zeta'(1-\rho-s);s=-\rho}=0$, since
$\zeta'$ is a derivative of a meromorphic function. Consequently,
$f^*(s)=f_1^*(s)+f_2^*(s)$ is meromorphic in $H_\eps$, with poles
only at $-\rho-\chi_k$, and

\begin{align}\label{mel14}
\Res\lrsqpar{f^*(s);s=-\rho-\chi_k}
\qquad k\in\bbZ.
\end{align}

Thus, by (<ref>) and (<ref>),

\begin{align}\label{mel16a}
\hP(k)&
\end{align}

for $k\neq0$,
(<ref>) yields,
with $s=-\rho-\chi_k$,

\begin{align}\label{mel16}
\hP(k)\log 2
+ f(1)\int_0^1 \lrpar{1+(2^\rho-1)\gf(u)}
(1+u)^{-\rho-\chi_k-1}\dd u,
\end{align}

which is (<ref>) in this case. For $k=0$, by the expansions
$\zeta(z)=(z-1)\qw+O(1)$ and $1-2^z=-(\log
2)z-\frac{(\log2)^2}{2}z^2+O(z^3)$ for small $|z|$, and the relations
(<ref>) and (<ref>), we have

\begin{align}
\notag\\
&=-z\qw-\frac{\log 2}{2}+\hP(0)\log 2+O(z).\label{mel20}
\end{align}

\begin{align}\label{mel21}
\hP(0) = \frac12+\frac{1}{\log 2}
\bigpar{(1-2^{\rho+s})f^*(s)\big|_{s=-\rho}}\FV,
\end{align}

and (<ref>) for $k=0$ follows from (<ref>).

Finally, (<ref>) shows that $g_1(n)=n^\rho$ for even $n$. Hence,
if $g(n)=n^\rho$ for even $n$, then $g_2(2m)=0$ for $m\ge1$, and then
$Q(n)=Q_2(n)=0$ for $n\ge1$ by Eodd, which proves the final
claim of the theorem.

Let $\ga=\gb=2$ (so $\rho=2$), and let $g(n)=n^2$, $n\ge2$. Then
TB1 applies, with $g_0(n)=0$. Hence, $Q(n)=0$, and (<ref>)

\begin{equation}\label{ve0}
n^\rho\log_2 n +n^{\rho }P(\log_{2}n), \qquad n\ge1.
\end{equation}

Furthermore, (<ref>) shows that, at least for $\Re s<-2$,

\begin{align}
g^*(s) = \frac{D(-s-1)}{s(s+1)},
\end{align}

where (<ref>) yields, recalling $g(0)=g(1)=0$,

\begin{align}
\qquad \Re s >1.
\end{align}

\begin{align}\label{ve1}
\end{align}

This shows that $g^*(s)$ is meromorphic in $\bbC$, with poles at
$0,-1$ and $-2$. The finite value $g^*(-2)\FV=\gamma-\frac{3}4$,
where $\gamma$ is Euler's constant. Hence, if for simplicity $f(1)=0$,
the Fourier series of $P(t)$ is, by (<ref>) and (<ref>),

\begin{align}\label{ve2}
P(t) = \frac{4\gamma-3}{4\log 2}+\frac12
+\frac{1}{2\log 2}\sum_{k\neq 0}\frac{4\zeta(1+\chi_k)+3}
{(1+\chi_k)(2+\chi_k)}\,e^{2k\pi i t}.
\end{align}

This sequence $f(n)$ with $g(n)=n^2$ and $f(1)=0$ is not in OEIS, but
the sequence $f_{\A022560}(n)=\A022560(n-1)$ discussed in
E22-ns satisfies the recurrence $\Lambda_{2,2}[f]=
\tr{\frac14n^2}$ with $f_{\A022560}(1)=0$. Hence,
$4f_{\A022560}(n)=f(n)+n^2-S_{2,2}(n)$ with $f(n)$ as above. The
formulas (<ref>)–(<ref>) now follow from (<ref>)
and (<ref>) together with (<ref>).

§ A SERIES REPRESENTATION FOR $\HAT{P}(K)$

We prove (<ref>) in Remark <ref>. For notational
simplicity we consider the case $k=0$; the general case is the same,
with $\rho$ replaced by $\rho+\chi_k$ below. (Note that
$2^{\rho+\chi_k}=2^\rho=\ga+\gb$ for all $k\in\bbZ$.) Consider the

\begin{align}
J := \int_0^1 \frac{1+(\alpha+\beta-1)\varphi(t)}
{(1+t)^{\rho+1}}\,\dd t
= \frac{1-2^{-\rho}}{\rho} + (\alpha+\beta-1)J(1)
\end{align}

\begin{align}
J(m) := \int_0^1 \frac{\varphi(t)}
{(m+t)^{\rho+1}}\,\dd t.
\end{align}

We now express $J(1)$ in a series form as follows. First, by applying
the recursive definition (<ref>) of $\varphi$:

\begin{align}
\varphi(t)
= \begin{cases}
\frac{\beta}{\alpha+\beta}\,\varphi(2t),
&\text{if }0\le t\le\frac12;\\
\frac{\alpha}{\alpha+\beta}\,\varphi(2t-1)
&\text{if }\frac12\le t\le1,
\end{cases}
\end{align}

we obtain, for any $m\ge1$,

\begin{align}
J(m) &:= \int_0^1 \frac{\varphi(t)}
{(m+t)^{\rho+1}}\,\dd t\notag\\
&= \beta\int_0^1
\frac{\varphi(t)}{(2m+t)^{\rho+1}}\,\dd t
+ \alpha\int_0^1
\frac{\varphi(t)}{(2m+1+t)^{\rho+1}}\,\dd t
+ \frac{\gb}{\ga+\gb}
\int_{\frac12}^1\frac1{(m+t)^{\rho+1}}\,\dd t\notag\\
&=: \beta J(2m)+\alpha J(2m+1) +K(m),
\end{align}

say. Iterating this gives, for any $N\ge0$,

\begin{align}\label{cc4}
J(1) = \sum_{0\le m<N}
\sum_{0\le j<2^m}\alpha^{\nu(j)}\beta^{m-\nu(j)}K(2^m+j)
+\sum_{0\le j<2^N}\alpha^{\nu(j)}\beta^{N-\nu(j)}J(2^N+j).
\end{align}

where $\nu(j)$ denotes the number of 1's in $j$'s binary expansion.
Since $|J(m)|,|K(m)| = O(m^{-\rho-1})$, the last sum in (<ref>)
is $O\bigpar{(\ga+\gb)^N2^{-(\rho+1)N}}=O\bigpar{2^{-N}}$ and
similarly the inner sum in the double sum is $O(2^{-m})$; hence we
can let $N\to\infty$, which yields

\begin{align}
J(1) = \sum_{m\ge0}\sum_{0\le j<2^m}
\alpha^{\nu(j)}\beta^{m-\nu(j)}K(2^m+j).
\end{align}

\begin{align}
\int_{\frac12}^1\frac1{(m+t)^{\rho+1}}\,\dd t
= \frac{\beta}{\rho(\alpha+\beta)}
\Lpa{\frac1{(m+\tfrac12)^{\rho}}
\end{align}

\begin{align}
&\int_0^1 \frac{1+(\alpha+\beta-1)\varphi(t)}
{(1+t)^{\rho+1}}\,\dd t \notag\\
&\quad = \frac{1-2^{-\rho}}{\rho}
+ \frac{\beta(\alpha+\beta-1)}{\rho(\alpha+\beta)}
\sum_{m\ge0}\sum_{0\le j<2^m}
\alpha^{\nu(j)}\beta^{m-\nu(j)}
\Lpa{\frac1{(2^m+j+\tfrac12)^{\rho}}
&\quad= \frac{\alpha+\beta-1}{\rho(\alpha+\beta)}
\llpa{1+\beta \sum_{m\ge0}\sum_{0\le j<2^m}
\alpha^{\nu(j)}\beta^{m-\nu(j)}
\Lpa{\frac1{(2^m+j+\tfrac12)^{\rho}}
\end{align}

The double sum can be converted into a single one as follows.

\begin{align}
&\sum_{m\ge0}\sum_{0\le j<2^m}
\alpha^{\nu(j)}\beta^{m+1-\nu(j)}
\Lpa{\frac1{(2^m+j+\tfrac12)^{\rho}}
&\qquad= 2^\rho \sum_{m\ge0}\sum_{0\le j<2^m}
\alpha^{\nu(j)}\beta^{m+1-\nu(j)}
\Lpa{\frac1{(2^{m+1}+2j+1)^{\rho}}
&\qquad= 2^\rho \sum_{m\ge1}\sum_{0\le j<2^{m-1}}
\alpha^{\nu(j)}\beta^{m-\nu(j)}
\Lpa{\frac1{(2^{m}+2j+1)^{\rho}}
-\frac1{(2^{m}+2j+2)^{\rho}}} \notag\\
&\qquad= 2^\rho\sum_{k\ge2}\frac{(-1)^k}{(k+1)^{\rho}}
\,\alpha^{\nu(\tr{2^{L_k-1}\{k/2^{L_k}\}})}
\beta^{{L_k}-\nu(\tr{2^{L_k-1}\{k/2^{L_k}\}})}.
\end{align}

§ RECURRENCES WITH MINIMISATION OR MAXIMISATION

Consider the class of sequences satisfying recurrences of the form

\begin{align}\label{ol1}
=\min_{1\le k\le
\left\lfloor \frac{n}{2}\right\rfloor}
\left\{ \alpha u(k)+\beta u(n-k)\right\}
\qquad (n\ge 2),
\end{align}

with $u(1)=1$.

This was studied in [5], where it was shown that if
$\alpha$ and $\beta$ are positive integers with $\beta \ge
\alpha$, then the minimum in (<ref>) is reached at
$k=\left\lfloor \frac{n}{2}\right\rfloor $ and the solution is given

\begin{equation}\label{e1}
=1+(\alpha +\beta -1)\sum_{1\le j<n}w(j),
\end{equation}

\begin{equation}\label{e10}
w(j):=\alpha ^{L_{j}-\nu _{0}(j)}\beta^{\nu _{0}(j)},
\end{equation}

and $\nu_{0}(j)$ denotes the number of zeros in the binary expansion
of $j$. We will extend this result and the arguments in
[5], and prove the following.

Let $\ga,\gb>0$ be real numbers such that either

$\beta \ge \alpha $ and $\beta \ge 1$, or

$\alpha \ge\beta $ and $\alpha +\beta \le1$.

Then the minimum in (<ref>) is reached at $k=\left\lfloor
\frac{n}{2}\right\rfloor $. Hence, (<ref>) reduces to
$\Lambda_{\ga,\gb}[u]=0$, and thus (<ref>) is solved by

We note first that $S_\gab(n)$ is given by the formula in (<ref>)
for any $\ga,\gb$. This follows by (<ref>), (<ref>), and
(<ref>) in S:recurrence, or by the proof in

It thus remains to show that if $u(n)$ is defined by (<ref>),

\begin{align}\label{ol2}
\alpha u(k)+\beta u(n-k)
\ge \alpha u\lpa{\left\lfloor \tfrac{n}{2}
\right\rfloor}
+\beta u\lpa{\left\lceil \tfrac{n}{2}\right\rceil}
\end{align}

for $1\le k<\left\lfloor \frac{n}{2}\right\rfloor $. By (<ref>),
the difference between the two sides of (<ref>) is

\begin{align}\label{ol3}
(\alpha +\beta -1)\left(\beta
\sum_{\left\lceil \frac{n}{2}\right\rceil\le j<n-k}
w(j)-\alpha \sum_{k\le j<\left\lfloor\frac{n}{2}\right\rfloor}
\end{align}

To prove that this is non-negative, we will show that

\begin{equation}\label{e2}
\beta \sum_{\left\lceil \frac{n}{2}\right\rceil\le j<n-k}
w(j)\ge \alpha \sum_{k\le j<\left\lfloor
\frac{n}{2}\right\rfloor}w(j)
\end{equation}

if <ref> holds, and that (<ref>) holds with the inequality
reversed if <ref> holds. The key is the following claim:

If $\beta \ge \alpha $ and $\beta \ge 1$, then

\begin{align}\label{paaa}
\beta w(n+2^{j})
\ge \alpha w(n),\quad \text{for all }
n\ge 1\text{ and }j\ge 0.
\end{align}

If $\alpha \ge \beta $ and $\beta \le 1$, then

\begin{align}\label{paab}
\beta w(n+2^{j})\le \alpha w(n),\quad
\text{for all }n\ge 1\text{ and }j\ge 0.
\end{align}

Proof of the claim. From (<ref>), we have

\begin{align}\label{e5}
\frac{\beta w(n+2^{j})}{\alpha w(n)}&
=\left( \frac{\beta }{\alpha }\right)^{
1+\nu _{0}(n+2^{j})-\nu _{0}(n)}
\alpha^{L_{n+2^{j}}-L_{n}} \notag\\
&=\left( \frac{\beta }{\alpha }\right)^{
1+\nu _{0}(n+2^{j})-\nu _{0}(n)-(L_{n+2^{j}}-L_{n})}
\beta^{L_{n+2^{j}}-L_{n}}.
\end{align}

It is easily seen that

\begin{equation}\label{e6}
1+\nu _{0}(n+2^{j})-\nu _{0}(n)
\ge L_{n+2^{j}}-L_{n}\ge 0.
\end{equation}

Both parts of the claim thus follow from (<ref>).

To show (<ref>), or its converse in case <ref>, which will
then complete the proof of the proposition, we combine the claim
above with a pairing between the sets, with

\begin{equation}
\left\{\left\lfloor \frac{n}{2}\right\rfloor -m,
\cdots ,\left\lfloor \frac{n}{2}\right\rfloor -1\right\}
\quad \text{and}\quad \left\{ \left\lceil
\frac{n}{2}\right\rceil,
\cdots ,\left\lceil \frac{n}{2}\right\rceil +m-1\right\},
\end{equation}

such that the difference between the elements of each pair is a power
of $2$. In other words, the proof is completed by applying the
following lemma (with a translation).

Let $n\ge1$ and

\begin{align}
A &:=\{1,2,\cdots ,n\}, \\
C_{1} &:=\{n+1,n+2,\cdots ,2n\}, \\
C_{2} &:=\{n+2,n+3,\cdots ,2n+1\}.
\end{align}

There exist one-to-one mappings

\begin{equation}
h_{1}:A\rightarrow C_{1}\quad
\text{and}\quad h_{2}:A\rightarrow C_{2},
\end{equation}

such that for each $k\in A$ there exist $j_{1},j_{2}$ with

\begin{align}
\text{for some } j_1\ge0, \label{Ah1}\\
&=k+2^{j_{1}},\quad \text{for some } j_2\ge0.\label{Ah2}
\end{align}

We prove the existence of $h_1$ by induction. Write $n=2^{\ell}+j$
for some $\ell\ge 0$ and $0\le j< 2^{\ell}$. (Thus,
$\ell=L_n$.) We want to show that there exist a one-to-one mapping

\begin{equation}
h_{1}:\{1,\cdots ,2^{\ell}+j\}
\rightarrow \{2^{\ell}+j+1,\cdots ,2^{\ell+1}+2j\},
\end{equation}

such that for each $1\le k \le 2^{\ell}+j$, (<ref>)

If $j=0$, we simply define $h_1(k):=k+2^\ell$.

If $j\ge1$, we first define $h_1(k)$ for $k\le 2j$ by
$h_1(k):=k+2^{\ell+1}$. This gives a mapping from $\{1,2,\cdots
,2j\}$ to $\{2^{\ell+1}+1,\cdots ,2^{\ell+1}+2j\}$. We remove these
two blocks, and it remains to define a one-to-one mapping satisfying
(<ref>) between $\{2j+1,\cdots ,2^{\ell}+j\}$ and
$\{2^{\ell}+j+1,\cdots ,2^{\ell+1}\}$. By subtracting $2j$ from each
term, it is equivalent to showing that there exist such a mapping
from $\{1,\cdots ,m\} $ to $\{m+1,\cdots ,2m\}$, where
$m=2^{\ell}-j$; this is true by the induction hypothesis.

The proof of the existence of $h_2$ is similar (but with
$n+1=2^\ell+j$); we omit the details.

If we replace $\min$ by $\max$ in (<ref>), we obtain a similar
result; see also [32] for $\ga>\gb=1$.

Let $\ga,\gb>0$ be real numbers such that $\alpha \ge \beta$,
$\beta \le 1$, and $\alpha+\beta \ge 1$. Then the maximum
in the recursion

\begin{align}\label{olmax}
u(n)=\max_{1\le k\le
\left\lfloor \frac{n}{2}\right\rfloor}
\left\{ \alpha u(k)+\beta u(n-k)\right\}
\qquad (n\ge 2),
\end{align}

with $u(1)=1$, is attained at $k=\left\lfloor
\frac{n}{2}\right\rfloor $. Hence, (<ref>) reduces to
$\Lambda_{\ga,\gb}[u]=0$, and thus (<ref>) is solved by

The proof above shows that under these conditions, (<ref>)
holds, and hence (<ref>) holds in the opposite direction. Thus the
difference in (<ref>) is $\le0$. (We do not obtain a second case
with $\ga+\gb<1$; we then would need (<ref>), but the condition
$\gb\ge1$ in <ref> above is incompatible with $\ga+\gb<1$.)

The results above can be extended to the recurrences of the form

\begin{align}
u(n)=\min_{1\le k\le
\left\lfloor \frac{n}{2}\right\rfloor
}\left\{ \alpha u(k)+\beta u(n-k)\right\}
+c\qquad (n\ge 2),
\end{align}

with the same conditions as above on $\ga$ and $\gb$; see
[22, 28] for more general versions. The solution
for the recurrence

\begin{align}
=\ga u\left( \left\lfloor \frac{n}{2}\right\rfloor \right)
+\gb u\left( \left\lceil \frac{n}{2}\right\rceil \right) +c
\end{align}

is, with $w(j)$ as in (<ref>)–(<ref>) above (see also

\begin{equation}
=u(1)+\bigpar{(\alpha +\beta -1)u(1)+c}
\sum_{1\le j<n}w(j).
\end{equation}

§ NOWHERE DIFFERENTIABILITY OF $P_{\TEXT{A006581}}(T)$

We prove in this appendix the fractal nature of the periodic function
$P(t)$ arising from A006581 (discussed in E22-ns with Fourier
expansion given in (<ref>)), namely,
$\Lambda_{2,2}[f]=g$ where $g(n) := \{\frac{n}2\}(n-1)$:

\begin{align}
f(n) = n^2P(\log_2n),
\end{align}

where, with $\bar{g}(x) := g(x)/x^2$,

\begin{align}
P(t) = \sum_{m\in\mathbb{Z}} 4^{-m-t}g(2^{m+t})
=\sum_{m\in\mathbb{Z}} \bar{g}(2^{m+t}).
\end{align}

Here $g(x)$ is extended from $g(n)$ as in (<ref>) with $\varphi(t)

\begin{align}
\bar{g}(x) = \frac1{2x^2}\times\begin{cases}
\{x\}\tr{x},
&\text{if }\tr{x} \text{ is even};\\
&\text{if }\tr{x} \text{ is odd}.
\end{cases}
\end{align}

Since $\bar{g}(2^{m+t})=0$ for $t\in[0,1]$ and $m\le 0$, we have

\begin{align}
P(t) = \sum_{m\ge1} \bar{g}(2^{m+t})
\qquad (0\le t\le 1).
\end{align}

The method of proof used here to prove the nowhere differentiability
of $P$ is standard and similar to that for the Takagi function given
in the survey paper [1].

Let $t\in[0,1)$, and define

\begin{align}
\tau_n := \log_2\frac{\tr{2^{n+t}}}{2^n}
\eqtext{and}
\tau_n' := \log_2\frac{\tr{2^{n+t}}+1}{2^n}.
\end{align}

To prove that $P$ does not have a finite derivative at $t$, it
suffices to show that the sequence

\begin{align}\label{E:P-ratio}
\frac{P(\tau_n)-P(\tau_n')}
= \sum_{1\le m\le n}\frac{\bar{g}(2^{m+\tau_n})-
\bar{g}(2^{m+\tau_n'})}{\tau_n-\tau_n'}
\end{align}

does not converge to a finite limit. Here we used the relation
$\bar{g}(2^{m+\tau_n})=\bar{g}(2^{m+\tau_n'})=0$ for $m>n$. Now
for $\theta\in[0,1)$

\begin{align}
\left\{
\begin{aligned}
\bar{g}(2k+\theta)
&= \frac{k\theta}{(2k+\theta)^2}\\
\bar{g}(2k+1+\theta)
&= \frac{k(1-\theta)}{(2k+1+\theta)^2}
\end{aligned}
\right.
\end{align}

so that, taking the right derivative at integer points here and below,

\begin{align}\label{ccg'}
\left\{
\begin{aligned}
\bar{g}'(2k+\theta)
&= \frac{k(2k-\theta)}{2(2k+\theta)^3}\\
\bar{g}'(2k+1+\theta)
&= -\frac{k(2k+3-\theta)}{(2k+1+\theta)^3}
.	\end{aligned}
\right.
\end{align}

If $1\le m\le n$, then

\begin{align}\label{cc5}
\floor{2^{m+\tau_n}} \le 2^{m+\tau_n}
\le 2^{m+t}
< 2^{m+\tau'_n}
\le \floor{2^{m+\tau_n}}+1.
\end{align}

It follows that $h(x):=\bar g(2^x)$ is infinitely differentiable on
$[m+\tau_n,m+\tau_n']$, and it is easily seen from (<ref>) that
$h''(x)=O(1)$ (uniformly in $m$ and $n$). We have, for some

\begin{align}
\frac{\bar{g}(2^{m+\tau_n})-
\bar{g}(2^{m+\tau_n'})}{\tau_n-\tau_n'}
= h'(m+\tau''_n)\notag\\
\end{align}

Hence, (<ref>) implies that

\begin{align}
\frac{P(\tau_n)-P(\tau_n')}{\tau_n-\tau_n'}
= \sum_{1\le m\le n}h'(m+t)+O(n2^{-n}),
\end{align}

and thus, if $P$ is differentiable at $t$, then the sum

\begin{align}\label{cc6}
\sum_{ m\ge 1}h'(m+t)
\end{align}

converges (and equals $P'(t)$).

On the other hand, it follows easily from (<ref>) that if
$x\ge2$, then $x\bar g'(x)\ge \frac19$ for even $\tr{x}$ and $x\bar
g'(x)\le -\frac14$ for odd $\tr{x}$. Hence, $|h'(m+t)|=|2^{m+t}\bar
g'(2^{m+t})\log 2|\ge \frac{\log2}9$ for all $m\ge1$. Consequently,
the sum (<ref>) diverges for any $t$, and thus $P$ is nowhere

Moreover, we note that if $2^t$ is a dyadic rational, then for all
large $m$, $2^{m+t}$ is an even integer, and thus
$h'(m+t)\ge\frac{\log 2}{9}$. Hence, in this case the sum (<ref>)
diverges to $+\infty$, and thus so does
$\bigpar{P(\tau_n)-P(\tau'_n)}/(\tau_n-\tau'_n)$ in
(<ref>); consequently, $P$ is not Lipschitz.

We do not know whether $P$ is Hölder continuous, and leave that
as an open problem. Note that LProp1 does not apply since
(<ref>) does not hold.

[1]
P. C. Allaart and K. Kawamura.
The Takagi function: a survey.
Real Anal. Exchange, 37(1):1–54, 2011/12.

[2]
J. L. Bentley, D. Haken, and J. B. Saxe.
A general method for solving divide-and-conquer recurrences.
ACM SIGACT News, 12(3):36–44, 1980.

[3]
R. C. Bose and R. J. Nelson.
A sorting problem.
J. ACM, 9(2):282–296, Apr. 1962.

[4]
N. Campbell.
Recurrences, 2020.
Unpublished treatise available at

[5]
K.-N. Chang and S.-C. Tsai.
Exact solution of a minimal recurrence.
Inform. Process. Lett., 75(1-2):61–64, 2000.

[6]
L. H. Y. Chen, H.-K. Hwang, and V. Zacharovas.
Distribution of the sum-of-digits function of random integers: a
Probab. Surv., 11:177–236, 2014.

[7]
Y.-G. Chen and C. Ji.
The number of multinomial coefficients not divided by a prime.
Acta Sci. Math. (Szeged), 64(1-2):37–48, 1998.

[8]
Y. G. Chen and C. G. Ji.
On a function related to multinomial coefficients. I.
Acta Math. Sin. (Engl. Ser.), 18(4):647–660, 2002.

[9]
M. Cohn.
Affine $m$-ary Gray codes.
Information and Control, 6(1):70–78, 1963.

[10]
T. H. Cormen, C. E. Leiserson, and R. L. Rivest.
Introduction to Algorithms.
MIT press, 1990.

[11]
T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein.
Introduction to Algorithms.
MIT press, 2022.

[12]
H. Delange.
Sur la fonction sommatoire de la fonction “somme des chiffres”.
Enseign. Math. (2), 21(1):31–47, 1975.

[13]
M. Drmota and W. Szpankowski.
A master theorem for discrete divide and conquer recurrences.
Journal of the ACM, 60(3):1–49, 2013.

[14]
K. Ellul, B. Krawetz, J. Shallit, and M.-W. Wang.
Regular expressions: new results and open problems.
J. Autom. Lang. Comb., 10(4):407–437, 2005.

[15]
J. Erickson.
Sowing games.
In Games of No Chance, pages 287–297. Cambridge University
Press, 1996.

[16]
S. R. Finch.
Mathematical Constants.
Cambridge University Press, 2003.

[17]
N. J. Fine.
Binomial coefficients modulo a prime.
Amer. Math. Monthly, 54:589–592, 1947.

[18]
P. Flajolet and M. Golin.
Mellin transforms and asymptotics. The mergesort recurrence.
Acta Inform., 31(7):673–696, 1994.

[19]
P. Flajolet, X. Gourdon, and P. Dumas.
Mellin transforms and asymptotics: Harmonic sums.
Theoret. Comput. Sci., 144(1-2):3–58, 1995.

[20]
P. Flajolet, P. Grabner, P. Kirschenhofer, H. Prodinger, and R. F. Tichy.
Mellin transforms and asymptotics: digital sums.
Theoret. Comput. Sci., 123(2):291–314, 1994.

[21]
Z. M. Franco.
Distribution of binomial coefficients modulo three.
Fibonacci Quart., 36(3):272–275, 1998.

[22]
M. L. Fredman and D. E. Knuth.
Recurrence relations based on minimization.
J. Math. Anal. Appl., 48:534–559, 1974.

[23]
O. Garet.
A simple master theorem for discrete divide and conquer recurrences.
North-Western European Journal of Mathematics, 8:91–101, 2022.

[24]
P. J. Grabner.
Completely $q$-multiplicative functions: the Mellin transform
Acta Arith., 65(1):85–96, 1993.

[25]
P. J. Grabner and H.-K. Hwang.
Digital sums and divide-and-conquer recurrences: Fourier expansions
and absolute convergence.
Constr. Approx., 21(2):149–179, 2005.

[26]
C. Heuberger, D. Krenn, and G. F. Lipnik.
Asymptotic analysis of $q$-recursive sequences.
Algorithmica, 84(9):2480–2532, 2022.

[27]
H.-K. Hwang, S. Janson, and T.-H. Tsai.
Exact and asymptotic solutions of a divide-and-conquer recurrence
dividing at half: theory and applications.
ACM Trans. Algorithms, 13(4):Art. 47, 43 pages, 2017.

[28]
H.-K. Hwang and T.-H. Tsai.
An asymptotic theory for recurrence relations based on minimization
and maximization.
Theoret. Comput. Sci., 290(3):1475–1501, 2003.

[29]
A. Ivić.
The Riemann Zeta-Function.
John Wiley & Sons, Inc., New York, 1985.

[30]
Y. Katznelson.
An Introduction to Harmonic Analysis.
John Wiley & Sons, Inc., New York-London-Sydney, 1968.

[31]
D. E. Knuth.
The Art of Computer Programming, Volume 2: Seminumerical
Addison-Wesley, Boston, third edition, 1997.

[32]
D. E. Knuth and O. P. Lossers.
Solving a recurrence by binary expansion. Solutions 11320 [2007,
Amer. Math. Monthly, 116(7):649, 2009.

[33]
W. Kuszmaul and C. E. Leiserson.
Floors and ceilings in divide-and-conquer recurrences.
In Symposium on Simplicity in Algorithms (SOSA), pages
133–141. SIAM, 2021.

[34]
OEIS Foundation Inc.
The On-Line Encyclopedia of Integer Sequences
(OEIS), 2022.
Published electronically at <http://oeis.org>.

[35]
S. Roura.
Improved master theorems for divide-and-conquer recurrences.
Journal of the ACM, 48(2):170–205, 2001.

[36]
B. D. Sharma and R. K. Khanna.
On $m$-ary Gray codes.
Inform. Sci., 15(1):31–43, 1978.

[37]
A. H. Stein.
Binomial coefficients not divisible by a prime.
In Number theory (New York, 1985/1988), volume 1383 of Lecture Notes in Math., pages 170–177. Springer, Berlin, 1989.

[38]
K. B. Stolarsky.
Power and exponential sums of digital sums related to binomial
coefficient parity.
SIAM Journal on Applied Mathematics, 32(4):717–730, 1977.

[39]
E. C. Titchmarsh.
The Theory of the Riemann Zeta-Function.
The Clarendon Press, Oxford, 1951.

[40]
N. A. Volodin.
Multinomial coefficients modulo a prime.
Proc. Amer. Math. Soc., 127(2):349–353, 1999.

[41]
H. von Koch.
Sur une courbe continue sans tangente, obtenue par une construction
géométrique élémentaire.
Arkiv för Matematik, Astronomi och Fysik, 1:681–704, 1904.

[42]
B. Wilson.
Asymptotic behavior of Pascal's triangle modulo a prime.
Acta Arith., 83(2):105–116, 1998.

[43]
A. Zygmund.
Trigonometric Series. Vol. I, II.
Cambridge University Press, Cambridge, second edition, 1959.