Neelabh Sinha

Michal Balazia

Franc ¸ois Bremond

FLAME: Facial Landmark Heatmap Activated Multimodal Gaze Estimation

3D gaze estimation is about predicting the line of sight of a person in 3D space. Person-independent models for the same lack precision due to anatomical differences of subjects, whereas person-specific calibrated techniques add strict constraints on scalability. To overcome these issues, we propose a novel technique, Facial Landmark Heatmap Activated Multimodal Gaze Estimation (FLAME), as a way of combining eye anatomical information using eye landmark heatmaps to obtain precise gaze estimation without any person-specific calibration. Our evaluation demonstrates a competitive performance of about 10% improvement on benchmark datasets ColumbiaGaze and EYEDIAP. We also conduct an ablation study to validate our method.

Introduction

Eye gaze is an important non-verbal cue of humans and is capable of describing human emotion and behaviour. It is used in large cross-domain applications such as study of human behavior [24], clinical diagnosis [START_REF] Vidal | Wearable eye tracking for mental health monitoring[END_REF], humancomputer and human-robot interaction [START_REF] Goldberg | Eye tracking in web search tasks: Design implications[END_REF][START_REF] Michalowski | A spatial model of engagement for a social robot[END_REF]. To have precision in such applications, it is very important to have an accurate estimation of eye gaze. 3D gaze estimation is about inferring the actual line-of-sight of a person in 3D space and is the main focus of this paper.

Owing to the importance of study of this problem, vast research has been conducted. With the evolution in the field of deep learning, this problem has also moved significantly towards using it to predict gaze [START_REF] Fischer | RT-GENE: Real-Time Eye Gaze Estimation in Natural Environments[END_REF][START_REF] Yu | Improving few-shot userspecific gaze adaptation via gaze redirection synthesis[END_REF]33,[START_REF] Zhang | Mpiigaze: Real-world dataset and deep appearance-based gaze estimation[END_REF]. Many techniques focus on person-independent gaze estimation [33,[START_REF] Zhang | Mpiigaze: Real-world dataset and deep appearance-based gaze estimation[END_REF] (Figure 1(a)) where the model is evaluated on un-encountered subjects. But gradually, it was proved that it is not possible to reduce error below a certain point due to heavy dependency of gaze on the anatomy of the eye of different subjects [START_REF] Guestrin | General theory of remote gaze estimation using the pupil center and corneal reflections[END_REF]. Few classical model-based techniques [START_REF] Guestrin | General theory of remote gaze estimation using the pupil center and corneal reflections[END_REF][START_REF] Hansen | In the eye of the beholder: A survey of models for eyes and gaze[END_REF] were able to incorporate these anatomical features with few calibration samples, but they are not robust to variation in the image under different settings. As a result, subject-specific personalized models [START_REF] Yu | Improving few-shot userspecific gaze adaptation via gaze redirection synthesis[END_REF]22] (Figure 1(b)) got introduced which were able to further reduce errors using additional calibration steps on the test subjects to incorporate anatomical features.

However, although these person-specific models tend to give a low value of error, using them to design practical systems is difficult. This is because, in real-world scenarios, systems keep on encountering new subjects. So, for additional calibration, there have to be steps to take test samples with accurate ground truth and feed it to the model to calibrate before the actual task can be carried out. In different settings with varying surroundings and distance between subject and camera, it is challenging to procure test samples with accurate ground truth to calibrate the model. Further, for applications like clinical diagnosis of psychological dis-orders, for example autism and schizophrenia, where gaze is a very important non-verbal cue, the affected people cannot follow precise instructions to be performed for acquiring calibration samples. This makes it practically impossible to rely on person-specific methods for many applications.

In order to address this issue, and also use the anatomical features while predicting the gaze, we propose the Facial Landmark Heatmap Activated Multimodal Gaze Estimation (FLAME), which focuses on precise gaze estimation without calibration by incorporating the anatomical information contained in eye landmarks along with RGB features (Figure 1(c)). Landmarks corresponding to the outline of the eye hold key information about its anatomical features, and as also shown by [29], depict strong correlation with gaze direction. However, since the extraction of eye landmarks using other algorithms are also prone to errors, to make the model robust to these, we take a heatmap based approach inspired by [START_REF] Duan | Revisiting skeleton-based action recognition[END_REF] for skeleton-based action recognition. Instead of using absolute coordinates of landmark points, we use a Gaussian probability distribution heatmap for each landmark. This allows us not to pose strong constraints on the obtained value of the landmark coordinates and rather to provide it in the form of a probability distribution, while also allowing the system to extract salient anatomical features from it.

The high-level idea is depicted in Figure 2, which follows a two-stream CNN based neural network, working on RGB and eye landmark heatmap modalities to predict gaze angles. These angles can easily be converted to directional vector using basic trigonometry. To the best of our knowledge, this is the first work in gaze estimation using the information from eye landmarks in the form of a heatmap along with the RGB image.

Key contributions of this work are as follows:

• A novel approach for precise gaze estimation without calibration by incorporating anatomical features of the eye seamlessly with only RGB image as input.

• A method for obtaining, and using eye landmarks as heatmaps to extract anatomical features and merge with RGB feature map.

• Extensive experimentation on public benchmarks to obtain accurate gaze estimations in personindependent settings.

In the rest of the paper, we discuss our proposed approach in detail. Section 2 describes the past work in this field, followed by complete details of the proposed methodology in Section 3. In Section 4 we discuss the implementation details, which is succeeded by results and ablations in Section 5. Section 6 concludes the work and discusses its future scope.

Related Work

3D gaze estimation can be broadly divided into 2 parts: geometric methods [START_REF] Hansen | In the eye of the beholder: A survey of models for eyes and gaze[END_REF] and appearance-based methods [START_REF] Tan | Appearance-based eye gaze estimation[END_REF]. Geometric methods [START_REF] Guestrin | General theory of remote gaze estimation using the pupil center and corneal reflections[END_REF][START_REF] Hansen | In the eye of the beholder: A survey of models for eyes and gaze[END_REF] focus on modelling the eye, extracting features in terms of geometrical parameters and use it to predict gaze. They have several constraints like near-frontal head-pose and high-resolution images, and thus have limited application as compared to appearancebased methods, which can directly estimate gaze from eye images, and thus have gained significant popularity in present days.

Among person-independent methods, [33] proposed a LeNet-based architecture using a single eye patch, [START_REF] Fischer | RT-GENE: Real-Time Eye Gaze Estimation in Natural Environments[END_REF] used parallel VGG-16 based architecture and fed eyes and face to estimate head pose and gaze angle, [4] used Dilated Convolutions, while [START_REF] Zhang | Mpiigaze: Real-world dataset and deep appearance-based gaze estimation[END_REF] used a more complex CNN backbone to predict gaze. [34] did full-face gaze estimation by generating an attention map, and [START_REF] Krafka | Eye tracking for everyone[END_REF] implemented three parallel CNN backbones to process left eye, right eye and full face to predict gaze. [START_REF] Ranjan | Light-weight head pose invariant gaze tracking[END_REF] depicted a more robust use of head-pose for gaze estimation. Further, [29] proposed a multi-task CNN to predict eye-gaze and facial landmarks, and concluded that eye-landmarks are strongly co-related to gaze. As a developing trend, to remove person-specific bias of the models as shown by [START_REF] Guestrin | General theory of remote gaze estimation using the pupil center and corneal reflections[END_REF], person-specific models [START_REF] Yu | Improving few-shot userspecific gaze adaptation via gaze redirection synthesis[END_REF]22] have been proposed for gaze estimation which calibrate the model on few samples of test subjects before evaluation. Unsupervised learning has also come into the picture with [START_REF] Yu | Unsupervised representation learning for gaze estimation[END_REF] proposing unsupervised gaze estimation model using gaze redirection mechanism. Recently, [6] proposed gaze estimation using transformers and [START_REF] Biswas | Appearance-based gaze estimation using attention and difference mechanism[END_REF] implemented gaze estimation by attention mechanism and also using difference layer to remove unwanted features from both eyes.

Multimodal fusion has been of significant interest across domains. With time, late fusion has gained significant popularity where each mode is processed separately and combined later. It can be achieved by element-wise summation, concatenation, bi-linear product [3], weighted average [21], rank minimization [START_REF] Deng | Retinaface: Single-shot multi-level face localisation in the wild[END_REF]. [START_REF] Jacobs | Adaptive mixture of local expert[END_REF] involves using attention to pick the best mode for each input. Multimodal fusion module [16] was proposed which allows slow modality fusion in features of different dimensions.

