Pierrick Coupé
email: pierrick.coupe@gmail.com

Simon F Eskildsen

José V Manjón

Vladimir Fonov

Louis Collins

Keywords: hippocampus, hippocampus volume, hippocampus grading, patient's classification, nonlocal means estimator, Alzheimer's disease, entorhinal cortex

come

Simultaneous Segmentation and Grading of Anatomical Structures for Patient's Classification: Application to Alzheimer's Disease

Pierrick Coupé 1 , Simon F. Eskildsen 1 , José V. Manjón 2 , Vladimir S. Fonov 1 , D. Louis Collins 1 and the Alzheimer's Disease Neuroimaging Initiative *

Introduction

The atrophy of medial temporal lobe structures, such as the hippocampus (HC) and entorhinal cortex (EC), is potentially specific and may serve as an early biomarker of Alzheimer's disease (AD) [START_REF] Frisoni | The clinical use of structural MRI in Alzheimer disease[END_REF]. In particular, atrophy of the HC and the EC can be used as a marker of AD progression since changes in these structures are closely related to changes in cognitive performance of the subject [START_REF] Frisoni | The clinical use of structural MRI in Alzheimer disease[END_REF]. The evaluation of structure atrophy is usually estimated by volumetric studies on anatomical MRI, requiring a segmentation step that can be very time consuming when done manually. This limitation can be overcome by using automatic segmentation methods.

In recent years, numerous methods have been proposed to automatically segment the hippocampus [START_REF] Barnes | A comparison of methods for the automated calculation of volumes and atrophy rates in the hippocampus[END_REF][START_REF] Bishop | Novel Fast Marching for Automated Segmentation of the Hippocampus (FMASH): method and validation on clinical data[END_REF]Chupin et al., 2009b;[START_REF] Collins | Towards accurate, automatic segmentation of the hippocampus and amygdala from MRI by augmenting ANIMAL with a template library and label fusion[END_REF]Coupe et al., 2011b;[START_REF] Gousias | Automatic segmentation of brain MRIs of 2-year-olds into 83 regions of interest[END_REF][START_REF] Khan | Optimal weights for local multi-atlas fusion using supervised learning and dynamic information (SuperDyn): validation on hippocampus segmentation[END_REF][START_REF] Lotjonen | Fast and robust multi-atlas segmentation of brain magnetic resonance images[END_REF][START_REF] Morey | A comparison of automated segmentation and manual tracing for quantifying hippocampal and amygdala volumes[END_REF][START_REF] Pohl | Automatic Segmentation Using Non-Rigid Registration[END_REF][START_REF] Van Der Lijn | Hippocampus segmentation in MR images using atlas registration, voxel classification, and graph cuts[END_REF][START_REF] Wang | A learning-based wrapper method to correct systematic errors in automatic image segmentation: consistently improved performance in hippocampus, cortex and brain segmentation[END_REF]. Among these methods, several have been used to classify AD patients using HC volume (Chupin et al., 2009a;[START_REF] Colliot | Discrimination between Alzheimer disease, mild cognitive impairment, and normal aging by using automated segmentation of the hippocampus[END_REF][START_REF] Morra | Comparison of AdaBoost and support vector machines for detecting Alzheimer's disease through automated hippocampal segmentation[END_REF][START_REF] Mueller | Hippocampal atrophy patterns in mild cognitive impairment and Alzheimer's disease[END_REF]. Despite the high segmentation accuracy of the new HC segmentation approaches, using the HC volume enables a separation between AD and cognitively normal (CN) subjects with a success rate only around 72-74% over the entire Alzheimer's Disease Neuroimaging Initiative (ADNI) database [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF]. This limited capability to classify AD patients using the HC volume only, may be due to a simplification of the complex atrophy patterns to a volume -a simple scalar. Recently, several shape analysis methods have been proposed [START_REF] Csernansky | Preclinical detection of Alzheimer's disease: hippocampal shape and volume predict dementia onset in the elderly[END_REF][START_REF] Gerardin | Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging[END_REF][START_REF] Gutman | Disease classification with hippocampal shape invariants[END_REF] to capture detailed HC structural modifications in order to obtain a more accurate classification. At 77% in the comparison proposed by [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF], the approach proposed in [START_REF] Gerardin | Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging[END_REF] yields slightly better classification than a volumetric approach. Therefore, development of new methods capable of estimating subtle anatomical modifications of HC appears to be a critical point to obtain better classification rate. Longitudinal approaches to the AD classification problem have also been investigated by estimating the HC atrophy rate over time [START_REF] Henneman | Hippocampal atrophy rates in Alzheimer disease: added value over whole brain volume measures[END_REF][START_REF] Wolz | Measurement of hippocampal atrophy using 4D graph-cut segmentation: application to ADNI[END_REF]. In [START_REF] Wolz | Measurement of hippocampal atrophy using 4D graph-cut segmentation: application to ADNI[END_REF], the authors reported a correct classification rate of 82% on 568 images of the ADNI dataset. However, this type of approach requires several time-points for a given patient. Finally, an emerging method is to segment subfields of the hippocampus [START_REF] Van Leemput | Automated segmentation of hippocampal subfields from ultra-high resolution in vivo MRI[END_REF][START_REF] Yushkevich | Nearly automatic segmentation of hippocampal subfields in in vivo focal T2-weighted MRI[END_REF]. This approach seems promising since it is potentially able to detect more detailed atrophic patterns. However, ultra-high resolution MRI is required, which is not yet the standard in clinical practice and thus limits the practical applicability of this approach for the moment.

The EC volume has also been investigated as a possible biomarker to detect AD [START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF][START_REF] Du | Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and Alzheimer's disease[END_REF][START_REF] Frisoni | Hippocampal and entorhinal cortex atrophy in frontotemporal dementia and Alzheimer's disease[END_REF][START_REF] Juottonen | Comparative MR analysis of the entorhinal cortex and hippocampus in diagnosing Alzheimer disease[END_REF][START_REF] Xu | Usefulness of MRI measures of entorhinal cortex versus hippocampus in AD[END_REF]. EC atrophy seems to appear slightly earlier in AD progression than HC atrophy, and thus could be used as a more specific biomarker in the initial stages of the disease [START_REF] Frisoni | The clinical use of structural MRI in Alzheimer disease[END_REF]. However, the high intersubject variability of the EC and the difficulty to define EC boundary in anatomical MRI make volumetric studies on EC very challenging [START_REF] Juottonen | Comparative MR analysis of the entorhinal cortex and hippocampus in diagnosing Alzheimer disease[END_REF][START_REF] Kenny | A volumetric magnetic resonance imaging study of entorhinal cortex volume in dementia with lewy bodies. A comparison with Alzheimer's disease and Parkinson's disease with and without dementia[END_REF][START_REF] Xu | Usefulness of MRI measures of entorhinal cortex versus hippocampus in AD[END_REF]. Therefore, studies based on EC volume have been limited to comparison of manual segmentations. Patient classification accuracy using EC volume greatly varies according to the dataset, from 67% [START_REF] Frisoni | Hippocampal and entorhinal cortex atrophy in frontotemporal dementia and Alzheimer's disease[END_REF] up to 87% [START_REF] Juottonen | Comparative MR analysis of the entorhinal cortex and hippocampus in diagnosing Alzheimer disease[END_REF]. Depending on the study, EC volume can be more sensitive than HC volume to separate CN versus AD [START_REF] Juottonen | Comparative MR analysis of the entorhinal cortex and hippocampus in diagnosing Alzheimer disease[END_REF], less sensitive [START_REF] Frisoni | Hippocampal and entorhinal cortex atrophy in frontotemporal dementia and Alzheimer's disease[END_REF] or similarly sensitive [START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF][START_REF] Xu | Usefulness of MRI measures of entorhinal cortex versus hippocampus in AD[END_REF]. As noticed in [START_REF] Xu | Usefulness of MRI measures of entorhinal cortex versus hippocampus in AD[END_REF], it seems that the theoretical advantage of EC measurements over HC measurements is adversely impacted by the difficulty to segment EC due to the ambiguity in defining its boundary in MRI. The development of automatic methods to segment EC is challenging. However, an accurate and consistent EC segmentation method could have an important impact on the use of this structure on large datasets and in a more systematic manner within the study of AD.

In this paper, we propose a new approach designed to address the two problems described earlier: i) to obtain a more detailed detection of structural changes caused by the disease and ii) to perform the automatic segmentation of complex structures such as EC. Inspired by work in image denoising [START_REF] Buades | A non-local algorithm for image denoising[END_REF]Coupe et al., 2008), we have recently proposed a new nonlocal patch-based label fusion method to segment anatomical structures (Coupe et al., 2011b). By taking advantage of pattern redundancy present within the subject's image, as well as the redundancy across training subjects, the nonlocal means scheme enables robust use of a large number of samples during estimation. In Coupé et al. (2011b), we applied this approach to label fusion for the segmentation of anatomical structures such as HC of healthy subjects and lateral ventricles of patients with AD. In this paper, we propose an extension of this patch-based segmentation method in order to evaluate the similarity (in the nonlocal means sense) of the intensity content of one MRI compared to several training populations. By using training populations with different clinical status (e.g., healthy CN subjects and patients with AD), a nonlocal means estimator is used to evaluate the proximity (i.e., the grade of the disease or the degree of anatomical change consistent with disease in the case of AD) of each voxel of the MRI under study compared to the training populations (see Fig. 1). Since the grade estimation and the label fusion steps require the same patch comparison step, simultaneous segmentation and grading of the studied structure can be achieved in one pass without extra computation. In SNIPE, the nonlocal patch-based comparison is used to i) efficiently fuse the labels of MRI in a training database in order to segment EC and HC, and simultaneously to ii) aggregate the clinical status of the populations constituting the training database in order to detect the presence (or not) of the disease. Finally, the average grading value obtained over the segmented structures is proposed as a new biomarker to estimate the clinical status of the subject under study as a potential computerized aid to diagnosis. The contributions of the paper are: i) the introduction of an innovative approach to better characterize the patterns of structural modification caused by the disease (e.g., anatomical changes such as atrophy in case of AD) through the new concept of grading, ii) the presentation of a method to automatically and simultaneously perform the segmentation and the grading of EC and HC, and iii) the demonstration that the proposed approach can be used as a novel biomarker to efficiently achieve patient classification in the context of AD.