Inspired by these works, we developed our own technique which we describe in-depth in the following sections.

Method Overview

The first part of our method is to extract the eye landmarks from the RGB image and generate the heatmap, followed by a two stream network having multiple components. From here on, for all discussions, r refers to RGB image and h represents eye landmark heatmap.

Eye Landmark Heatmap

To get the heatmap, we first extract the 2D eye landmarks from corresponding RGB images using facial landmark submodule [START_REF] Baltrusaitis | Constrained local neural fields for robust facial landmark detection in the wild[END_REF][START_REF] Zadeh | Deep constrained local models for facial landmark detection[END_REF] of OpenFace 2.0 [START_REF] Baltrusaitis | Openface 2.0: Facial behavior analysis toolkit[END_REF]. We prefer this because it can provide landmarks for the outline of the cornea and pupil along with outer eye, and is significantly precise and easy to use. OpenFace gives 28 2D landmarks for each eye as shown in Figure 3 in pixel coordinates, capturing the outline of the complete eye, cornea and pupil. These three components not only allow us to capture the shape of the outer eye but also the boundaries of inner regions, giving a large variety of descriptive information about the anatomy of the eye.

P (X) = 1 (2π) d |Σ| exp - 1 2 (X -µ) Σ -1 (X -µ) (1)
with mean µ and covariance matrix Σ. We use it to obtain the Gaussian probability distribution heatmap of the eye landmark coordinates having the same dimensions as the input image.

Transfer Function

The transfer function T (r, h) is the function that defines the exchange of information between the two streams. It is based on Multimodal Transfer Module (MMTM) [16], which is a slow modality fusion block used to re-calibrate channel-wise features based on squeeze and excitation mechanism between any two feature maps of arbitrary dimension. The complete set of steps can be visualized from Figure 4. Assume an input feature vector of height n and width m. Let r i ∈ R n×m×cr i be the i-th RGB feature map of c ri channels and h i ∈ R n×m×c h i the i-th heatmap feature vector of c hi channels. We first squeeze them into S r i ∈ R cr i and S h i ∈ R c h i , the channel-wise descriptors, using global average pooling across the spatial dimensions:

S ri (c) = 1 n • m n p=1 m q=1 r i (p, q, c) (2)
S hi (c) = 1 n • m n p=1 m q=1 h i (p, q, c) (3)
Then, we concatenate these two, and pass it through a dense layer to obtain a joint representation Z i . This Z i is passed through separate dense layers to obtain excitation signals E r i ∈ R cr i and E h i ∈ R c h i for both RGB and heatmap modalities. This hybrid feature map along with head pose is passed to a fully-connected regression network to predict gaze angles.

The final re-calibrated feature map can be obtained by taking sigmoid σ of the excitation signal and performing channel-wise multiplication with the actual feature map:

r i (c) = 2 • σ (E ri (c)) • r i (c) (4)
h i (c) = 2 • σ (E hi (c)) • h i (c) (5)
We take dim(Z i ) = 1 4 (dim (S r i ) + dim (S h i )), as suggested in [16]. In our case, dim (S r i ) = dim (S h i ).

To design the complete transfer function T (r, h), we use this MMTM block for feature reactivation at the last and second last feature maps, as given in Figure 5. This is because we want the network to learn some initial representation and then, the higher-level features of both streams can be utilized by each other. At the deeper layers, using this transfer function can help both individual unimodal streams to learn better representation by benefiting from the information extracted by the other stream.

CNN Backbone

We believe that gaze is a relatively low-level feature and thus, a shallow architecture is effective to learn relevant features. So, we designed a custom CNN backbone based on Residual Blocks [START_REF] He | Deep residual learning for image recognition[END_REF] as given in Figure 5.

It starts with a 3×3 convolution layer followed by Batch Normalization, ReLU and maxpool with kernel size (2, 2) and stride (2, 2). This is followed by 3 modules, each having 2 residual blocks followed by a maxpool layer. The intermediate feature output at the end of each module has 64, 128 and 256 channels respectively. Both RGB and heatmap modalities are processed using identical backbones to maintain symmetry, as it creates identical feature dimensions to ease and balance other operations in the network.

Aggregation Function

The aggregation function A (r, h) is used to fuse the final feature map of both streams in order to get a hybrid representation. This is done by channel-wise concatenation of the individual feature maps, and then passing this through a residual block. The output of this is a hybrid feature map having relevant information from both RGB and heatmap modalities.

FLAME

The complete architecture of FLAME, as depicted in Figure 5, is designed from the previously described individual units. The two separate CNN backbones take 120×120 RGB eye patch and 120×120×28 eye landmark heatmap. Both learn features independent of each other until the last two modules where both the features are fused into one another by the transfer function T (r, h), which consists of two MMTM blocks. Thereafter, using the two final feature maps, a hybrid feature map is obtained by the aggregation function A (r, h). This feature map contains the information of the gaze direction as well as the anatomical information of the eye and is a better descriptor of the gaze direction.

It is then flattened and concatenated with the head pose angles and is passed through two dense layers of 512 neu-rons each, followed by a two-dimensional output layer for the pitch and yaw gaze angles directly in the World Coordinate System (WCS). In the dense layers, dropout is used with p dropout = 0.2 and ReLU is used as activation function. The output layer, however, is kept as linear with no activation function.

Loss Function

This model is end-to-end trainable using the gaze angle ground truth. The 3D angular loss can be defined as

L ang ( ĝp , ĝt ) = arccos ĝp . ĝt | ĝp || ĝt | , (6)
where ĝp and ĝt represent the predicted and true gaze vectors respectively. This loss, given by Equation 6, directly gives the angular deviation of predicted gaze vector ĝp from the true gaze vector ĝt in 3D space. However, training directly on this loss function, as used by a few models [START_REF] Yu | Improving few-shot userspecific gaze adaptation via gaze redirection synthesis[END_REF], is not ideal because the gradient of this loss, given by

∂L ang (x) ∂x = 1 √ 1 -x 2 where x = ĝp . ĝt | ĝp || ĝt | , (7)
increases sharply with decrease in loss and approaches ∞ as the loss approaches 0 (x approaching 1). So, training on this loss makes the network unstable with decrease in loss. Thus, we define a vector difference loss for training, which is the sum of squared errors of each individual component of the gaze vector g k = (g x k , g y k , g z k )

where k = {p, t}, given by

L vec ( ĝp , ĝt ) = g x p -g x t 2 + g y p -g y t 2 + g z p -g z t 2 . (8)
For evaluation, we use the 3D angular loss itself as per Equation 6.

Implementation Details

Datasets

We perform the experiments with two benchmark datasets: ColumbiaGaze [START_REF] Smith | Gaze Locking: Passive Eye Contact Detection for Human?Object Interaction[END_REF], EYEDIAP [START_REF] Mora | Eyediap: a database for the development and evaluation of gaze estimation algorithms from rgb and rgb-d cameras[END_REF]. EYEDIAP is a video dataset with 16 subjects, with different head motion and gaze target settings. We select every 20-th frame of the HD screen target videos (continuous and discrete) under both static and moving head pose of all subjects for our experiments. ColumbiaGaze is a dataset with discrete values of head pose and targets, having 56 subjects, each with a combination of 7 horizontal and 3 vertical gaze targets. We select all the images from this dataset. For all these images, the ones on which OpenFace outputs could not be obtained were left out of the final dataset. For our experiments, we use the raw images without any re-alignment.

Preprocessing

We first extract face crops of the ColumbiaGaze and EYEDIAP images using RetinaFace [8], and zero-pad it to a uniform size of 384×480 (4:5 ratio). As ColumbiaGaze dataset was of higher resolution, we rescale the padded 4:5 crops to this dimension. Then, we run OpenFace 2.0 to obtain eye landmarks, and generate the heatmap for all the data with dimension d = 2, Σ = I 2 (2D identity matrix) and µ = (x i , y i ) , i = 1, 2, . . . , 28 where (x i , y i ) are the obtained eye landmark coordinates, following the method discussed in Section 3.1. We then divide the complete dataset into cross-subject split of 8:1:1 for train, validation and test sets, respectively.

Training and Optimization

The network is trained with randomly initialized weights for 200 epochs with a batch size of 8, using Adam Optimizer with an initial learning rate of 10 -4 , and is decreased by a factor of 0.5 after 85-th, 120-th and 175-th epoch. During training, 120×120 pixels eye patch of either left or the right eye is picked randomly for each image sample along with the corresponding landmark heatmap of that eye having same resolution, and fed to the network. The eye crop is obtained by calculating the centre of the eye using the four corner landmarks (8, 14, 11, 17 in Figure 3), and taking a 120×120 patch across it with the obtained centre. The same operation is done with landmark heatmap to maintain alignment between both inputs. Head pose is used as available with the datasets, but can also be obtained from facial analysis toolkits like OpenFace 2.0 with sufficient precision. The model is then trained against ground truth gaze angles using the vector difference loss given by Equation 8.

Experimental Settings

In order to conduct an ablation study of the proposed method, we implemented three other models in addition to our main approach of FLAME. All four experimental settings are described below: FLAME: Our original approach as proposed in Section 3.5 with settings as described in Section 4.3 on all the datasets. FLAME-aggregation-only (F-AO): The original architecture without the transfer function T (r, h). This is to evaluate the impact of using multimodal transfer. FLAME-additive-fusion (F-AF): Another approach to fuse both modalities, where aggregation function is the element-wise addition of the final feature map of both modalities. This is to analyse the effectiveness of aggregation function A (r, h). There is no T (r, h) here. FLAME-baseline (F-B): This is the baseline approach that performs gaze estimation using only the RGB images and proposed CNN backbone. It is to see the contribution of the entire heatmap modality and proposed novelty.

Results

The performance of our method can be evaluated by Equation 6, the angular error between true and predicted gaze vectors in 3D space. Results are given in Table 1 for both cross-subject and cross-dataset validation to see how well the model performs on unencountered subjects (cross-subject evaluation), and also on completely unknown setup (distance of subjects, camera settings, etc.) and range of ground truth (cross-dataset evaluation). In Table 2 we can see that when the model was not provided with any previously encountered subjects, it was able to extract the information from eye landmarks and reduce the mean error by 1.29 • and 0.7 • as compared to RGB only representing the baseline (F-B). However, to still analyse the factors contributing to errors, it can be understood that we are using 2D facial landmarks, and not 3D. This might fail to capture the precise orientation and shape of the eye in 3D space, and may not describe accurate anatomical structure. Additionally, there can be significant error in eye landmark extraction, particularly when the eye is not fully visible, or when the person is looking directly into the cam-era. These can be larger for the cornea and pupil landmarks which have finer boundaries.

Cross-dataset error is predominantly because the two datasets are very different to one another in many parameters such as range of gaze angles, distance of subject, or camera parameters. However, although these values are not directly comparable to other methods due to difference in datasets used for cross-dataset evaluation, the performance is competitive to them.

Ablation Study

It is also important to understand the contribution of all the components proposed in our methodology. For this, we implemented three other architectures as described in Section 4.4. The results obtained are given in Table 2. We can see that FLAME outperforms all other experimental settings. Figure 7 depicts the distribution of error for EYEDIAP and ColumbiaGaze datasets on all four experimental settings through a box plot. There is a difference in mean error of 0.18 • for EYE-DIAP and 0.42 • for ColumbiaGaze between F-AO and FLAME, where the absence of transfer function was the only difference. This shows its contribution.

Further, F-AO and F-AF differ in aggregation function A (r, h). We compare our approach with element-wise addition method in F-AF, where A (r, h) = r 4 + h 4 . The error difference highlights the effectiveness of our aggregation function for creating a hybrid feature map.

At last, the contribution of our complete methodology can be seen through the difference of FLAME with F-B.

We also compared the performance of FLAME with different input resolutions. In practical settings, having highresolution eye images is a constraint due to limited camera resolution, distance of subject from it, etc. On Columbi-aGaze, we obtained error of 4.79 • in 60×60 pixels and 5.5 • on 30×30 pixels image resolution. Based on this, we would propose 120×120 pixels as optimum and 60×60 pixels as the minimum threshold of input resolution with FLAME.

Comparison with State-of-the-art

We also compare our method with various state-of-theart. The results are given in Table 3. All of them are for leave-one-subject-out settings (or cross-subject evaluation). It can be seen that our method performs significantly better than all of the person-independent methods for gaze estimation. This shows the impact of using eye landmark heatmaps to utilize anatomical information of the eye along with the RGB image. When comparing to the models of [START_REF] Yu | Unsupervised representation learning for gaze estimation[END_REF], a person-specific method (based on few-shot calibration), our approach performs better than all 3 methods on both datasets, except for one model in ColumbiaGaze. A larger error than person-specific models can be because when the model gets the test samples directly, it can extract additional information like orientation in 3D space, PERCLOS (% of eye closed), and other salient anatomi-cal features, which are not possible to capture using 2D eye landmarks. For ColumbiaGaze, SOTA benchmarks are not available for many models.

Thus, with all these analyses and experimentation, we can justify that our approach holds key to a precise gaze estimation.

Conclusion

We presented an approach of using 2D eye landmarks along with RGB images to perform 3D gaze estimation by incorporating anatomical features of the eye without any person-specific calibration. In our approach, we obtained a 2D Gaussian probability distribution heatmap of eye landmarks from the RGB image itself and then used a twostream CNN-based network to effectively extract relevant features. We used this along with head pose to predict gaze angles. Our method gives better results than all state-of-theart methods on two benchmark datasets. This proves that eye landmarks can play a vital role in incorporating anatomical information to predict gaze more accurately.

We believe more precise systems can be designed if we can incorporate more information related to the anatomical aspects of the eye, and use eye landmarks in other ways. As of now, we leave that for further study.

Figure 1 :

 1 Figure 1: Different types of gaze estimation methods: (a) person-independent technique, (b) person-specific technique, (c) FLAME. Training subjects are green and test subjects are yellow. r stands for RGB image and h stands for eye landmark heatmap.

Figure 2 :

 2 Figure 2: Method Summary. Landmark heatmap is first obtained from cropped eye images. Then, both are processed through a two-stream CNN-based network to extract gaze and anatomical information separately, which is further used along with head pose to predict gaze. Dotted line marks the pre-training extraction and solid line marks the training pipeline.

Figure 3 :

 3 Figure 3: Eye landmarks from Openface 2.0

Figure 4 :

 4 Figure 4: Multimodal Transfer Module [16]

Figure 5 :

 5 Figure 5: Complete model architecture. Two identical CNN backbones process the RGB and heatmap modality separately, with exchange of information at intermediate layers as per transfer function and later fusing them using aggregation function.This hybrid feature map along with head pose is passed to a fully-connected regression network to predict gaze angles.

Figure 6 :

 6 Figure 6: Analysis of predictions versus true labels. Top heatmaps are evaluating the yaw angle and bottom heatmaps are evaluating the pitch angles.

Figure 7 :

 7 Figure 7: Distribution of errors of all experimental settings

Table 1 :

 1 Results on cross-subject, cross-dataset evaluation.

	Train	Test ColumbiaGaze EYEDIAP
	ColumbiaGaze	4.64 •	12.53 •
	EYEDIAP	12.83 •	4.62 •
	Figure 6 represents the relation between prediction and
	truth values for both datasets. They pretty much follow a
	linear relationship which further proves the correctness of
	our methodology.

Table 2 :

 2 Mean and standard deviation (std) of angular error for different experimental settings

	Experiment Setup ColumbiaGaze	EYEDIAP
		mean	std	mean	std
	F-B	5.93 • 3.20 •	5.32 • 3.08 •
	F-AF	5.88 • 3.06 •	5.30 • 3.03 •
	F-AO	5.06 • 3.13 •	4.80 • 3.02 •
	FLAME	4.64 • 2.86 •	4.62 • 2.93 •

Table 3 :

 3 Mean 3D angular error as reported by different SOTA methods: person-specific methods (top), personindependent methods (bottom) error metric as reported by[START_REF] Cheng | A coarseto-fine adaptive network for appearance-based gaze estimation[END_REF], # error metric as reported by[6]

	Method	ColumbiaGaze EYEDIAP
	Yu et al. (sup. HCS) [31]	5.25 •	7.09 •
	Yu et al. (sup. WCS) [31]	3.54 •	6.79 •
	Yu et al. (unsup.) [31]	7.15 •	8.2 •
	Zhang et al. [34]	-	6.0 •
	Zhang et al. [33] #	-	7.37 •
	Zhang et al. [35] #	-	6.79 •
	Cheng et al. [5]	-	5.3 •
	Cheng et al. [6]	-	5.17 •
	Fischer et al. [10] *	-	6.4 •
	Chen et al. [4] *	-	5.9 •
	Kellnhofer et al. [17] #	-	5.36 •
	Yu et al. [29]	-	8.54 •
	FLAME (ours)	4.64 •	4.62 •