Materials and Methods

Method overview

• The nonlocal means estimator:

The nonlocal means filter was first introduced by [START_REF] Buades | A non-local algorithm for image denoising[END_REF] for the purpose of image denoising. In nonlocal means-based approaches [START_REF] Buades | A non-local algorithm for image denoising[END_REF]Coupe et al., 2008), the patch P(x i ) surrounding the voxel x i under study is compared with all the patches P(x j ) of the image Ω (or a subpart of the image) whatever their spatial distance to P(x i ) (i.e., this is the meaning of the term "nonlocal"). According to the patch similarity between P(x i ) and P(x j ), estimated by the sum of squared differences (SSD) measure, each patch receives a weight w(x i , x j ):

(1) where ||.|| 2 is the L2-norm computed between each intensity of the elements of the patches P(x i ) and P(x s,j ), and h 2 is the smoothing parameter of the weighting function. This weighting function is designed to give a weight close to 1 when the SSD is close to zero and a weight close to zero with the SSD is high. Finally, all the intensities u(x j ) of the central voxels of the patches P(x j ) are aggregated through a weighted average using the weights w(x i , x j ). In this way, the denoised intensity û(x i ) of the voxel x i can be efficiently estimated:

(2) Despite its simplicity, the nonlocal means filter has been demonstrated to have excellent denoising performance. This filter is currently one of the most studied denoising filters and many improvements have been proposed since its introduction (see [START_REF] Buades | Image Denoising Methods. A New Nonlocal Principle[END_REF] for a review of these improvements). The efficiency of the nonlocal means filter relies on two intuitive aspects: the pattern redundancy present in an image (i.e., its self-similarity) and the robust detection of samples derived from the same population by using local context (i.e., patch-based comparison):

-First, to improve the accuracy of an estimator, it is possible to reduce the committed error by increasing the number of involved samples. By using an infinite number of samples derived from the same population, the error theoretically converges to zero. To drastically increase the number of samples used, the nonlocal means filter takes advantage of the redundancy of information by using all the similar voxels present over the entire image.

-Second, to ensure that the used samples are derived from the same population, the surrounding neighbor of a voxel can be used to robustly detect similar realizations of the same process. In the nonlocal means approach, this task is achieved by patch-based comparison using SSD. Two voxels with similar surrounding patches are considered as similar and to belong to the same population. More precisely, the nonlocal means filter performs patch comparison to estimate the degree of the similarity between two voxels. This way, each involved sample has a weight (see Eq. 1) reflecting its relevance.

Finally, a simple weighted average (see Eq. 2) is used to aggregate the samples according to their relevance. This way, the resulting estimator embodies the two interesting qualities described above: to build on a large number of samples and to ensure that the involved samples are derived from the same population.

• From denoising to segmentation:

In [START_REF] Coupé | Nonlocal patch-based label fusion for hippocampus segmentation[END_REF]Coupé et al. ( , 2011b)), we were the first to introduce the nonlocal means estimator in the context of segmentation by averaging labels instead of intensities. By using a training library of N subjects, whose segmentations of structures are known, the weighted label fusion is estimated as follows:

(3) where l(x s,j ) is the label (i.e., 0 for background and 1 for structure) given by the expert to the voxel x s,j at location j in training subject s. It has been shown that the nonlocal means estimator v(x i ) provides a robust estimation of the expected label at x i . With a label set of {0,1} voxels with value v(x i )≥0.5 are considered as belonging to the considered structure and the remaining voxels as background.

In [START_REF] Coupé | Nonlocal patch-based label fusion for hippocampus segmentation[END_REF]Coupé et al. ( , 2011b)), we showed that accurate segmentations of anatomical structures can be obtained using this simple patch-based label fusion framework. In addition, to take advantage of the self-similarity of the image as done for denoising, the nonlocal label fusion also relies on intersubject anatomical consistency. Therefore, many similar patches (self-similarity) can be found in every training subject (inter-subject consistency), thus improving the final estimation. Finally, compared to atlas-based methods using nonlinear registration, the nonlocal patch-based approach has the advantage of better handling the inter-subject variability problem. Contrary to the one-to-one correspondence assumed by nonlinear warping methods between the source and the target image, the nonlocal means estimator makes it possible to deal with one-to-many mappings, which better captures the link between subjects' anatomies. This interesting aspect of the nonlocal means estimator has been used to improve video super-resolution without explicit estimation of inter-frame motion [START_REF] Protter | Generalizing the nonlocal-means to super-resolution reconstruction[END_REF][START_REF] Takeda | Super-resolution without explicit subpixel motion estimation[END_REF].

• From segmentation to grading:

In this paper, we propose to extend this segmentation method to efficiently aggregate clinical status (CN or AD) in order to estimate the proximity (in the nonlocal means sense) of each voxel compared to both populations constituting the training library (see Fig. 1). To achieve this goal, we introduce the new concept of patch-based grading that reflects the similarity of the patch surrounding the voxel under study with all the patches present in the different training populations. In this way, the neighborhood information is used to robustly drive the search of anatomical patterns that are specific to a given subset of the training library. When the training populations include data from subsets of subjects in different clinical states, this approach provides an estimation of the grade (i.e., degree of closeness to one group or another) for each voxel:

(4) where p s is the clinical status of the training subject s. In our case, p s =-1 was used for AD status and p s =1 for CN status. A negative grading value (respectively, a positive grading value) g(x i ) indicates that the neighborhood surrounding x i is more characteristic of AD than CN (respectively, of CN than AD) (see Fig. 2). The absolute value |g(x i )| provides the confidence given to the grade estimation. When |g(x i )| is close to zero, the method indicates that the patch under study is similarly present in both populations and thus is not specific to one of the compared populations and provides little discriminatory information. When |g(x i )| is close to 1, the method detects a high proximity of the patch under study with the patches present in one of the training populations and not in the other. Finally, for each subject, an average grading value is computed over all voxels in the estimated structure segmentation (i.e., for all x i with v(x i )≥0.5) for each side (e.g., ! g HC "left or ! g EC "right ). Since the grading and the segmentation involve the same patch comparison step, these structures are extracted at the same time that their grade is estimated (see Fig. 2).

Several strategies can be used to fuse the average grading of the studied structures. First, each side of the structure can be used separately. Second, it is possible to assign the same weight to the left and right HC and EC (e.g., ! g HC = (g HC "left + g HC "right ) /2 ). This strategy of fusing both sides to be more robust to segmentation inaccuracy was used by Chupin et al. in a volumetric study (Chupin et al., 2009a). During our experiments, we found that these two strategies provided similar results for HC and EC. However, for the HC-EC complex, the best strategy was to compute left and right average grading values over HC-EC segmentation (what gives more importance to HC) and then to use the mean of both sides ( ! g HCEC = (g HCEC "left + g HCEC "right ) /2 ). Therefore, we decided to present all the results using the second strategy.

Training library construction

• Datasets

In this study, the ADNI database (www.loni.ucla.edu/ADNI) was used to validate the proposed approach. This database contains both 1.5T and 3.0T T1-w MRI scans. For our experiments, we randomly selected 120 MRI scans, 60 1.5T MRI baseline scans of CN subjects and 60 1.5T MRI baseline scans of patients with AD.

• Preprocessing

All the selected images were preprocessed as follows: 1) correction of inhomogeneities using N3 [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF], 2) registration to the stereotaxic space using a linear transform to the ICBM152 template (1x1x1 mm³ voxel size) [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF] and3) cross-normalization of the MRI intensity using the method proposed in [START_REF] Nyul | Standardizing the MR image intensity scales: making MR intensities have tissue specific meaning[END_REF]. After preprocessing, all the MRIs are coarsely aligned (linear registration), tissue intensities are homogeneous within each MRI volume (inhomogeneity correction) and across the training database (intensity normalization) (see Fig. 1).

• Label propagation