*


Fredrik Johansson

RAPID COMPUTATION OF SPECIAL VALUES OF DIRICHLET L-FUNCTIONS

Keywords: 2020 Mathematics Subject Classification. Primary 33F05, 11M06; Secondary 33B20, 65D20

We consider computing the Riemann zeta function ζ(s) and Dirichlet L-functions L(s, χ) to p-bit accuracy for large p. Using the approximate functional equation together with asymptotically fast computation of the incomplete gamma function, we observe that p 3/2+o(1) bit complexity can be achieved if s is an algebraic number of fixed degree and with algebraic height bounded by O(p). This is an improvement over the p 2+o(1) complexity of previously published algorithms and yields, among other things, p 3/2+o(1) complexity algorithms for Stieltjes constants and n 3/2+o(1) complexity algorithms for computing the nth Bernoulli number or the nth Euler number exactly.

Introduction

Let χ be a Dirichlet character modulo q ≥ 1. The associated Dirichlet L-function is the analytic continuation of

L(s, χ) = n=1 χ(n) n s , Re(s) > 1 (1)
to s ∈ C with the possible exception of a pole at s = 1. The Riemann zeta function ζ(s) is the Dirichlet L-function corresponding to the trivial character χ(n) = 1, which is the unique character modulo q = 1. If χ is a primitive character, then the function L(s, χ) is represented in the entire complex plane by a convergent expansion, the approximate functional equation [Coh19, Theorem 7.3] Γ s + δ 2 L(s, χ) = δ q,1 π s/2 α (s-1)/2 s -1

- α s/2 s + ∞ n=1 χ(n) n s Γ s + δ 2 , πn 2 α q + ω π q s-1/2 ∞ n=1 χ(n) n 1-s Γ 1 -s + δ 2 , πn 2 αq (2)
where δ ∈ {0, 1} is the parity χ(-1) = (-1) δ , ω is the root number of χ, which satisfies |ω| ≤ 1, and Γ(a, z) = ∞ z t a-1 e -t dt is the incomplete gamma function. The quantity α is a free positive parameter; we may take α = 1 to balance the rate of convergence of both series. If χ is not primitive, we can decompose L(s, χ) in terms of primitive functions.

The expansion (2) is useful for high-precision computation of L(s, χ) due to the super-exponential decay Γ(a, Cn 2 ) ≈ exp(-Cn 2 ) of the incomplete gamma functions. We need only O(p 1/2 ) terms for a desired bit precision p, which should be contrasted with Euler-Maclaurin summation [BC21, §4.2] [Joh14b] and similar methods which require O(p) terms. (An equally important advantage of (2) is that we only need O(q 1/2 ) terms as a function of the modulus q.)

The drawback of ( 2) is that we have to compute the nonelementary incomplete gamma functions. Our goal is to study this problem with attention to the bit complexity when p → ∞. We use "time" and "bit operations" synonymously, and recall that floating-point numbers with p-bit precision can be multiplied in time O(p log p) = p 1+o(1) [START_REF] Harvey | Integer multiplication in time O(n log n)[END_REF]. Using Euler-Maclaurin summation, for example, it is easy to show that we can compute L(s, χ) or any of its s-derivatives L (j) (s, χ) to p-bit accuracy in time p 2+o(1) for fixed s, χ and j. Our main observation is the following improved complexity bound for special values s.

Theorem 1. Let χ be a fixed Dirichlet character, and let s ∈ Q be an algebraic number of fixed degree such that the minimal polynomial of s over Z has coefficients bounded in absolute value by O(p). Then, for any fixed j ≥ 0, the value L (j) (s, χ) can be approximated with absolute error less than 2 -p in time p 3/2+o(1) using p 1+o(1) space. (When s = 1 and this point is a pole, the corresponding Laurent series coefficient is computed.)

Proof. The height condition implies that |s| = O(p). Since any terms and prefactors appearing in (2) and in the asymptotics of the incomplete gamma function are bounded by exp(|s| 1+o(1) ), it is sufficient to truncate both infinite series to N = p 1/2+o(1) terms and approximate the terms to p 1+o(1) bits.

The function y(z) = Γ(a, z) is holonomic, satisfying a linear differential equation Ay = 0 with A ∈ Z[a, z, d dz ]. We evaluate y (for two values of the parameter a) at N points z. We can compute each such function value in p 1+o(1) bit operations and p 1+o(1) space using the bit-burst algorithm [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], employing arithmetic in the number field Q(s). This bound holds uniformly for the required values of a and z, by the same argument as in [Mez12, Corollary 1], using the facts that |z| ≤ p 1+o(1) and that A as well as the defining polynomial of Q(s) have fixed degree and coefficients of bit size O(log p).

Using standard methods, the remaining operations (evaluation of Dirichlet characters, the gamma function, and elementary functions) fall within the same complexity bounds.

For derivatives L (j) (s, χ), and at removable singularities, the equivalent operations can be carried out using arithmetic on truncated formal power series.

We will provide additional details below. The only interesting point in the proof of Theorem 1 is the use of the bit-burst algorithm instead of naive summation, which allows us to compute the function Γ(a, z) in quasilinear rather than quadratic time. The bit-burst algorithm for holonomic functions has been known since the 1980s [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF][START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions D-finies[END_REF][START_REF] Mezzarobba | A note on the space complexity of fast D-finite function evaluation[END_REF], and since the 1970s in special cases [START_REF] Brent | The complexity of multiple-precision arithmetic[END_REF], yet we are not aware of a correct complexity bound of this kind in the literature for Dirichlet L-functions or even for the special case of the Riemann zeta function.

The use of the approximate functional equation for L-function computation has been studied in detail by Rubinstein [START_REF] Oded | Evidence for a spectral interpretation of the zeros of Lfunctions[END_REF] and several other authors [START_REF] Dokchitser | Computing special values of motivic L-functions[END_REF][START_REF] Andrew | Artin's conjecture, Turing's method, and the Riemann hypothesis[END_REF][START_REF] Molin | Intégration numérique et calculs de fonctions L[END_REF][START_REF] Belabas | Numerical Algorithms for Number Theory: Using Pari/GP[END_REF]. These works do not mention the bit-burst algorithm or address the bit complexity for large p, instead focusing on parameters relevant for numerical testing of the generalized Riemann hypothesis, namely with fixed p and with large | Im(s)| and/or large q. For large | Im(s)|, one must use a "smoothed" version of (2) to avoid exponentially large cancellation, or Riemann-Siegel type expansions; the method in Theorem 1 is not competitive in this setting, where the best methods achieve O(| Im(s)| 1/2 ) or lower complexity for a fixed level of accuracy.

Borwein, Bradley and Crandall [START_REF] Borwein | Computational strategies for the Riemann zeta function[END_REF] and Crandall [START_REF] Crandall | Unified algorithms for polylogarithm, L-series, and zeta variants[END_REF][START_REF] David | Crandall's computation of the incomplete gamma function and the Hurwitz zeta function, with applications to Dirichlet L-series[END_REF] discuss the approximate functional equation in the context of high-precision zeta function computation, but do not mention the bit-burst algorithm or give a complexity bound of this type. Crandall [START_REF] Crandall | Unified algorithms for polylogarithm, L-series, and zeta variants[END_REF] writes that the incomplete gamma function can be computed using p 1+o(1) "operations", but this is referring to full-precision arithmetic operations, which would give us p 2+o(1) bit complexity for Γ(a, z) and p 5/2+o(1) bit complexity for L-functions. A refined algorithm that achieves p 2+o(1) bit complexity is described in [BBC00, §7]; see §5 below.

In a 1988 paper, Borwein and Borwein [START_REF] Borwein | On the complexity of familiar functions and numbers[END_REF] claim that ζ(s) can be computed in time p 1+o(1) if s is a fixed rational number, and in time p 3/2+o(1) if s is a fixed generic (computable) complex number. No explicit algorithm is given to justify these claims: the authors simply write "we truncate both the integral and the sum" with reference to the formula

ζ(s)Γ s 2 π -s/2 - 1 s(s -1) = ∞ 1 t (1-s)/2 + t s/2 t ∞ n=1 e -n 2 πt dt (3)
which, apart from minor differences in notation, is the approximate functional equation for ζ(s) in the form originally derived by Riemann [START_REF] Riemann | Ueber die Anzahl der Primzahlen unter einer gegebenen Grosse[END_REF] (we obtain the series in incomplete gamma functions by integrating term by term). Both bounds claimed by the Borweins are a factor p 1/2 better than all methods known to this author. Lacking evidence to the contrary, we believe that the Borweins had in mind some combination of the algorithms that will be described below and that their complexity analysis was erroneous. Our goal with this article is therefore in part to correct the record. 1 The rest of this paper is structured as follows: §2 discusses some consequences of Theorem 1, §3 gives a more detailed description of the algorithm, and §4 reports implementation results. Finally, §5 discusses alternative algorithms for use when s is not algebraic.

Applications and generalizations

The immediate application of Theorem 1 is that p 3/2+o(1) can be a significant improvement over p 2+o(1) for numerical evaluation to tens of thousands of digits. Such computations are not exclusively done to test algorithms; for example, integer relation searches employing 50,000-digit precision have been successful in discovering new identities involving special values of L-functions [START_REF] Bailey | Experimental mathematics: Recent developments and future outlook[END_REF].

2.1. Values at integers. The most famous special values, ζ(n) and L(n, χ) with n ∈ Z, can be expressed in terms of logarithmic derivatives of the gamma function at rational points, and can consequently be computed in quasilinear time p 1+o(1) using binary splitting [START_REF] Karatsuba | Fast evaluation of the Hurwitz zeta function and Dirichlet L-series[END_REF][START_REF] Johansson | Arbitrary-precision computation of the gamma function[END_REF]. 1 We mention that J. Borwein coauthored the 2000 survey paper [BBC00] on ζ(s) computation, which discusses the approximate functional equation prominently but does not mention the claims from 1988. This omission suggests that the Borweins were aware of the error (but perhaps did not consider it important enough to publish a correction). The other complexity results in [START_REF] Borwein | On the complexity of familiar functions and numbers[END_REF], for instance concerning Γ(s), are correct.

Alternatively (and often more efficiently), binary splitting can be applied directly to convergence-accelerated series for the Riemann zeta function [Joh14a, §4.7] or hypergeometric series for particular values such as

ζ(3) = 5 2 ∞ n=1 (-1) n+1 (n!) 2 n 3 (2n)! (4)
which have been used to compute billions of digits [START_REF] Sebah | The Apery's constant: ζ(3)[END_REF][START_REF] Yee | y-cruncher -a multi-threaded pi-program[END_REF]. However, these quasilinearity results all assume that n = p o(1) . For example, if n and p are proportional (or proportional up to logarithmic factors), then the complexity degenerates to p 2+o(1) or worse. The complexity is also softly quadratic with n ∝ p if we compute L(s, χ) directly using the L-series (1) or the corresponding Euler product. The p 3/2+o(1) complexity of Theorem 1 is then an improvement over previous algorithms.

2.2.

Bernoulli and Euler numbers. The Bernoulli numbers and Euler numbers are the rational numbers and integers respectively defined by

x e x -1 = ∞ n=0 B n n! x n , 1 cosh(x) = ∞ n=0 E n n! x n . (5)
The odd-index values are trivial, while the even-index values can be expressed in terms of Dirichlet L-functions as

B 2n = (-1) n+1 2(2n)! (2π) 2n ζ(2n), E 2n = (-1) n 4 n+1 (2n)! π 2n+1 β(2n + 1) (6)
where β(s) = L(s, χ 4.3 ) is the Dirichlet beta function, corresponding to the character modulo q = 4 with (χ 4.3 (n)) ∞ n=0 = (0, 1, 0, -1, . . .). There are Θ(n log n) bits in E n and in the numerator of B n (the denominator is easy to determine), so we can recover the exact values by evaluating the L-functions numerically to p = n 1+o(1) bits. As a corollary of Theorem 1, we have the following:

Theorem 2. The nth Bernoulli number B n and Euler number E n can be computed exactly in time n 3/2+o(1) using n 1+o(1) space.

The significance of this result is that all methods known until quite recently (for instance those employing the Euler product) require at least n 2+o(1) time.

Harvey [START_REF] Harvey | A subquadratic algorithm for computing the n-th Bernoulli number[END_REF] gave the first subquadratic algorithm for computing B n , which uses n 4/3+o(1) time and space or n 3/2+o(1) time when confined to n 1+o(1) space. We fail to improve on Harvey's bound, but the methods are independent: ours is numerical; Harvey's uses modular arithmetic and does not involve L-functions. For Euler numbers, no subquadratic algorithm has been published before ours, though it is plausible that Harvey's algorithm can be generalized to this case. 2.3. Sparse zeta-expansions. Many slowly converging series or products can be evaluated to high precision using zeta function acceleration [START_REF] Flajolet | Zeta function expansions of classical constants[END_REF], which is based on the formal rearrangement

n∈A f (1/n) = m f m ζ A (m), f (z) = m f m z m , ζ A (s) = n∈A n -s . ( 7
)
The complexity of approximating such a sum to p-bit accuracy is typically p 2+o(1) provided that the transformed series converges geometrically and that O(p) coefficients f m and zeta values ζ A (m) can be evaluated simultaneously to p-bit accuracy in time p 2+o(1) using FFT-based power series operations (this is the case if f is elementary and ζ A (s) = ζ(s), for example).

If the resulting zeta-expansion is sparse, then Theorem 1 may yield an improved complexity bound. An example is the Landau-Ramanujan constant

λ =   1 2 p≡3 mod 4 1 1 -p -2   1/2 ≈ 0.764 (8)
which appears in the asymptotic formula λx/ log(x) for the number of integers k ≤ x expressible as a sum of two squares. Flajolet and Vardi [START_REF] Flajolet | Zeta function expansions of classical constants[END_REF] obtain the sparse zeta-type expansion

λ = 1 √ 2 ∞ n=1 1 - 1 2 2 n ζ(2 n ) β(2 n ) 1/2 n+1 (
γ = log(4π) -2 + 2 √ π ∞ n=1 1 n Γ( 1 2 , πn 2 ) + 2 ∞ n=1 Γ(0, πn 2 ) (10)
along with a much more complex formula for γ 1 written in terms of series of 2 F 2 and 3 F 3 hypergeometric functions and digamma functions. Coffey notes that these formulas "may have some attraction for computation" due to the e -πn 2 type decrease of the terms. He also notes that the method can be generalized to Dirichlet L-function analogs of Stieltjes constants. As indicated in the proof of Theorem 1, implementing (2) with power series arithmetic provides such a generalization without requiring the derivation of unwieldy formulas for the higher derivatives.

Similarly, we obtain a p 3/2+o(1) complexity algorithm for the Glaisher-Kinkelin constant e 1/2-ζ ′ (-1) , Keiper-Li coefficients, etc.

The approximate functional equation has been used in a somewhat different way by Keiper [Kei92] to compute Stieltjes constants and other series coefficients related to the Riemann zeta function, with worse complexity than Theorem 1; we revisit this topic in §5. The number ζ(1/2) ≈ -1.4603545 makes an interesting appearance in a formula in Ramanujan's lost notebook (see [START_REF] Andrews | Ramanujan's Lost Notebook, Part IV[END_REF]§8.3], where generalizations to other values of ζ(s) and L(s, χ) with s ∈ Q are discussed as well). For any α, β > 0 such that