From the 120 processed MRI scans, 20 scans (10 CN and 10 AD) were randomly selected to be used as seed dataset for segmentation. The HC and the EC of this seed dataset were manually segmented by following the protocol defined in [START_REF] Pruessner | Volumetry of temporopolar, perirhinal, entorhinal and parahippocampal cortex from high-resolution MR images: considering the variability of the collateral sulcus[END_REF]. The manual segmentations of the seed dataset were then propagated to the 100 remaining scans constituting our test dataset using the method described in (Coupe et al., 2011b). After the segmentation propagation step, the test dataset was composed of 100 MRI (50 CN subjects and 50 patients with AD) with their corresponding automatic segmentations (see Fig. 1). In our test dataset, the average age of the populations is 74.8 (±4.8) for CN and 74.9 (±6.4) for AD. The age for the two populations is not significantly different (p=0.36, unpaired t-test). In addition, the Mini Mental State Evaluation (MMSE) is 29.1 (±1.2) for CN and 23.2 (±2.0) for AD.

Implementation details

In all experiments described here, the optimal parameters empirically found in Coupé et al. (2011b) for HC segmentation have been used and thus the patch size was fixed to 7x7x7 voxels and the preselection threshold set to th=0.95.

As done in Coupé et al. (2011b), Ω was replaced by a cubic volume V i centered on x i . First, this strategy to use a semi-local paradigm instead of a fully nonlocal paradigm makes the processing computationally practical. In the denoising literature, this approach is used in the majority of the papers and has been shown to produce near-optimal or optimal results except for images with repetitive textures [START_REF] Brox | Efficient nonlocal means for denoising of textural patterns[END_REF]. Second, as shown in Coupé et al. (2011b), in the case of HC segmentation, limitation of the search window provides better results (see left of Fig. 8 in Coupé et al. (2011b)). Since all the images are linearly registered, the patches belonging to HC are located within a restricted area. By using a larger search window, outliers are added that marginally degrade the segmentation and uselessly increases the computational time. While in Coupé et al. (2011b) the search window size was fixed, we used a locally adaptive search window size. The initialization of the search window was set to 9x9x9 voxels as suggested in Coupé et al. (2011b). However, in the case when no similar patches can be found in this search window (i.e., none of the patches pass through the pre-selection), its radius is increased by one voxel until at least one similar patch in each population is found (i.e., at least one patch in each population pass through the pre-selection step). For all the studied subjects, the largest search window size found was 15x15x15 voxels.

The automatic local adaptation of the smoothing parameter h 2 (x i ) (see Eq. 1) proposed in Coupé et al. (2011b) has been slightly modified. During all the experiments, the squared smoothing parameter was set proportional (with λ=0.5) to the minimal SSD:

(5)

The value of lambda slightly changes the segmentation results. When we validated our segmentation method on the ADNI dataset in Coupé et al. (2011a), using λ=0.5 instead of λ=1 changed the median Dice-Kappa values from 0.882 to 0.883 for CN and from 0.836 to 0.838 for AD.

Finally, a subject selection was also applied to reduce the number of training MRI required as suggested in [START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF]. For each structure, the N closest subjects (in terms of SSD over the initialization mask as done in Coupé et al. (2011b)) are equally selected from both populations (N/2 from the CN population and N/2 from the AD population) (see Fig. 1). This is done to ensure that the size of the "patch pool" from the AD population is coarsely similar to the size of the "patch pool" from the CN population.

For a given subject with N=20 (i.e., 10 AD training templates and 10 CN training templates), the segmentation and the grading maps were obtained in less than 4 minutes for left and right HC and less than 2 minutes for left and right EC using a single core of an Intel Core 2 Quad Q6700 processor at 2.66 GHz.

Validation framework

Our validation framework was designed to compare the capability of different SNIPE-based biomarkers to discriminate between patients and controls. The biomarkers studied were: HC volume, HC grade, EC volume and EC grade as well as their combination.

First, to obtain the segmentation and the grade of the subjects within the test dataset, a leave-oneout procedure was performed over the 100 subjects using their corresponding automatic segmentations resulting from the label propagation step (see Fig. 1). For each subject, the N closest training subjects were selected from the 99 remaining subjects in the library The average grading value was then estimated over the EC and the HC segmentations (for both left and right sides) obtained at the same time (see an example in Fig. 2). These segmentations were also used to measure the HC and EC volumes in the stereotaxic space.

Once all the subjects had a volume and a grade for each structure, a quadratic discriminant analysis (QDA) was performed. Each subject was classified by performing a QDA over the 99 remaining subjects. This approach was applied to volume-based classification, grade-based classification and the combination of both for HC, EC and HC + EC. We found that QDA slightly improved the results compared to linear discriminant analysis, especially when the subject's age was used as an additional parameter. The success rate (SR), the specificity (SPE), the sensitivity (SEN), the positive predictive value (PPV) and negative predictive value (NPV) are presented for each of the tested biomarkers (see [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF] for details on these quality metrics).

Results

Figure 2 shows the grading maps obtained for 2 test subjects (1 CN and 1 AD). The corresponding average grading values and the estimated volumes are also provided for left and right HC and for left and right EC. Visually, the CN subject clearly appears closer to the CN population (mainly red color related to values close to 1) while the AD patient is visually closer to the AD population (mainly purple and black colors related to values close to -1). In addition, Fig. 2 also provides a visual assessment of the quality of the segmentation and grading.

Volumetric study

The left column of Fig. 3 shows the volumes for the 100 subjects of the test dataset for HC and EC for N=80 (i.e., 40 CN and 40 AD). The volumetric approach provided a classification success rate of 80% for HC and 69% for EC. The use of both structures at the same time produced a success rate of 78% through our QDA-based classification. This result indicates that the estimated HC volume is more powerful than the EC volume to identify patients with AD. This observation is in accordance with [START_REF] Frisoni | Hippocampal and entorhinal cortex atrophy in frontotemporal dementia and Alzheimer's disease[END_REF]. Our result using only HC volume is slightly superior to a recently published method comparison [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF]. This might come from differences in the test dataset used here or due to a higher accuracy and consistency of the segmentation method used compared to Chupin et al. (2009b). The success rate obtained with EC volume is similar to the results reported in [START_REF] Frisoni | Hippocampal and entorhinal cortex atrophy in frontotemporal dementia and Alzheimer's disease[END_REF] but lower than the values reported in other studies [START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF][START_REF] Xu | Usefulness of MRI measures of entorhinal cortex versus hippocampus in AD[END_REF], all using manual segmentations. Figure 3 shows the higher variability of EC volume compared to HC volume. As mentioned in the introduction, this range of volumes comes from the high inter-subject variability of EC, but may also be due to the difficulty to distinguish EC structure boundaries on anatomical MRI (e.g., identification of the collateral sulcus and the sulcus semiannularis). Due to this last point, less accurate segmentations may be obtained for this structure and thus the introduction of segmentation errors may negatively impact the patient's classification. The use of both structures at the same time did not improve the result compared to the method based on HC only, while improvements have been observed in [START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF] doing similar experiments on manual segmentations.

Grading study

The right column of Fig. 3 shows the grading values for the 100 subjects of the test dataset for HC and EC for N=80. The success rate of the classification was 89% for HC, 78% for EC and 90% for the combination of both structures. For HC, the success rate obtained by using QDA is similar to thresholding the grading value at zero (4 false positives CN and 7 false negatives AD). In fact, in the perfect case, the 50 first subjects (CN) should have positive average grading values and the 50 last (AD) should have negative average grading values. This result indicates that the HC grade estimator is not biased and thus that the sign of the final grading value can be used directly to classify the patient. On the other hand, the EC grade estimator is biased in the sense that the optimal threshold obtained using QDA is superior to zero. As shown on Fig 3, the EC grades of AD are frequently superior to zero, thus indicating a higher similarity with the patches present in CN population. As we will show later, the normal age-related structural changes in the EC may disturb the detection of the disease-related anatomical changes. However, this bias, which depends on the training library used, can be partially compensated for by using QDA, yielding a success rate of 78%. Finally, by computing the average grade value over the HC and the EC improved the HC results and leads to a very high success rate of 90%.

Comparison of SNIPE anatomical biomarkers

In Tab. 1, the SEN, SPE, PPV and NPV obtained by the different SNIPE-based biomarkers considered are presented. These results show that for both structures studied, the classification based on grading provides significantly better results than the volumetric approach (89% vs. 80% for HC and 78% vs. 69% for EC). Moreover, while the combination of HC+EC tends to spoil the results of volumetric analysis, the combination of both slightly improves the results of the grading study. Three different combinations of biomarkers obtained a success rate of 90% during our experiments: HC volume and grade, HC + EC grade, HC + EC volume and grade. In the three cases, the HC grading was used, indicating a potential key role of this new imaging biomarker.

Table 1: Results of the patient classification (AD vs CN) for the different SNIPE-based biomarkers under investigation. These results were obtained by using discriminant analysis through a leaveone-out procedure on the test dataset with N = 80 (i.e.,40 CN and 40 AD)

Impact of the number of selected best training subjects