αβ = 4π 3 , ∞ n=1 1 e n 2 α -1 = π 2 6α + 1 4 + √ β 4π ζ( 1 2 ) + ∞ n=1 cos( √ nβ) -sin( √ nβ) -e - √ nβ √ n(cosh( √ nβ) -cos( √ nβ) . (11)
There is a parallel to the free parameter in (2): by varying α, we can force either the left or the right series to converge faster at the expense of the other. Setting α = O(1/p) in Ramanujan's formula minimizes the total number of terms for a given precision p; this leads to an algorithm with p 2+o(1) bit complexity to compute ζ(1/2), comparable to Euler-Maclaurin summation and inferior to Theorem 1. 2.6. Hurwitz zeta-type functions. The method behind Theorem 1 is not restricted to "proper" L-functions. Crandall [START_REF] Crandall | Unified algorithms for polylogarithm, L-series, and zeta variants[END_REF][START_REF] David | Crandall's computation of the incomplete gamma function and the Hurwitz zeta function, with applications to Dirichlet L-series[END_REF] has given a formula analogous to (2) (Crandall calls this a "Riemann-splitting representation") for the Lerch transcendent, which is the analytic continuation of the series

Φ(z, s, a) = ∞ n=0 z n (n + a) s , |z| < 1. ( 12
)
Combining Crandall's expansion with bit-burst evaluation of the incomplete gamma function should lead to p 3/2+o(1) algorithms for the following:

• The Lerch transcendent Φ(z, s, a) with s ∈ Q, z, a ∈ C and its s-derivatives.

• The Hurwitz zeta function ζ(s, a) with s ∈ Q, a ∈ C and its s-derivatives.

• The generalized Stieltjes constants γ n (a) with a ∈ C.

• The polylogarithm Li s (z) with s ∈ Q, z ∈ C and its s-derivatives.

• The Barnes G-function G(z) with z ∈ C.

We have not checked the details of these computations (validity of analytic continuations, possible exceptional points, explicit error bounds, uniform complexity with respect to parameters), and we leave this for a future study.

As in the case of integer zeta values, p 1+o(1) algorithms are already available for the above functions in some more restricted cases, e.g. for ζ(n, a) with a ∈ Q.

The algorithm

Since the proof of Theorem 1 above is quite terse and the bit-burst algorithm for generic holonomic functions requires much more complicated machinery than in the specialized case of computing Γ(a, z), we give a more explicit description here.

We may rely on ball arithmetic [START_REF] Van Der Hoeven | Ball arithmetic[END_REF][START_REF] Johansson | Arb: Efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF], which means that explicit error bounds need to be derived only for the truncation errors in infinite series; asymptotic estimates suffice for choosing the floating-point precision.

3.1. The outer series. To bound the tails of the infinite series in (2), the following formulas may be used. Similar bounds can also be found in [START_REF] Oded | Evidence for a spectral interpretation of the zeros of Lfunctions[END_REF].

Lemma 3. For real z > 0 and complex a = σ + τ i, the order j ≥ 0 parameter derivative of the incomplete gamma function Γ(a, z) satisfies the bound

Γ (j,0) (a, z) = z a-1 log j (z)e -z ∞ 0 e -t 1 + t z a-1 1 + log(1 + t/z) log(z) j dt ≤ z σ-1 log j (z)e -z ∞ 0 exp -t + (σ -1)t z + jt z log(z) dt ≤ z σ-1 log j (z)e -z 1 -B j (σ, z) , B j (σ, z) = 1 z max(σ -1, 0) + j log(z)
(13) provided that B j (σ, z) < 1, and assuming that z > 1 if j ≥ 1.

We illustrate how to bound the zeroth derivative of the first of the two infinite series in (2):

Lemma 4. Assume that N, q ≥ 1, s = σ + τ i with σ, τ ∈ R, δ ∈ {0, 1}. Define C = (σ + δ)/2 and D = πα/q. If DN 2 > C -1, then ∞ n=N χ(n) n s Γ s + δ 2 , πn 2 α q ≤ ∞ n=N 1 n σ (Dn 2 ) C-1 e -Dn 2 1 -B 0 (C, Dn 2 ) = D C-1 1 -B 0 (C, DN 2 ) ∞ n=N e -Dn 2 n 2-δ ≤ D C-1 1 -B 0 (C, DN 2 ) e -DN 2 N 2-δ (1 -e -D )
.

(

) 14
We can bound the tails of the s-derivatives as follows: we expand the power series product n -(s+X) Γ((s + X + δ)/2, παn 2 /q) symbolically, apply the bound (13) for each coefficient, and compute geometric series bounds similar to those in (14).

The bound for the other series in (2) is identical but with C = (1 -σ + δ)/2 and D = π/(qα).

Implementation remarks.

There is no need to derive a closed formula for choosing N ; we can simply evaluate the bound (14) for N = 1, 2, . . . and stop when the error meets a target tolerance.

For optimal performance, we should compute a tight estimate of the number of bits that each term contributes to the final result and only compute to that precision locally. It is useful to note that for a ∈ R and z > 0,

log(Γ(a, z)) ≈ (a -1) log(z) -z a < z a(log(a) -1) a ≥ z. ( 15
)
gives an accurate order-of-magnitude estimate of the incomplete gamma function.

Evaluation of the incomplete gamma function.

The idea of the bit-burst algorithm is to analytically continue a holonomic function y using the Taylor series method for ODEs, following a path

z initial z 1 z 2 . . . (16)
that approaches the target point z exponentially and thus converges in O(log(p)) steps. For example, for our application we may choose z 1 = ⌊z/2 32 ⌋2 32 , z 2 = ⌊z/2 64 ⌋2 64 , z 3 = ⌊z/2 128 ⌋2 128 , . . ., where successive steps double the number of leading bits extracted from the binary expansion of z.2 At each step, the Taylor series can be evaluated using binary splitting, and the exponentially converging steps balance the bit sizes of the coefficients against the number of terms in each Taylor series so that the overall bit complexity is quasilinear in p.

3.2.1. Hypergeometric series. Let (a) n = a(a+1) • • • (a+n-1). For the first Taylor step z initial z 1 , we may choose z initial = 0. Here we have the hypergeometric series

Γ(a, x) = Γ(a) - x a e -x a ∞ n=0 x n (a + 1) n , x = z 1 , (17)
which is valid for all x > 0 when a ∈ {0, -1, -2, . . .}. If the series is truncated after N terms where N > -Re(a) -1 and |a

+ N + 1| > |x|, then ∞ n=N x n (a + 1) n ≤ |x| N |(a + 1) N | 1 1 -C , C = |x| |a + N + 1| (18)
Truncation bounds for derivatives of this series with respect to a can be obtained similarly; see [Joh19, Theorem 1].

At the poles of the gamma function, a limit computation is needed; this can be done using power series arithmetic or explicitly using the formula [START_REF]Digital Library of Mathematical Functions[END_REF]8.4

.15] Γ(-n, x) = (-1) n n! (ψ(n + 1) -log(x)) -x -n n-1 k=0 + ∞ k=n+1 (-x) k k!(k -n) . (19)
When x = z 1 ≈ z is sufficiently large, we can also start from z initial = ∞ and use the asymptotic series

Γ(a, x) = x a-1 e -x N -1 n=0 (1 -a) n (-x) n + R N (a, x) (20)
for the first step. We do not need (20) in the proof of Theorem 1, but practically speaking it makes a significant difference for efficiency to choose this expansion whenever min N |R N (a, x)| is smaller than the target tolerance. The error term satisfies

|R N (a, x)| ≤ |(1 -a) N |/|x| N if a ∈ R and x > 0 provided that N ≥ a -1.
For error bounds with complex a, the formulas in [Nat13, §13.7] may be used.

Expansions at generic points.

For consecutive Taylor steps z k-1 z k , we write the expansion as

Γ(a, z k ) = ∞ n=0 c n (z k-1 )x n , x = z k -z k-1 . ( 21
)
where the coefficients c n (z k-1 ) need to be determined. We denote the local expansion point by u instead of z k-1 below to simplify the formulas. The function y(z) = Γ(a, z) satisfies the second-order differential equation zy ′′ + (z -a + 1)y = 0. We can translate this differential equation to the point u and convert it to the second-order linear recurrence relation

u(n + 1)(n + 2)c n+2 (u) + (n + 1)(n + 1 + u -a)c n+1 (u) + nc n (u) = 0. ( 22
)
We can also apply this recurrence to parameter derivatives. Explicitly, define

c j,n (u) = Γ (j,n) (a, u) j!n! (23) so that Γ (j,0) (a, z k ) = ∞ n=0 c j,n (u)
x n , and let S n = J j=0 c j,n (u)X j in the ring of truncated formal power series C[[X]]/ X J+1 . Then

u(n + 1)(n + 2)S n+2 + (n + 1)(n + 1 + u -a -X)S n+1 + nS n = 0. ( 24
)
Truncation bounds for the Taylor series ( 21) and its parameter derivatives can be obtained using the Cauchy integral formula: for n ≥ 1,

c j,n (u) = - d du n-1 u a-1 log j (u)e -u j!n! = - 1 2πi 1 j!n γ t a-1 log j (t)e -t (t -u) n dt. ( 25
)
Lemma 5. For a ∈ C, u > 0, j ≥ 0, n ≥ 1 and 0 < R < u, the coefficients c j,n (u) (and c n (u) = c 0,n (u)) satisfy the bound

|c j,n (u)| ≤ 1 j!n 1 R n-1 M R (u), M R (u) = max t:|t-u|=R |t a-1 log j (t)e -t |. ( 26
)
Consequently, for N ≥ 1 and |x| < R, tails of the Taylor series satisfy

∞ n=N c j,n (u)x n ≤ RM R (u) j!N C N 1 -C , C = |x| R . (27)
To bound M R (u), it suffices to bound t a-1 log j (t)e -t on the disk |t-u| ≤ R using naive upper bounds for the elementary functions, or using interval arithmetic. A simple algorithm to choose R is to start with R = (u + |x|)/2 and iterate R ← R/2 as long as this decreases the bound. The results can be improved slightly with a proper numerical minimization algorithm.

3.2.3. Overall algorithm. The final step is to rewrite the series expansions in matrix form and evaluate the matrix products using binary splitting. We sketch the complete algorithm to compute the incomplete gamma function.

We write x + [±ε] to express the use of an enclosure with midpoint x and radius ε to represent an exact quantity.

Algorithm 6 (Bit-burst evaluation of Γ(a, z), with z > 0).

• Choose initial number of bits b 1 = 32 and let x = z 1 = 2 b1 ⌊z/2 b1 ⌋.

• If the asymptotic series (20) is accurate enough:

-Choose N and denote by ε a bound for the remainder term in (20).

-Compute P = U N -1 • • • U 1 U 0 using binary splitting, where

U n = a-n-1 x 0 1 1 . -Now P 2,1 = N -1 n=0 (1-a)n (-x) n . Compute y 1 = x a-1 e -x (P 2,1 + [±ε]
), which equals Γ(a, z 1 ).

• Otherwise:

-Choose N and denote by ε a bound for the remainder term in (17).

-Compute P = U N -1 • • • U 1 U 0 using binary splitting, where

U n = x a+n+1 0 1 1 . -Now P 2,1 = N -1 n=0 x n
(a+1)n . Compute y 1 = Γ(a)-x a e -x (P 2,1 +[±ε])/a, which equals Γ(a, z 1 ).

-(At a pole of the gamma function, perform the formal limit computation in the above steps.) • For k = 2, 3, . . . with b k = 2b 1 , perform the following:

-

Let z k = 2 b k ⌊z/2 b k ⌋.
If this approximates z to within the target precision, set z k = z instead and make this the last iteration.

-Compute x = z k -z k-1 .
-Choose N and denote by ε a bound for the remainder term in (21).

-Compute P = U N -1 • • • U 1 U 0 using binary splitting, where

U n =    0 x 0 xn Q x(n+1)(n+1+z k-1 -a) Q 0 1 0 1    , Q = -z k-1 (n + 1)(n + 2). -Compute y ′ k-1 = -z a-1 k-1 e -z k-1 . -Compute y k = P 3,1 y k-1 + P 3,2 y ′ k-1 + [±ε], which equals Γ(a, z k ). • Return y k , which equals Γ(a, z).
As noted previously, derivatives up to order j with respect to a can be computed using the same algorithm by substituting a → a+X and working in C[[X]]/ X j+1 .

Let us elaborate on the technical details in the proof of Theorem 1. When computing the values Γ((s + δ)/2, z) and Γ((1 -s + δ)/2, z) used in the approximate functional equation, the recurrence matrices in Algorithm 6 will have entries in the number field Q(s), or in the power series ring R = Q(s)[[X]]/ X j+1 if we compute derivatives. When both j and the degree of s are fixed, this ring is a finite-dimensional vector space over Q, and if the minimal polynomial of s has height h, the product or sum of N entries in R with b-bit coefficients will have coefficients with O(N (b + log h))-bit numerators and denominators.

Summing over all bit sizes in the recursion trees for the binary splitting and the consecutive bit-burst steps and using the bound b 1+o(1) for the bit complexity of arithmetic on b-bit rational numbers, we obtain the p 1+o(1) complexity bound for each call to Algorithm 6.

Implementarion remarks.

In practice, we should clear denominators so that the matrices have integral entries in the binary splitting products. The products should then be computed using with truncation (rounding) to reduce memory usage and improve performance [START_REF] Mezzarobba | A note on the space complexity of fast D-finite function evaluation[END_REF]. Further constant-factor savings are possible by eliminating various redundant computations in the binary splitting process.

The gamma function Γ(a) with algebraic a can be computed in quasilinear time by evaluating Γ(a, N ) with a sufficiently large N using binary splitting [START_REF] Brent | Fast multiple-precision evaluation of elementary functions[END_REF][START_REF] Johansson | Arbitrary-precision computation of the gamma function[END_REF]. In any case, this only needs to be done once: the same Γ(a) value can be recycled for all evaluations of Γ(a, z).

When using (17), there can be significant cancellation between the gamma function and the series. This does not affect the absolute error when s is small, but when s is large, we need to increase the working precision to compensate. The precision increases with z up to the point where we can switch to the asymptotic series (20).

In the pseudocode for Algorithm 6, we evaluate the first derivative y ′ (z) = -z a-1 e -z in each Taylor step. There are several ways to do this: we can compute the elementary functions from scratch, we can perform bit-burst analytic continuation of the function y ′ (z), or we can perform bit-burst evaluation of the factors z a-1 and e -z separately using the standard binomial and exponential Taylor series. Which method performs better may depend on several factors, but either approach achieves quasilinear complexity.

In the approximate functional equation, we need to evaluate Γ(a, z) for successive values z = Cn 2 . It is tempting to reuse the computed Γ(a, z) values, starting the bit-burst evaluation at z initial = C(n -1) 2 . However, this appears to be a net slowdown, the main reason being that the recurrence matrices are much simpler for the hypergeometric series at the origin than for the expansions at generic points.

Implementation results

We have implemented the algorithm for L(s, χ) with s ∈ Q in Arb [START_REF] Johansson | Arb: Efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. We leave an implementation for s ∈ Q and for L (j) (s, χ) for future work. 4.1. Fixed rational points. Table 1 illustrates the precision-dependent scaling of the implementation of the approximate functional equation (AFE) when s is a fixed simple fraction. We also show timings for the Euler-Maclaurin implementation of Dirichlet L-functions in Arb (EM). 3We observe that the AFE is competitive from about 10 4 digits for computing the Riemann zeta function. The advantage is greater for L-functions with larger modulus q due to the O(q 1/2 ) scaling.

At high enough precision, the subquadratic asymptotic complexity of the AFE is evident since the measured time increases by (barely) less than a factor 10 when the precision is multiplied by 10 1/2 . Asymptotically for a p 3/2+o(1) complexity algorithm, the time should only increase by a factor 10 3/4 ≈ 5.6, but the tested precisions are small enough for the hidden logarithmic factors in the complexity bounds to influence the running time. The timings for the p 2+o(1) EM algorithm also increase by factors somewhat larger than 10 for the same reason.

There is roughly a factor two slowdown with both algorithms going from ζ(1/2) to ζ(4/3), for different reasons: the AFE is inherently twice as fast for ζ(1/2) since only one of the two series has to be computed; the slowdown with EM is an implementation artifact (the Arb code is not optimized for rational powers).

The million-digit computation of ζ(1/2), which takes less than five hours on a single core and requires negligible memory, appears to be a precision record for a zeta constant not at an integer. 4 4.2. Computation of Bernoulli numbers. Table 2 compares three algorithms to compute B n as an exact fraction:

• MM: Harvey's implementation of his n 2+o(1) multimodular algorithm [START_REF] Harvey | A multimodular algorithm for computing Bernoulli numbers[END_REF] (available in the bernmm module in SageMath). • EP: the classical n 2+o(1) zeta function algorithm using the Euler product implemented in Arb. • AFE: the approximate functional equation implemented in Arb.

The last two implementations also support computing Euler numbers. The MM and EP algorithms both scale superquadratically with n, the multimodular algorithm having a slight edge for n larger than 10 4 . This confirms the observations in [Har10, Table 1]. The AFE appears to scale weakly subquadratically, but for reasonably sized n, it is roughly 10 3 times slower than the Euler product. To explain this gap we need to consider the logarithmic and constant factor overheads that we have neglected in the complexity analysis so far.