This experiment presents the impact of the number of selected best training subjects on the studied biomarkers. Figure 4 presents the success rate for all the biomarkers from N=20 (10 CN and 10 AD) to N=80 (40 CN and 40 AD). As some authors have noted that using the age of subjects could increase the classification accuracy (Chupin et al., 2009a;[START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF], we used the subject's age as supplementary information during QDA.

Volume (see top of Fig. 4): For HC, the classification accuracy was quite stable from N=40 to N=80. In (Coupé et al., 2011b), we showed that a plateau in terms of segmentation accuracy was reached around N = 30. For EC, the best results were obtained for N=80. This result seems to indicate that a large library is required to achieve consistent segmentation of EC. Indeed, increasing the size of the "patch pool" and better address issues related to inter-subject variability. The addition of the age as parameter in QDA improved the results of the classification, especially for EC and HC+EC biomarkers. By performing the QDA only with age provided a success rate of 48% in the classification. Finally, at N=60, the HC volume combined with the age provided a success rate of 82%.

Grade (see middle of Fig. 4): For HC, the best classifications were obtained by using high N values (N=60 and N=80). For EC, the best classification rate was obtained for the smallest value of N=20, a result that was not expected. However, by also using age, the best results were obtained for N=80 for EC. For HC and for HC+EC, using the age improved the results of the classification. In these cases, HC-based classification yielded a success rate of 92% and HC+EC a success rate of 93% at N=40 and N=60.

Volume + Grade (see bottom of Fig. 4): By combining the volume and the grade of the SNIPE biomarkers, we obtained slightly better results than by using only the grade, except for EC. By using the age of the subjects, the volume and the grade over HC (with N=60) provided 92% classification accuracy. Combining all the parameters (i.e., volume, grade and age) slightly decreased the results for the biomarkers involving EC compared to use only grade and age.

Relationship between SNIPE grade and age

As shown in the previous experiment, using the subject's age improved the classification based on the grading measure, except for EC with N=20. This supplementary information seems to help distinguish age-related MRI changes from those related to AD pathology. Figure 5 shows the grade values as a function of age on HC + EC with N = 60 (the case with the highest classification accuracy: 93%). It appears that the grading values decrease with age in both populations. This variation indicates that the grading measure captures the age-related anatomical changes (possibly related to atrophy), and thus this observation may explain the better results obtained using age for all the biomarkers except for EC with N=20. As previously mentioned, QDA provides slightly better results than LDA during classification (between 0 to 2% depending on the biomarker studied). This slightly better fitting is assessed by Pearson's coefficient and corresponding p-value of the linear and quadratic regressions presented in Fig. 5. While for CN, the traditional linear model and quadratic provided similar results, it seems that for AD a quadratic model fits better than a linear model. The nonlinear nature of the atrophy related to AD has recently been studied [START_REF] Frisoni | The clinical use of structural MRI in Alzheimer disease[END_REF][START_REF] Frisoni | In vivo mapping of incremental cortical atrophy from incipient to overt Alzheimer's disease[END_REF][START_REF] Jack | Atrophy rates accelerate in amnestic mild cognitive impairment[END_REF]. As noticed in [START_REF] Frisoni | The clinical use of structural MRI in Alzheimer disease[END_REF], while the majority of studies are based on the linear assumption of the AD progression, brain atrophy during AD is not a linear process. In addition, the grade measure is correlated with age while the volume does not appear to be statistically correlated with age since similar regressions provided correlation of r = 0.31 for CN and r = 0.34 for AD with respective p-values of 0.09 and 0.06.

Relationship between SNIPE grade and MMSE score

Finally, the link between the mini mental state examination (MMSE) score and the grade is studied.

The MMSE is a test evaluating the cognitive function of the patient. As noticed in [START_REF] Black | The search for diagnostic and progression markers in AD: so near but still too far?[END_REF], a useful imaging biomarker should have a link with the cognitive decline of the patient with AD usually estimated by using MMSE. Several studies have investigated the relationship between the MMSE score and the volume or the shape of key structures such as HC [START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF][START_REF] Du | Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and Alzheimer's disease[END_REF][START_REF] Thompson | Mapping hippocampal and ventricular change in Alzheimer disease[END_REF] and EC [START_REF] Devanand | Hippocampal and entorhinal atrophy in mild cognitive impairment: prediction of Alzheimer disease[END_REF][START_REF] Du | Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and Alzheimer's disease[END_REF][START_REF] Kenny | A volumetric magnetic resonance imaging study of entorhinal cortex volume in dementia with lewy bodies. A comparison with Alzheimer's disease and Parkinson's disease with and without dementia[END_REF]. As done in [START_REF] Du | Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and Alzheimer's disease[END_REF], we investigated the correlation between MMSE score and anatomical measurements (i.e., volume and grade) for HC and EC. Fig. 6 shows the plots of the grade and the volume as functions of the MMSE score. For both structures, the coefficient of correlation for grade was higher (r = 0.75 for HC and r = 0.58 for EC) than for the volume (r = 0.55 for HC and r = 0.28 for EC). A statistically significant correlation has been found in all cases, consistent with previous literature [START_REF] Du | Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and Alzheimer's disease[END_REF][START_REF] Kenny | A volumetric magnetic resonance imaging study of entorhinal cortex volume in dementia with lewy bodies. A comparison with Alzheimer's disease and Parkinson's disease with and without dementia[END_REF][START_REF] Thompson | Mapping hippocampal and ventricular change in Alzheimer disease[END_REF].

Another trend was that the HC measurements were more consistent with MMSE scores than EC measurements (see Fig. 6). Finally, the HC grade was the biomarker most consistent with MMSE with a high coefficient of correlation (r = 0.75).

In [START_REF] Du | Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and Alzheimer's disease[END_REF], the authors obtained a correlation coefficient of r = 0.48 for HC and r = 0.48 for EC volume based on manual segmentations with a p-value less than 0.001 in both cases. In our experiment, slightly higher correlation was obtained for HC, but a significantly lower value was obtained for EC as assessed by our higher p-value=0.005. However, our correlation coefficient between EC volume and MMSE score is similar to the correlation presented in [START_REF] Kenny | A volumetric magnetic resonance imaging study of entorhinal cortex volume in dementia with lewy bodies. A comparison with Alzheimer's disease and Parkinson's disease with and without dementia[END_REF] (r = 0.34). It should be noted that the estimation of correlation on discrete functions such as MMSE can bias the significance of correlation. However, we wanted to compare our results with previously published studies using this metric.

Discussion

During our experiments we showed that: i) HC-based measures were more discriminant than ECbased measures, ii) the grading had a higher discriminatory capability than the volume, iii) by adding the age, the classification rate improved, especially when using the HC-grade-based metrics, iv) by computing the grade over a larger area (HC+EC) tended to slightly improve results, especially when the subjects' ages were used within the classification model, and v) the optimal size of the number of selected training subjects were N=60 (60% of the full library) in the majority of the situations studied. A balance appears to be required between using a large enough training population and potentially introducing outlier subjects by using all the available subjects. According to the structure of interest, a different number of training subjects could be used. Moreover, by using a larger library, it could be possible to select a higher number of subjects without introduction of outliers. The difficult segmentation of EC due to inter-subject variability could be partially compensated by using non-linear registration of training subjects instead of linear registration. However, this type of approach is more computational intensive. The introduction of shape priors [START_REF] Hu | Appearance-based modeling for segmentation of hippocampus and amygdala using multi-contrast MR imaging[END_REF] could also be a possibility to deal with ambiguity of the EC boundaries.

The SEN, SPE, PPV and NPV obtained by SNIPE are competitive compared to the ten methods compared in [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF] involving voxel-based morphometry (VBM) [START_REF] Ashburner | Voxel-based morphometry--the methods[END_REF], cortical thickness [START_REF] Fischl | Cortical surface-based analysis. II: Inflation, flattening, and a surface-based coordinate system[END_REF], HC volume (Chupin et al., 2009b) and HC shape [START_REF] Gerardin | Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging[END_REF]. In that comparison paper, the best VBM-based approach obtained 89% accuracy; the best method based on cortical thickness obtained 85% accuracy, the best approach using HC volume 74% accuracy and the method using HC shape 77% accuracy. However, during our experiment, only a subset of the entire ADNI database has been used, contrary to the experiments done in [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF]. Moreover, the classification algorithm used in [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF] was a support vector machine while we used a quadratic discriminant analysis approach. Despite these differences, the classification results obtained by using grade only are competitive to the best results reported in [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF]. Moreover, by adding the subjects' age yielded an accuracy of 93%. This result is similar to the highest classification accuracy 93.3% reported on a similar sized subset of ADNI (51 AD and 52 CN) in [START_REF] Zhang | Multimodal classification of Alzheimer's disease and mild cognitive impairment[END_REF]. However, [START_REF] Zhang | Multimodal classification of Alzheimer's disease and mild cognitive impairment[END_REF] used a multimodal approach involving positron emission tomography (PET) and cerebro-spinal fluid (CSF) markers to reach this degree of accuracy. By using only MRI, their method based on volumetric features provided an accuracy of 86.2%.

It appears that using a larger area of analysis by grading several structures tended to improve the grading estimation. The extension of grading to other key structures impacted by AD seems to be an interesting path to follow for further research. Structures such as parahippocampal cortex and perihinal cortex [START_REF] Dickerson | Differential effects of aging and Alzheimer's disease on medial temporal lobe cortical thickness and surface area[END_REF] or fornix and mammillary body [START_REF] Copenhaver | The fornix and mammillary bodies in older adults with Alzheimer's disease, mild cognitive impairment, and cognitive complaints: a volumetric MRI study[END_REF] could be valuable anatomical structures to improve AD detection. Moreover, further work should investigate the spatial distribution of grade maps over the populations. This information could help to detect more discriminant areas for classification and might provide information on the AD progression. Finally, the application of the proposed grading measure to other diseases has a great potential. Moreover, the difficult problem of clinical differentiation (such as AD and frontal lobe dementia for instance) should also be investigated.