An analysis with Stirling's formula shows that the cutoff in the Euler product for computing ζ(n) to log 2 |B n | bits of accuracy is N ≈ n/(2πe), and there are about N/ log N primes up to this cutoff. A similar analysis for the AFE gives the cutoff N ≈ n 1/2 (2π) -1/2 log(n/(2e)), where all terms are needed, and there are two such series to compute. Considering these facts alone, the AFE thus saves at most a factor n 1/2 /(2e √ 2π log 2 n) ≈ n 1/2 /(13.6 log 2 n) asymptotically, which means we need n ≈ 10 7 to break even assuming that the terms have unit cost.

The last assumption is obviously false: series of matrix products are roughly O(log 2 n) slower than integer powers. An asymptotic speedup of n 1/2 /(C log 4 n), C > 10 1 , is consistent with the observed three-orders-of-magnitude slowdown for n ≈ 10 6 and suggests that we may need n larger than 10 15 for the AFE to win, though there is too much uncertainty to extrapolate reliably.

An interesting question is whether Harvey's subquadratic algorithm for Bernoulli numbers [START_REF] Harvey | A subquadratic algorithm for computing the n-th Bernoulli number[END_REF] can perform better, but unfortunately no implementation exists.

Evaluation for non-algebraic s

If s is not algebraic and instead must be represented by a p-bit floating-point approximation, then the bit-burst algorithm does not offer any improvement over naive series evaluation since the recurrence matrices will not have small entries, and we only obtain a p 5/2+o(1) algorithm to compute L(s, χ) via (2). However, there at least four independent ways to reduce the complexity to p 2+o(1) :

(1) We evaluate each Γ(a, z) via (17) (optionally together with (20)) using a baby-step giant technique, exploiting the hypergeometric structure of the terms: if the truncated series is represented as a matrix product N -1 n=0 U n of length N = m 2 , we expand m-1 n=0 U x+n as a matrix of rational functions in x and evaluate at m points using fast multipoint evaluation [START_REF] Borwein | Reduced complexity evaluation of hypergeometric functions[END_REF]. This achieves p 3/2+o(1) complexity for each incomplete gamma function. This is the approach described in [START_REF] Borwein | Computational strategies for the Riemann zeta function[END_REF].

(2) We expand a truncation of the series (17) (optionally together with (20)) as a polynomial in z of degree p 1+o(1) and evaluate it simultaneously at the requisite p 1/2+o(1) values of z using fast multipoint evaluation. (3) Instead of using the series in incomplete gamma function, we use the integral form (3) and its analog for Dirichlet L-functions. Using a standard numerical integration method with geometric rate of convergence for analytic functions, for example Gaussian, Clenshaw-Curtis or double exponential quadrature, we need p 1+o(1) evaluations of the integrand. Evaluating the theta function in the integrand using the q-series costs p 1/2+o(1) multiplications, resulting in an p 5/2+o(1) algorithm for L(s, χ). This is the method used by Keiper [START_REF] Keiper | Power series expansions of riemann's ξ function[END_REF]. However, we can compute the theta function in quasilinear time using arithmetic-geometric mean iteration instead [START_REF] Labrande | Computing Jacobi's theta in quasi-linear time[END_REF], and this achieves p 2+o(1) complexity for L(s, χ).

(4) As above, but we expand the truncated q-series for the theta function as a polynomial and evaluate it at all the integration nodes using fast multipoint evaluation.

The techniques for fast evaluation of theta functions used in methods (3) and (4) have previously been used in the context of computing class polynomials via numerical approximations of the roots [START_REF] Enge | The complexity of class polynomial computation via floating point approximations[END_REF].

We also mention a version of method (2) that is asymptotically slower but may be superior at realistic levels of precision. We can expand a truncation of the series (17) (and optionally (20)) for Γ(a, πn 2 α/q) as a polynomial in n 2 . This reduces the computation to multipoint evaluation at the small integers n = 1, 2, . . ., which may be performed using repeated applications of Horner's rule instead of fast multipoint evaluation. This results in a p 5/2+o(1) algorithm but with very little overhead since the Horner evaluations only involve additions and p-by-1-word multiplications which are orders of magnitude cheaper than full p-by-p multiplications.

Yet another option is the Booker-Molin method, which employs a Fourier series that can be precomputed for efficient evaluation at many values of s [BC21, §9.4].

It is not clear a priori which of the above methods will perform better, so further implementation studies are needed.

Did Borwein and Borwein [START_REF] Borwein | On the complexity of familiar functions and numbers[END_REF] have one of the methods above in mind for nonalgebraic s, or did they have an entirely different algorithm? It is likely that they considered (1) or (2) since their paper discusses the same multipoint evaluation techniques for other functions, although their wording is more suggestive of an algorithm along the lines of (3) or (4). Either way, there seems to be no obvious way to obtain a subquadratic algorithm for ζ(s) or L(s, χ) for non-algebraic s, and this remains an open problem along with the problem of finding a quasilinear algorithm for s ∈ Q \ Z.

2. 5 .

 5 Values at rational points. Theorem 1 implies p 3/2+o(1) complexity for computing the values L(s, χ) with s ∈ Q. These constants have various applications; ζ(n + 1/2) and β(n + 1/2) appear in thermodynamics (Bose-Einstein statistics) and in connection with lattice sums describing the electrostatic potentials in crystals (Madelung constants) [Fin03, §1.10]. It is a famous open problem whether L(1/2, χ) = 0 for all primitive characters χ [Pla11, §7.6].

Table 1 .

 1 Time in seconds to compute values of L-functions at fixed simple rational points, using Euler-Maclaurin summation (EM) and the approximate functional equation (AFE).

	Digits Number EM	AFE Number	EM AFE
	10 3	ζ(1/2) 0.0076 0.037 L(1/2, χ23.19) 0.15 0.18
	⌊10 3.5 ⌋	0.19	0.29	2.3 1.5
	10 4	2.7	2.7	38 14
	⌊10 4.5 ⌋	52	27	621 131
	10 5	887	262	1282
	⌊10 5.5 ⌋		2175
	10 6		17004
	10 3	ζ(4/3) 0.014 0.083 L(4/3, χ23.19) 1.5 0.38
	⌊10 3.5 ⌋	0.34	0.68	39 3.5
	10 4	6.0	7.2	795 30
	⌊10 4.5 ⌋	100	66	295
	10 5	1808 618	2821
	⌊10 5.5 ⌋		5242

Table 2 .

 2 Time in seconds to compute the Bernoulli number Bn and Euler number En using Harvey's multimodular algorithm (MM), the Euler product (EP), and the approximate functional equation (AFE). Timings marked * were estimated based on the time to evaluate a sparse subsequence (1/10 or 1/100) of the terms, giving an accurate estimate of the time for the full computation without performing it. We indicate the number of digits in the numerator of Bn and in En.

	Number	n	Digits	MM	EP	AFE
	Bn	10 3	1779	0.0066	0.00010	0.067
		⌊10 3.5 ⌋	7180	0.025	0.0011	0.83
		10 4	27691	0.10	0.012	11
		⌊10 4.5 ⌋	103330	0.47	0.18	142
		10 5	376772	2.7	1.9	1707
		⌈10 5.5 ⌉	1349518	22	21	16578
		10 6	4767554	222	224	159945*
		⌈10 6.5 ⌉	16657389	2329	2567	1587800*
	En	10 3	2372		0.00026	0.19
		⌊10 3.5 ⌋	9076		0.0026	2.0
		10 4	33699		0.033	24
		⌊10 4.5 ⌋	122367		0.49	293
		10 5	436962		5.9	2874
		⌈10 5.5 ⌉	1539903		68
		10 6	5369590		726

The initial number of bits is a tuning parameter; instead of the constant

32, we may start with O(log p) bits, for example.

The benchmarks were run on a 1.90 GHz Intel i5-4300U CPU.

In 2013, the author computed the first nontrivial zero 1 2 + 14.134 . . . i of ζ(s) to 303,000 digits using Euler-Maclaurin summation. This took 20 hours and used 62 GB of memory, the high memory usage being the main obstacle to reaching higher precision[START_REF] Johansson | Rigorous high-precision computation of the Hurwitz zeta function and its derivatives[END_REF] (this is an implementation problem that can be avoided).

Acknowledgements

The author was supported in part by the ANR grant ANR-20-CE48-0014-02 NuSCAP.