Using SSD as similarity metric, our approach is sensitive to inaccuracy in inter-subject intensity normalization. In Coupé et al. (2011a;2011b), we demonstrated that the proposed preprocessing pipeline involving [START_REF] Nyul | Standardizing the MR image intensity scales: making MR intensities have tissue specific meaning[END_REF] provides a sufficiently robust normalization to obtain accurate segmentations. In this paper, we also showed that the preprocessing pipeline used yields high classification accuracy. Nevertheless, any improvements on the inter-subject normalization should yield further improvements in grading estimation. The use of other similarity metrics less sensitive to intensity normalization should be studied in future work. However, according to our experiments, there is no trivial solution since cross-correlation or correlation ratio cannot distinguish constant areas with different means (e.g, in CSF and white matter), mutual information requires a higher number of samples (bigger patch) and introduces the binning problem for histogram construction, and finally the SSIM index also requires matching of intensity. The use of hybrid metrics based on intensity and derivatives could be further investigated.

As for VBM-based approaches, SNIPE requires several scans of each population to be usable. The construction of a large enough training library might be an issue for trials based on a small number of subjects. However, the number of training subjects required by SNIPE is similar to the number required by VBM studies. As noted by [START_REF] Pell | Selection of the control group for VBM analysis: influence of covariates, matching and sample size[END_REF], a group size of 30 to 50 subjects per population is typical in a VBM study while a group size of 70-90 subjects per population is optimal for detection of HC volume loss. In our experiment, we found that 30 subjects from each population is sufficient to provide very high classification rates.

In this proof of concept study of our proposed grading technique, we focused on the problem of AD vs. CN classification. However, the prediction of conversion from prodromal AD (also known as mild cognitive impairment or MCI) to clinically definite AD is more useful from a clinical and diagnostic point of view. The prediction of patients with MCI who will convert to AD and those who will stay stable is an extremely complex task for which no method has yet provided satisfactory classification results [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF][START_REF] Davatzikos | Prediction of MCI to AD conversion, via MRI, CSF biomarkers, and pattern classification[END_REF]. Proposed methods based on structural MRI have been focusing on gray matter loss as markers for prediction. Our proposed grading and segmentation method SNIPE may add valuable information for the problem of prediction.

In this paper, a new method is proposed to robustly detect the patterns of anatomical change in the hippocampus and entorhinal cortex caused by AD. Based on a nonlocal means estimation framework, the proposed novel grading measure (i.e., anatomical change possibly related to atrophy in the context of AD) enables an accurate distinction between CN subjects and patients with AD leading to a classification success rate of 90%. When the subject's age is combined with the grading measure, a success rate of 93% was obtained. These results are competitive compared to the AD detection performance of VBM, cortical thickness, HC volume and HC shape methods extensively compared in [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF]. In contrast to these approaches, SNIPE has the advantage of: i) simplicity (it can be coded in few hundred lines of code), ii) low computational cost (as it does not require non-rigid registration), iii) robustness of the process (all the subjects get final grading maps) and iv) the possibility to achieve individual classifications based on a MRI data from a single time point (contrary to group classifications or longitudinal studies). These first results are promising and indicate that this new structure grading approach could be a useful biomarker to efficiently detect AD. Further work will investigate the possibility to discriminate populations of patients with MCI compared to AD or CN and furthermore, even the possibility of predicting AD.

Figure 1 :

 1 Figure1: Global overview of the proposed grading and segmentation method SNIPE. The N/2 closest subjects in the training library are selected from both populations. All these subjects have automatic segmentations, resulting from the label propagation step, of the structures studied (the left HC segmentations are displayed in this example). For each patch of the subject under study (an AD patient in this example), a comparison is performed with all the patches constituting the selected training subjects (some examples of similar patches are displayed). In this way, the simultaneous segmentation and grading is obtained for the structure studied. The final grading value corresponds to the average value over the estimated segmentation. This procedure is carried out for each structure studied: left and right HC and left and right EC.

  .

	AD vs. CN	SR	SEN	SPE	PPV	NPV
	HC volume	80%	78%	82%	81%	79%
	HC grading	89%	86%	92%	91%	87%
	HC volume and grading	90%	88%	92%	92%	88%
	EC volume	69%	66%	72%	70%	68%
	EC grading	78%	74%	82%	80%	76%
	EC volume and grading	78%	74%	82%	80%	76%
	HC + EC volume	78%	76%	80%	79%	77%
	HC + EC grading	90%	86%	94%	93%	87%
	HC + EC volume and grading	90%	88%	92%	92%	88%

Acknowledgments

This study was supported by the Canadian Institutes of Health Research (CIHR, MOP-84360 and MOP-111169) and CDA (CECR)-Gevas-OE016. This work was also partially supported by the Spanish Health Institute Carlos III through the RETICS Combiomed, RD07/0067/2001. This work benefited from the use of ITK-SNAP from the Insight Segmentation and Registration Toolkit (ITK) for 3D rendering.

CN subject (ID 34)

The success rates are obtained using quadratic discriminant analysis between both populations in a leave-one-out fashion.

Figure 4: Impact of N on the classification accuracy using volume, grade and age as well as their combination. The success rates are obtained using QDA between both populations in a leave-oneout fashion.


Christophe Arene
email: arene@iml.univ-mrs.fr

Romain Cosset
email: romain.cosset@loria.fr

CONSTRUCTION OF A k-COMPLETE ADDITION LAW ON ABELIAN SURFACES WITH RATIONAL THETA CONSTANTS

Keywords: theta functions, Jacobian, genus 2 curve, addition law, completeness, embedding, line bundle, finite field

come

Construction of a k-complete addition law on abelian

surfaces with rational theta constants

Introduction

Cryptographic protocols using abelian varieties, specifically elliptic curves and abelian surfaces, are a promising way of research. They are based on the discrete logarithm problem for which the computation of the addition of two points is central. In particular, one pays attention to two aspects. Obviously, the number of operations needed to compute the equations must be as small as possible. It appears that their domain of definition has also to be taken into account. Indeed, with the development of embedded cryptosystems, the theoretical resistance of the discrete logarithm problem is no longer sufficient to ensure the protocol security, we also have to deal with physical attacks. For instance the implementation of the usual formulae on the Weierstraß model of an elliptic curve is vulnerable against side-channel attacks due to the use of different formulae for a generic addition or a doubling (see [START_REF] Lange | SCA resistant parallel explicit formula for addition and doubling of divisors in the Jacobian of hyperelliptic curves of genus 2[END_REF] for a possible alternative on genus 2 curve cryptosystems).

In this paper we only consider this second problem. Lange and Ruppert [START_REF] Lange | Complete systems of addition laws on abelian varieties[END_REF] first considered complete sets of addition laws, i.e. for all P , Q in A k there is an addition law defined at (P, Q) (see Definition 1.4). An addition law is said to be k-complete if it is defined over A × A(k). Examples of k-complete addition laws in genus 1 included the Edwards curves [START_REF] Edwards | A normal form for elliptic curves[END_REF][START_REF] Bernstein | Faster Addition and Doubling on Elliptic Curves[END_REF] or the twisted Hessian curves [START_REF] Bernstein | Twisted Hessian curves[END_REF][START_REF] Farashahi | Efficient Arithmetic on Hessian Curves[END_REF]. See also [START_REF] Kohel | Addition law structure of elliptic curves[END_REF] for a large study of the structure of the space of addition laws on elliptic curves and the completeness of addition laws acted on by a torsion subgroup. Our aim is to find a k-complete addition law on the Jacobian of genus 2 hyperelliptic curves.

In the first section we introduce the theta coordinates on the Jacobian of an hyperelliptic curve and explain the link with the classical Mumford coordinates. We then sketch the theory of addition laws. In section 2 we explain how to construct in practice an F q -complete addition law.

1.1. Theta functions of level 4, Link with genus 2 curves. In this subsection, we are interested in arithmetic aspects of the Jacobian of a genus 2 curve. We work over C to simplify the introduction and the use of theta functions. But the results remain true over a non binary finite field (see Remark 1.3). For the classical theory of theta functions, the reader is referred to [START_REF] Mumford | Tata lectures on theta. I[END_REF][START_REF] Mumford | Jacobian theta functions and differential equations[END_REF]. Let Ω be an element of the Siegel half-space:

Ω ∈ Mat 2×2 (C), t Ω = Ω, ℑ(Ω) > 0 , the classical Riemann theta function is defined by ϑ (z, Ω) = n∈Z 2 exp iπ t nΩn + 2iπ t n z .
For all elements a et b of Q 2 , the theta function with characteristics a, b is defined by

ϑ [ a b ] (z, Ω) = exp iπ t aΩa + 2iπ t a (z + b) ϑ (z + Ωa + b, Ω) = n∈Z 2 exp iπ t (n + a)Ω(n + a) + 2iπ t (n + a) (z + b) .
The characteristics are considered modulo Z 2 since for all α, β in Z 2 we have

ϑ a+α b+β (z, Ω) = ϑ [ a b ] (z, Ω) exp(2iπ t aβ)
We will consider theta functions of level 4 which means that the characteristics live in 1 2 Z 2 /Z 2 . A classical result of Lefschetz states that the theta functions of level 4 give an embedding of C 2 /ΩZ 2 + Z 2 into P 15 (C). For a proof see [START_REF] Mumford | Abelian varieties[END_REF]p. 29].

For the sake of readability we use the following notations:

Notation 1.1 ([Gau07, Section 7.1]). We index the sixteen theta functions of level 4 as follow:

ϑ 1 (z) = ϑ t (0 0) t (0 0) (z, Ω) , ϑ 2 (z) = ϑ t (0 0) t ( 1 2 1 2 ) (z, Ω) , ϑ 3 (z) = ϑ t (0 0) t ( 1 2 0) (z, Ω) , ϑ 4 (z) = ϑ t (0 0) t (0 1 2 ) (z, Ω) , ϑ 5 (z) = ϑ t ( 1 2 0) t (0 0) (z, Ω) , ϑ 6 (z) = ϑ t ( 1 2 0) t (0 1 2 ) (z, Ω) , ϑ 7 (z) = ϑ t (0 1 2 ) t (0 0) (z, Ω) , ϑ 8 (z) = ϑ t ( 1 2 1 2 ) t (0 0) (z, Ω) , ϑ 9 (z) = ϑ t (0 1 2 ) t ( 1 2 0) (z, Ω) , ϑ 10 (z) = ϑ t ( 1 2 1 2 ) t ( 1 2 1 2 ) (z, Ω) , ϑ 11 (z) = ϑ t (0 1 2 ) t (0 1 2 ) (z, Ω) , ϑ 12 (z) = ϑ t (0 1 2 ) t ( 1 2 1 2 ) (z, Ω) , ϑ 13 (z) = ϑ t ( 1 2 0) t ( 1 2 0) (z, Ω) , ϑ 14 (z) = ϑ t ( 1 2 1 2 ) t ( 1 2 0) (z, Ω) , ϑ 15 (z) = ϑ t ( 1 2 0) t ( 1 2 1 2 ) (z, Ω) , ϑ 16 (z) = ϑ t ( 1 2 1 2 ) t (0 1 2 ) (z, Ω) .
Remark that the first ten theta functions are the even ones and the last six are the odd ones. For simplicity, we drop the Ω. The evaluation at 0 of these functions are called theta constants. We write them ϑ i instead of ϑ i (0).

Consider an hyperelliptic curve C of genus 2. Associated to this curve is its period matrix Ω which is an element of the Siegel half-space. The Abel-Jacobi map is an analytic isomorphism between Jac(C) and C 2 /ΩZ 2 + Z 2 . The Thomae formulae [Tho70] (see also [START_REF] Mumford | Jacobian theta functions and differential equations[END_REF]III.8]) link the 4th power of the theta constants with the parameters of the curve. Up to isomorphisms, we can recover the theta constants by taking well chosen roots [START_REF] Cosset | Computing (ℓ, ℓ)-isogenies in polynomial time on jacobians of genus 2 curves[END_REF]. Assume that the curve is in Rosenhain form:

C : y 2 = f (x) = x(x -1)(x -λ)(x -µ)(x -ν),
then the ordering {0, 1, λ, µ, ν} leads to the following relations:

ϑ 5 ϑ 1 4 = µ λν , ϑ 7 ϑ 1 4 = µ (ν -1) (λ -µ) µ ν (µ -1) (λ -ν) , ϑ 3 ϑ 1 4 = µ (ν -1) (λ -1) λν (µ -1) , ϑ 4 ϑ 1 4 = µ (λ -1) (ν -µ) λ (µ -1) (ν -λ) .
We can take a square root of the preceeding quotients in an arbitrary way.

The other squares of theta constants of level 4 are given by the formulae:

ϑ 2 6 = 1 ν ϑ 2 1 ϑ 2 4 ϑ 2 5 , ϑ 2 8 = 1 λ ϑ 2 1 ϑ 2 7 ϑ 2 5 , ϑ 2 2 = (ν -1) ϑ 2 5 ϑ 2 6 ϑ 2 3 , ϑ 2 9 = (λ -1) ϑ 2 5 ϑ 2 8 ϑ 2 3 , ϑ 2 10 = ϑ 2 1 ϑ 2 2 -ϑ 2 3 ϑ 2 4 ϑ 2 8
, where arbitrary square roots can be taken. Note that we have to take a field extension to take these roots.

We need to have an explicit algebraic morphism between Jac(C) and the image in P 15 (C) of the embedding by the theta functions of level 4. These formulae can be found in [START_REF] Cosset | Computing (ℓ, ℓ)-isogenies in polynomial time on jacobians of genus 2 curves[END_REF] for the genus 2 case and in [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF] for the general case. Let {a 1 , . . . , a 5 } be the ordered roots of f and let

η 1 := t 1 2 , 0; 0, 0 , η 2 := t 1 2 , 0; 1 2 , 0 , η 3 := t 0, 1 2 ; 1 2 , 0 , η 4 := t 0, 1 2 ; 1 2 , 1 2 , η 5 := t 0, 0; 1 2 , 1 2 , η ∞ := t [0, 0; 0, 0] .
For a subset S in {1, . . . , 5, ∞}, we set

η S = i∈S η i ,
and we define η ′ S and η ′′ S to be the first and second part of η S . This notation comes from the fact that if we denote ∞ the point at infinity of C and A i the point with affine coordinate equal to (a i , 0) for i = 1, . . . , 5 and

A ∞ = ∞, then the divisor i∈S (A i ) -#S(∞) is mapped to Ωη ′ S + η ′′ S by the Abel- Jacobi map.
Let • denote the symmetric difference of two sets. All theta functions of level 4 can be written as ϑ [η U •V ] with U := {1, 3, 5} and a subset V of {1, . . . , 5} of odd cardinality. For each such subset, Van Wamelen [vW98] defines the function t V (z) to be t

V (z) = f V ϑ [η U oV ] (z), where f V is a con- stant which is f V = ϑ [0] /ϑ [η U •V ]
for the even functions (i.e. #V = 3) and which is, for the others,

f 1 = -1 √ a 2 -a 1 ϑ 1 ϑ 5 ϑ 6 ϑ 8 ϑ 2 ϑ 3 ϑ 9 ϑ 10 , f 2 = -1 √ a 2 -a 1 ϑ 5 ϑ 6 ϑ 8 ϑ 4 ϑ 7 ϑ 10 , f 3 = -1 √ a 2 -a 1 ϑ 1 ϑ 6 ϑ 2 ϑ 4 , f 4 = 1 √ a 2 -a 1 ϑ 5 ϑ 3 , f 5 = -1 √ a 2 -a 1 ϑ 1 ϑ 8 ϑ 7 ϑ 9 , f ∅ = f {1,2,3,4,5} = -1 √ a 2 -a 1 3 ϑ 2 5 ϑ 2 6 ϑ 2 8 ϑ 2 ϑ 3 ϑ 4 ϑ 7 ϑ 9 ϑ 10 .
The following theorem is a sum up of results from Van Wamelen.

Theorem 1.2. Let D = (P 1 ) + (P 2 ) -2(∞) be a non theta divisor which corresponds to a vector z ∈ C 2 /(ΩZ 2 + Z 2 ). Let (x i , y i ) be the coordinates of the point P i , i = 1, 2. Write (u, v) for the Mumford's polynomials of D. For k ∈ {1, . . . , 5}, and l, m two distinct elements of {1, . . . , 5} \ {k} we have

u(a k ) = t 2 k (z) t 2 ∅ (z) , v(a k ) = Y k,m -Y k,l a l -a m , Y l,m := y 1 (x 2 -a l )(x 2 -a m ) -y 2 (x 1 -a l )(x 1 -a m ) x 2 -x 1 = c 1,2 t l (z)t m (z)t {l,m} (z) t 3 ∅ (z) , Y := y 1 y 2 = 5 l=1 t l (z) t ∅ (z) ,
where c 1,2 is just a sign ±1.

By evaluating u at the roots of f , we obtain formulae for computing all the ϑ i (z) 2 /ϑ 16 (z) 2 with 1 ≤ i ≤ 16. To get the theta functions of level 4, we will use the doubling formulae [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF]:

4ϑ [ a b ] (2z) ϑ [ a b ] ϑ [ 0 0 ] 2 = α,β∈ 1 2 Z 2 /Z 2 exp -4iπ t aβ ϑ a+α b+β (z) 2 ϑ [ α β ] (z) 2 , 4ϑ [ a b ] (2z) ϑ [ a 0 ] ϑ 0 b ϑ [ 0 0 ] = α,β∈ 1 2 Z 2 /Z 2 exp -4iπ t aβ ϑ a+α b+β (z) ϑ a+α β (z) ϑ [ α b+β ] (z) ϑ [ α β ] (z) .
The first formula allows to recover the even theta functions. For the odd theta functions, we will use the second formula. The products on the right side can be expressed in terms of the constants f V and the functions Y l,m , Y and u(a i ). Since we need to divide by some u(a i ), we make the hypothesis that the divisor is not of 2-torsion. For instance, the second formula gives

ϑ 16 (2z)ϑ 1 ϑ 4 ϑ 8 = ϑ 1 (z)ϑ 4 (z)ϑ 8 (z)ϑ 16 (z) -ϑ 9 (z)ϑ 12 (z)ϑ 13 (z)ϑ 15 (z) + ϑ 5 (z)ϑ 6 (z)ϑ 7 (z)ϑ 11 (z) -ϑ 2 (z)ϑ 3 (z)ϑ 10 (z)ϑ 14 (z), ϑ 16 (2z)ϑ 1 ϑ 4 ϑ 8 = t 2,4 (z)t 2,3 (z)t 3,4 (z)t ∅ (z) f 2,4 f 2,3 f 3,4 f ∅ + t 1,5 (z)t 2 (z)t 4 (z)t 3 (z) f 1,5 f 2 f 4 f 3 + t 3,5 (z)t 4,5 (z)t 2,5 (z)t 1 (z) f 3,5 f 4,5 f 2,5 f 1 + t 1,3 (z)t 1,4 (z)t 1,2 (z)t 5 (z) f 1,3 f 1,4 f 1,2 f 5 , ϑ 16 (2z)ϑ 1 ϑ 4 ϑ 8 t 4 ∅ (z) = Y 2,4 Y 2,3 Y 3,4 u(a 2 )u(a 3 )u(a 4 ) 1 f 2,4 f 2,3 f 3,4 f ∅ + Y 1,5 Y u(a 1 )u(a 5 ) 1 f 1,5 f 2 f 3 f 4 + Y 2,5 Y 3,5 Y 4,5 Y u(a 2 )u(a 3 )u(a 4 )u(a 5 ) 2 1 f 2,5 f 3,5 f 4,5 f 1 + Y 1,2 Y 1,3 Y 1,4 Y u(a 1 ) 2 u(a 2 )u(a 3 )u(a 4 ) 1 f 1,2 f 1,3 f 1,4 f 5 .
Remark 1.3. Although we have defined our theta function over C, our results apply to other fields (of characteristic different from 2). To prove this over a finite field (the relevant case in cryptography), we can use Lefschetz's principle and reduction to prove all the results for ordinary varieties. In general we can always use the algebraic theory of theta functions [Mum66, Mum67a

1.2. Addition laws. Now, we focus on the notion of addition law. Let k be a field and A/k be an abelian variety of dimension g. We assume an embedding of A in some projective space P r is fixed and given by a very ample line bundle L = L(D) for D an effective divisor. We denote by ι : A ֒→ P r the corresponding morphism and also assume in the sequel that this embedding is projectively normal. Recall that by definition this means that for every n ≥ 1 the restriction map

H 0 P r , O P r (n) → H 0 (A, L n ) is
surjective. This is the case in the classical settings where L = L n 0 0 with L 0 an ample line bundle and n 0 ≥ 3 [BL04, p.187].

Let I 1 (resp. I 2 ) be the homogeneous ideal in k[X 0 , . . . , X r ] (resp. in k[Y 0 , . . . , Y r ]) defined by A. The group law

µ : A × A → A, (X, Y ) → X + Y
can be locally described by bihomogenous polynomials. More precisely, an addition law p of bidegree (m, n) on ι(A) ⊂ P r is the data of r+1 polynomials

p 0 , . . . , p r ∈ k[X 0 , . . . , X r ]/I 1 ⊗ k[Y 0 , . . . , Y r ]/I 2 ,
not all zero, bihomogeneous of degree m in X 0 , . . . , X r and degree n in Y 0 , . . . , Y r such that we have

ι • µ (X, Y ) = p 0 ι(X), ι(Y ) : . . . : p r ι(X), ι(Y )
for all points (X, Y ) ∈ A × A k where these polynomials are not all zero. The set of points where an addition law is not defined is called its exceptional subset. It will be convenient for our purpose in Section 2 to use the structure of k-vector space of addition laws having fixed bidegree. In this sense we need to define the zero addition law (independent of the bidegree) given by zero polynomials. It is denoted by 0 and its exceptional subset is A × A k .

In this paper we are interested in the construction of a single addition law which describes the group morphism µ on A × A (k) where k is a non binary finite field and A/k an abelian surface (i.e. the Jacobian of a genus 2 curve) embedded in P 15 . Definition 1.4. A set S of addition laws is said to be k-complete if for any point (X, Y ) ∈ A × A (k) there is an addition law in S defined on an open subset containing (X, Y ). The set S is said to be complete if the previous property is true over k. If S = {p} is a singleton, we say the addition law p is k-complete (or complete when k = k).

Given m, n ≥ 2 the following proposition interprets the addition laws of bidegree (m, n) as global sections of a certain line bundle M m,n . A more explicit description of this link can be found in [AKR11, Section 2].

Proposition 1.5 ([LR85, Lemma 2.1]). Let π 1 , π 2 : A × A → A be the projection maps on the first and second factor. There is an addition law (respectively a complete set of addition laws) of bidegree (m, n) on A with respect to the embedding in P r determined by L if and only if

H 0 (A × A, M m,n ) = 0
(respectively the linear system | M m,n | is base point-free), where

M m,n = µ * L -1 ⊗ π * 1 L m ⊗ π * 2 L n .
The next lemma is specific to the biquadratic case (m=n=2) and gives a nice description of the line bundle M 2,2 which have no equivalent statement for others bidegrees that we currently know.

Lemma 1.6 ([LR85, Propositions 2.2 and 2.3]). Let L be an ample line bundle on A and δ :

A × A → A be the difference map (X, Y ) → X -Y . 1) if L is not symmetric then H 0 (A × A, M 2,2 ) = 0. 2) if L is symmetric, then M 2,2 = δ * L. Moreover M 2,2 is base point-free and h 0 (A × A, M 2,2 ) = h 0 (A, L).
We end this section with the statement of the existence of the addition law we want to construct.

Proposition 1.7 ([AKR11, Statement and Proof of Theorem 4.8]). Let k = F q , q ≥ 7, be a finite field and C/F q be a genus 2 curve. There exists an F q -complete biquadratic addition law on the embedding of Jac(C) in P 15 by 4Θ. Moreover its exceptional subset is explicitely determined.

Construction

2.1.

A basis of biquadratic addition laws on Jac(C) ֒→ P 15 . Riemann's addition formulae are widely known and common in the litterature. We use the general formulae given by Baily [START_REF] Baily | On the Theory of θ-Functions, the Moduli of Abelian Varieties, and the Moduli of Curves[END_REF] and apply it to obtain the following formulae for theta function of level 4.

Proposition 2.1 ([Bai62, Section 2.2, Formulae (9)]). Let a k , b l ∈ 1 2 Z 2 /Z 2 , k, l = 1, . . . , 4. Assume we have -a 1 + a 2 + a 3 + a 4 = 2a, -b 1 + b 2 + b 3 + b 4 = 2b
with a and b in 1 2 Z 2 /Z 2 then for all z 1 , z 2 in C 2 we have 4 ϑ

a 1 b 1 (z 1 + z 2 ) ϑ a 2 b 2 (z 1 -z 2 ) ϑ a 3 b 3 (0) ϑ a 4 b 4 (0) = α,β∈ 1 2 Z 2 /Z 2 ϑ a 1 +a+α b 1 +b+β (z 2 ) ϑ a 2 -a+α b 2 -b+β (z 2 ) ϑ a 3 -a+α b 3 -b+β (z 1 ) ϑ a 4 -a+α b 4 -b+β (z 1 ). For all a 1 , a 2 , b 1 , b 2 in 1 2 Z 2 /Z 2 , there exists a 3 , a 4 , b 3 , b 4 in 1 2 Z 2 /Z 2
verifying the condition of the proposition and such that the constant ϑ

a 3 b 3 ϑ a 4
b 4 is non zero. We now go back to notation 1.1.

Remark 2.2. The embedding Jac(C) ֒→ P 15 is given by the line bundle L = L(4Θ). From now on we cosider the functions ϑ i as global sections of this line bundle.

Remark 2.3. The formulae above express for i, j = 1, . . . , 16 the product ϑ i (z 1 + z 2 )ϑ j (z 1 -z 2 ) as a biquadratic bihomogeneous polynomial in the level 4 theta functions ϑ 1 (z 1 ), . . . , ϑ 16 (z 1 ) and ϑ 1 (z 2 ), . . . , ϑ 16 (z 2 ) . Note also that they are defined over the field of definition of the theta constants.

Next, fixing the index j, if z 1 , z 2 are such that ϑ j (z 1 -z 2 ) = 0, there exists biquadratic bihomogeneous polynomials p i,j such that

ϑ i (z 1 + z 2 )ϑ j (z 1 -z 2 ) = p i,j (ϑ 1 (z 1 ), . . . , ϑ 16 (z 1 )), (ϑ 1 (z 2 ), . . . , ϑ 16 (z 2 )) .
This allows us to construct an addition law p j = (p 1,j , . . . , p 16,j ) defined outside the exceptional subset δ * (ϑ j ) 0 . Indeed, let

X k = ϑ 1 (z k ) : • • • : ϑ 16 (z k ) ∈ Jac(C), k = 1, 2, be two points such that X 1 -X 2 ∈ (ϑ j ) 0 (or satisfying ϑ j (z 1 -z 2 ) = 0). We have ι • µ (X 1 , X 2 ) = ϑ 1 (z 1 + z 2 ) : • • • : ϑ 16 (z 1 + z 2 ) = ϑ j (z 1 -z 2 )ϑ 1 (z 1 + z 2 ) : • • • : ϑ j (z 1 -z 2 )ϑ 16 (z 1 + z 2 ) = p j (X 1 , X 2 ).
Notation 2.4. For j = 1, . . . , 16, we denote p j the addition law on Jac(C) whose exceptional subset is δ * (ϑ j ) 0 presented above.

Clearly the set of addition laws {p 1 , . . . , p 16 } is complete. Moreover we have the following proposition: Proposition 2.5. Let C be a genus 2 curve. The set {p 1 , . . . , p 16 } is a basis of the set of biquadratic addition laws on Jac(C) ֒→ P 15 .

Proof. We have dim k L(4Θ) = 16, so by Lemma 1.6 case 2) we only need to show that the family is free. Let assume there exists a linear relation

λ j p j = 0. (1)
Let denote by O J the neutral element of Jac(C), then for all X = ϑ 1 (z) : • • • : ϑ 16 (z) ∈ Jac(C), the relation λ j p j (X, O J ) = 0 gives λ j p i,j (X, O J ) = 0 for all i = 1, . . . , 16. Moreover there exists a k 0 such that ϑ k 0 (z) = 0, so

0 = λ j p k 0 ,j (X, O J ) = λ j ϑ k 0 (z + 0)ϑ j (z -0) = ϑ k 0 (z) λ j ϑ j (z).
The dependance in k 0 being eliminated we finally get

λ j ϑ j = 0
which is wrong because the family {ϑ j , j = 1, . . . , 16} is a basis for the theta functions of level 4. Hence the assumption of the existence of the relation (1) is not true, and {p 1 , . . . , p 16 } is a free family.

2.2. Idea of the construction. From now on and without mention of the contrary we assume that k = F q , where q is greater than 7 and is odd. In the previous subsection we built a basis for the space of addition laws we are interested in. We want here to construct the addition law announced in Proposition 1.7. We denote it p. According to Lemma 1.6 case 2) its exceptional subset is of the form δ * D with D ∈ Div(Jac(C)) and D ∼ 4Θ.

We recall below the expression of D but suppose here we know it. Our aim is to find a projective solution (λ 1 : • • • : λ 16 ) for the relation p = λ j p j (2) using an interpolation method. To get this, let X ∈ D. As we want p not to be defined on δ * D, in particular at (X, O J ), we search for solutions of the linear system 0 = λ j p j (X, O J ).

Varying X ∈ D we expect to get a linear system of rank 15. Note that the solution is projective because the exceptional subset of an addition law defines it up to scalar multiplications. We recall here the construction of the divisor D ∈ Div k Jac(C) introduced above. The assumption q ≥ 7 implies the existence of a degree 4 closed point of the form P 0 , P σ 0 , P 0 , P 0 σ , with σ the Frobenius over F q . Let define α 0 := (P 0 ) + (P σ 0 ) -2(∞) and for l = 0, 1, 2, α l+1 := α σ l . Then the divisor D = Θ α l has the desired properties to induce an F q -complete biquadratic addition law (namely it is k-rational without k-rational points and linearly equivalent to 4Θ [AKR11]), where Θ α l , l = 0, . . . , 3, is the translation by α l of the theta divisor Θ on Jac(C). In the sequel, D is taken of this form.

The following proposition allows to avoid the computation of the last six coefficients (the odd ones) and then to reduce significantly the running time. Recall that the addition laws p j have δ * (ϑ j ) 0 as exceptional subset (Notation 2.4).

Theorem 2.6. Assume k = F q , with q ≥ 7 and (2, q) = 1. Let p be the addition law introduced above and p = λ j p j the desired linear relation. We have

λ 11 = • • • = λ 16 = 0.
Proof. Remark that by construction for all X ∈ D we have -X ∈ D. We have

∀X ∈ D, p(X, O J ) = p(O J , X) = 0. (3)
Using the parity of the theta functions ϑ j we get that the second equality becomes

p(O J , X) = λ j p j (O J , X) = 10 j=1 λ j p j (X, O J ) - 16 j=11 λ j p j (X, O J ).
We use it in the formulae (3) and are led to consider the two next equations ∀X ∈ D, Remark 2.7. We did not get more information on the coefficients λ j using that p(-X, O J ) = p(O J , -X) = 0 for X ∈ D.

Numerical results.

AVIsogenies is a Magma package for working with genus 2 curves (and more generally with abelian varieties) using theta functions1 . Using some already implemented functions, we wrote codes to compute the coefficients λ i given an hyperelliptic curves. This code is now part of the AVIsogenies package. This compution took less than a minute. It is possible to check that the addition law is F 10007 -complete by an exhaustive computation. Note that it is enough to check p(D, O J ) = D for all divisor ±D of Jac(C) (F 10007 ). This verification took almost a week.

Concerning the efficiency of these addition laws, it is clearly not to their advantage when we look at p 1 expressed below in Example A.1. One verifies in this appendix that the total cost to compute the desired addition law p is 736m + 32s + 124m ϑ , where m denotes a multiplication, s is for a squaring and m ϑ represents a multiplication by a coefficient that only depends on the theta constants, which can be precomputed.

In comparison, the classical representation of points in Jac(C) as elements of the divisor class group of C and the use of Mumford's representation and Cantor's algorithm provides extremly cheaper costs, e.g. 47m + 4s for a general addition in even characteristic (see [START_REF] Lange | Formulae for arithmetic on genus 2 hyperelliptic curves[END_REF]). There also exist pseudoaddition laws on the Kummer surface of the variety that can be computed much faster [START_REF] Duquesne | Montgomery scalar multiplication for genus 2 curves[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF].

10 j=1λ

 10 j p j (X, O J ) = 0, and 16 j=11 λ j p j (X, O J ) = 0. (4) Let us define the two biquadratic addition laws appearing here p 1 := 10 j=1 λ j p j , p 2 := 16 j=11 λ j p j . Let δ * D 1 , δ * D 2 , with D 1 , D 2 be two divisors on Jac(C), be their respective exceptional subsets. We want to prove that p 2 is zero. They verify for k = 1, 2 either D k ∼ 4Θ or p k = 0. The formulae (4) imply D ≤ D k , hence either D = D k and then p k = λp for some λ ∈ F q or p k = 0 for k = 1, 2. But p 2 (O J , O J ) = 0 because the theta constants involved are zero, moreover the k-rational point (O J , O J ) is not an element of δ * D, hence the second addition law p 2 is zero. A fortiori λ 11 = • • • = λ 16 = 0 and p = p 1 .

Example 2. 8 .

 8 Consider the curveC : y 2 = f (x) = x 5 + 5782x 4 + 2517x 3 + 2312x 2 + 9402xdefined over F 10007 . The non-zero associated theta constants areϑ 1 = 1, ϑ 2 = 5242, ϑ 3 = 7727, ϑ 4 = 678, ϑ 5 = 3926, ϑ 6 = 7092, ϑ 7 = 5628, ϑ 8 = 7556, ϑ 9 = 3666, ϑ 10 = 904. Let K = F 10007 [X]/ X 2 + 1 ≃ F 10007 2and x 0 = 8310 + 2164 √ -1. The point P 0 = x 0 , f (x 0 ) is a point of the curve C (F 10007 4 ) which doesn't belong to C (F 10007 2 ). The corresponding non-zero λ i are given byλ 1 = 1, λ 2 =6924, λ 3 = 1940, λ 4 = 9380, λ 5 = 5155, λ 6 = 1278, λ 7 = 7239, λ 8 = 1761, λ 9 = 6859, λ 10 = 5891.

The authors acknowledge the financial support by grant ANR-09-BLAN-0020-01 from the French ANR and the AXA Research Fund for the PhD grant of the first author. 1

It can be found at http://avisogenies.gforge.inria.fr/.

m ϑ = 45m ϑ . For the remaining p i,j with 11 ≤ i ≤ 16 and 1 ≤ j ≤ 10, we point out that if a monomial X i 0 X j 0 Y k 0 Y l 0 appears, so does X k 0 X l 0 Y i 0 Y j 0 with the same sign. We use then the relation

Appendix A. Operation count

We start by computing the addition laws p i , i = 1, . . . , 10 and then use the Formula (2). We remark that there are eight bihomogeneous monomials appearing in p i,j , i = j. Also, p i,j and p j,i are defined by the same monomials up to a sign; this is the case for the p i,i , i = 1, . . . , 10, too. Now we describe the cost of their computation. We do not take into account additions or sign changes costs. Given two points X 1 : • • • : X 16 and Y 1 : • • • : Y 16 we first compute all the products X i X j and Y i Y j , this costs 240m + 32s and the products X i X j Y i Y j in 256m. These monomials are exactly the one included in the ten first polynomials of the addition laws p i (see Example A.1), so the polynomials p i,i are calculated in 10m ϑ and the p i,j , 1 ≤ i, j ≤ 10, in

to calculate each p i,j with 4m+1m ϑ . Hence, the ten addition laws p 1 , . . . , p 10 can be computed in 736m + 32s + 115m ϑ . Finally the computation of the k-complete addition law p requires the 9 multiplications by the coefficients λ i which also can be precomputed, so we count them as 9m ϑ .

Example A.1. As an illustrative exemple, we present the addition law p 1 .