Hugues De Lassus Saint-Geniès

David Defour
email: david.defour@univ-perp.fr

Guillaume Revy
email: guillaume.revy@univ-perp.fr

Range Reduction Based on Pythagorean Triples for Trigonometric Function Evaluation

Software evaluation of elementary functions usually requires three steps: a range reduction, a polynomial evaluation, and a reconstruction step. These evaluation schemes are designed to give the best performance for a given accuracy, which requires a fine control of errors. One of the main issues is to minimize the number of sources of error and/or their influence on the final result. The work presented in this article addresses this problem as it removes one source of error for the evaluation of trigonometric functions. We propose a method that eliminates rounding errors from tabulated values used in the second range reduction for the sine and cosine evaluation. When targeting correct rounding, we show that such tables are smaller, and make the reconstruction step cheaper than existing method. This approach relies on Pythagorean triples generators. Finally, we show how to generate tables indexed by up to 10 bits in a reasonable time and with little memory consumption.

I. INTRODUCTION

The representation formats (binary16, binary32, and binary64) and the behavior of floating-point arithmetics available in computers are defined by the IEEE 754-2008 standard [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]. For the five basic operations (+, -, ×, /, and √ ), this standard requires the system to return the rounding of the exact result, according to one of the four rounding modes (to nearest, toward -∞, toward +∞, and toward 0). This property is called correct rounding, and it warranties the quality of the result. However, due to the Table Maker's Dilemma [START_REF] Lefèvre | The table maker's dilemma[END_REF] and the difficulties to develop accurate and efficient evaluation schemes, correct rounding is only recommended for elementary functions in the IEEE 754-2008 standard.

The algorithms used for the evaluation of elementary functions, like sine, cosine, logarithm, and exponential, can be classified into at least two categories. The first category concerns algorithms based on small and custom operators combined with tabulated values that target small accuracy, less than 30 bits [START_REF] Dassarma | Faithful bipartite ROM reciprocal tables[END_REF], [START_REF] De Dinechin | Multipartite table methods[END_REF], [START_REF] Defour | A new scheme for tablebased evaluation of functions[END_REF]. The second category concerns algorithms that usually target single or double precision (24 or 53 bits) with an implementation on general processors that rely on the available hardware units [START_REF]CR-Libm, a library of correctly rounded elementary functions in double-precision[END_REF], [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF].

Implementations of those functions that target correct rounding [START_REF] Muller | Elementary functions -algorithms and implementation[END_REF] are usually divided into two or more phases. A quick phase based on a fast approximation which ensures a few extra bits than the targeted format. With this phase, rounding is possible most of the time at a reasonable cost. When rounding is not possible, a much slower but more accurate phase is used. The quick phase uses operations with a precision similar to the input and output precision, while the accurate phase is based on extended precision. For example, in the correctly rounded library CR-Libm, the quick phase for the sine and cosine function in double precision targets 66 bits while the accurate phase corresponds to 200 bits [6, § 7]. The accuracy used for the accurate phase is linked with the search for worst cases for the Table Makers Dilemma [START_REF] Lefevre | Hardest-to-round cases[END_REF]. This corresponds to 142 bits for trigonometric functions in rounding to nearest in double precision over the range [2 -25 , (1 + 25317/2 16 )2 9 ].

Given a floating-point number x, the internal design of each phase is based on the following 4-step process:

1) A first range reduction, based on mathematical properties, narrows the domain of the function to a smaller one.

For the sine function, this consists in using the following property:

sin(x + k • π/2) = f k (x * ) with f k (x * ) = ± sin(x * ) or f k (x * ) = ± cos(x * )
depending on k mod 4. This step results in a reduced argument x * in the range [0, π 4 ]. 2) A second range reduction, based on tabulated values, reduces further the range [0, π 4 ]. The argument x * is split into two parts, x * h and x * l , such that:

x * = x * h + x * l with |x * l | < 2 -p . (1)
The term x * h corresponds to the p leading bits of x * , which are used to address a table of π/4 × 2 p entries with precomputed values of sin h = sin(x * h ) and cos h = cos(x * h ). In the case of sine, we have

sin(x * ) = sin h • cos(x * l ) + cos h • sin(x * l ).
3) Then, a polynomial evaluation is used to approximate the remaining terms over the restricted range. In the case of sine function, this corresponds to the evaluation of two polynomials of small degree that approximate sin(x * l ) and cos(x * l ) over the range [0, 2 -p ]. 4) Finally, all terms computed during steps 2 and 3 are merged together during the reconstruction step.

Satisfactory solutions exist to address the first range reduction, the generation and evaluation of precise and efficient polynomial evaluation schemes, and the reconstruction step. The interested reader can find more details in [START_REF] Muller | Elementary functions -algorithms and implementation[END_REF].

In this article, we address the second range reduction based on tabulated values and more specifically for the sine and cosine functions. The proposed method relies on precomputed sine and cosine values with remarkable properties that simplify and accelerate the evaluation of these functions. These properties are threefold: (1) First, each value has to represent exactly both the sine and cosine of a given number, that is, without any rounding error. For this purpose, we use points that are rational numbers. (2) Second, each numerator of rational values should be exactly representable in a representation format available in hardware, that is, as an integer or a floating-point number. (3) Third, each of them has to share the same denominator. These three properties put all together lead to tabulated values for sine and cosine that are exact and exactly representable with rational values sharing a common denominator. Rational numbers that fulfill these three properties are linked to Pythagorean triples.

This article is organized as follows: Section II gives some background on the second range reduction step for the sine and cosine functions. Section III details properties of the proposed tables and shows how they remove two sources of error involved in this step. Section IV presents our approach to build tables of exact points for sine and cosine using Pythagorean triples. Then, Section V presents some experimental results which show that we can compute tables up to 10 bits of index in a reasonable execution time and memory consumption. Finally, some comparison results with other classical approaches are given in Section VI, before concluding in Section VII.

II. PRELIMINARIES

The second range reduction is usually implemented using table lookup methods. This section presents three solutions that address this step. Note that in the sequel of this article, without loss of generality, we will consider the sine function evaluation, while cosine evaluation can be easily derived.

Here and hereafter we denote by x the rounded value of x to a given precision and according to a given rounding mode.

To characterize rounding error we use -p that corresponds to a number y such that |y| ≤ 2 -p . For example, in the case of rounding to the nearest in double precision, we have

x ≤ x(1 + -54 ).

A. Tang's Tables

Tang proposed a general method used to implement elementary functions that relies on hardware-tabulated values [START_REF] Tang | Table-lookup algorithms for elementary functions and their error analysis[END_REF]. Given the reduced argument x * as in (1), Tang's method uses the upper part x * h to address a table. It retrieves two values sin h and cos h which are good approximations of sin(x * h ) and cos(x * h ), respectively, rounded to the destination format. Hence, we have:

sin h = sin(x * h ) • (1 + -54n ) and cos h = cos(x * h ) • (1 + -54n ),
where n is the number of floating-point numbers used to represent those quantities.

Then, these values are combined with the results of the evaluation of the two polynomials S(x) and C(x), defined as:

S(x) = sin(x) -x and C(x) = cos(x) -1.
Finally, the value sin(x * ) is reconstructed as follows:

sin(x * ) = sin h • cos(x * l ) + cos h • sin(x * l ) = sin h • 1 + C(x * l ) + cos h • x * l + S(x * l )
. In order to reach a given accuracy, this method requires the use of a higher precision than the targeted one. This extra precision usually involves the use of error-free algorithms such as Fast2Sum or exact multiplication [START_REF] Knuth | The Art of Computer Programming[END_REF].

Tang's method is well suited for hardware implementations on modern architectures. It takes advantage of the capability on these architectures to access tabulated values in memory and to perform floating-point computations concurrently. Once the argument x * is split into two parts x * h and x * l , the memory units can provide the two tabulated values sin h and cos h , while floating-point units (FPUs) evaluate the polynomials C(x * l ) and S(x * l ). As the degree of the polynomials decreases when the table size increases, the objective is to find parameters such that the polynomial evaluations take as long as memory accesses, in average [START_REF] Defour | Cache-optimised methods for the evaluation of elementary functions[END_REF].

B. Gal's Accurate Tables

In Tang's method, sin h and cos h are approximations of sin(x * h ) and cos(x * h ), respectively. They are rounded according to the format used in the table and the targeted accuracy for the final result.

To increase the accuracy of these tabulated values, Gal proposed a method to transfer some of the errors due to rounding over the reduced argument [START_REF]Computing elementary functions: a new approach for achieving high accuracy and good performance[END_REF]. This consists in introducing a small "corrective" term on the values x * h , denoted by corr here and hereafter. For each input entry x * h of the table, this term is carefully chosen to ensure that both sin(x * h +corr) and cos(x * h + corr) are very close to a floating-point number. In [START_REF] Gal | An accurate elementary mathematical library for the IEEE floating point standard[END_REF], Gal and Bachelis were able to find corr such that

sin h = sin(x * h + corr) • (1 + -(10+53n) ) and cos h = cos(x * h + corr) • (1 + -(10+53n)
). This corresponds to 10 extra bits of accuracy for both tabulated values of sine and cosine compared to Tang's tabulated values, thanks to a small perturbation corr on the input number x * . However, this imposes to store the corrective term corr along with the values sin h and cos h as the unevaluated sum of n floating-point numbers such that 10 + 53n is larger than the targeted accuracy. The value sin(x * ) is thus reconstructed as follows:

sin(x * ) = sin h • cos(x * l -corr) + cos h • sin(x * l -corr).
Gal's solution requires an exhaustive search in order to determine a "good" corrective term for each entry x * h . The search space grows exponentially with the number of extra bits for sin h and cos h . Stehlé and Zimmermann proposed an improvement based on LLL algorithm [START_REF] Stehlé | Gal's accurate tables method revisited[END_REF] to speed up the search. They were able to increase the accuracy of Gal's table by 11 extra bits.

C. Brisebarre et al.'s (M, p, k)-Friendly Points

In 2012, Brisebarre, Ercegovac, and Muller proposed a new method for sine and cosine evaluation with a few table lookups and additions in hardware [START_REF] Brisebarre | k)-friendly points: A table-based method for trigonometric function evaluation[END_REF], [START_REF] Wang | k)-friendly points: A table-based method to evaluate trigonometric function[END_REF]. Their approach consists in tabulating four values a, b, z, and x, defined as:

z = 1/ a 2 + b 2 and x = arctan(b/a),
where a and b are small particular integers. The reconstruction then corresponds to:

sin(x * ) = b • cos(x * -x) + a • sin(x * -x) • z.
The values (a, b, z, x) are chosen among specific points with integer coordinates (a, b) called (M, p, k)-friendly points. These points are recoded using canonical recoding which minimizes the number of non-zero digits in table entries. More precisely, a and b must be positive integers lower than M , such that the number z = 1/ √ a 2 + b 2 has less than k non-zero bits in the first p bits of its canonical recoding. These conditions make hardware multiplications by a, b, and z cheaper. Compared to other hardware evaluation schemes, this solution reduces by 50% the area on FPGAs for 24-bit accuracy.

III. A TABLE OF EXACT POINTS

Tang's tables store tabulated value of equally spaced distinguished points rounded to the destination format. This rounding error is problematic when seeking an accurate evaluation scheme since at least three double precision floatingpoint numbers are required. We have seen in Section II how Gal, Stehlé, and Brisebarre improved the accuracy of those tabulated values. This consists in finding input numbers whose images by the function are close to machine representable numbers. This artificially increases the accuracy of stored values by a few extra bits, but those values still embed some errors.

The proposed improvement is based on the following fact: After the first range reduction, the reduced number x * has to be considered as an irrational number, and therefore it has to be rounded to a given precision. We look for almost regularly spaced points whose images by the considered function will be exactly representable. The table will store exact values for sin h and cos h in a machine friendly format, plus a truncated corrective term to apply to the reduced argument as in Gal's method. This way, the error will be concentrated in the reduced number used in the polynomial evaluation x * l . We will now describe what good points for this method are and how the evaluation scheme can benefit of such points. For this purpose, let x be the input floating-point number and x * the reduced argument after the first range reduction such that

x * = x -k • π 2 , with k ∈ N * .
As π/2 is an irrational number, x * is irrational as well, and it has to be rounded to some precision j such that

x * = x * • (1 + -j ).
We should mention that x * is generally represented as the unevaluated sum of two or more floating-point numbers to reach an accuracy of j bits. As seen in Section I, the second range reduction splits x * into two parts, x * h and x * l , such that

x * = x * h + x * l and |x * l | < 2 -p .

The first p bits of x *

h are used to address the table T made of π/4 × 2 p entries. This table has to store values such that sin h and cos h can be recovered without any rounding error. This means that the following properties hold:

1) The value sin h and cos h have to be rational numbers such that

sin h = s n s d and cos h = c n c d , with s n , s d , c n , c d ∈ N * .
2) In order to make the reconstruction step easier, the denominators s d and c d have to be equal, that is,

s d = c d . 3)
In order to avoid unnecessary division during reconstruction, the common denominator has to be equal to the same k for every input number x * h . This value k can be seen as the lowest common multiple (LCM) of the denominator of all entries. 4) In order to minimize the table size, stored values s n and c n have to be representable as a machine word.

With numbers satisfying those properties, the reconstruction step corresponds to:

sin(x * ) = s n • P cos (x * l -corr) + c n • P sin (x * l -corr)
, where

• P cos (x) and P sin (x) are polynomial approximations, defined as follows for x ∈ [0, 2 -p ]:

P cos = cos(x) k and P sin = sin(x) k ,
• The tabulated values s n , c n , and corr are such that

x * h = arcsin s n k -corr = arccos c n k -corr with corr = corr•(1+ -53n
) corresponding to an approximation of the corrective term stored as the unevaluated sum of n floating-point numbers, • And s n and c n are stored exactly in a machine word (i.e double precision floating-point number).

As we can observe, we impose sin h and cos h to be rational numbers. This should requires to store both numerator and denominator, and to perform two divisions in the reconstruction step. However, we avoid the need to store the denominator and perform the associated division by integrating it in the polynomial approximation. This eliminates the cost of the division, and the error coming from this operation.

IV. BUILDING THE TABLE OF EXACT POINTS

The proposed table used to perform the second range reduction brings several benefits over existing solutions. However, building such a table of exact points is not trivial and it involves results from Pythagorean triples. These objects are described in the first part of this section. Then, we present a method to efficiently build tables of exact points for the sine and cosine functions.

A. Pythagorean Triples

Pythagorean triples are mathematical objects that have been widely studied in number theory. By definition a triple (a, b, c) of integers is a Pythagorean triple if and only if: A Pythagorean triple for which the fractions a/c and b/c are irreducible is called a primitive Pythagorean triple (PPT). This notion is essential since a PPT and its multiples refer to similar triangles, and thus are associated with the same angle θ. For example, the well known PPT (3, 4, 5) and all its multiples are associated to the angle θ = arcsin(3/5) ≈ 0.6435 rad. Recall that we are looking for rational numbers which approximate the sine and cosine for various angles between 0 and π/4. Therefore we will focus only on primitive Pythagorean triples.

a 2 + b 2 = c 2 ,

B. Construction of Subsets of Primitive Pythagorean Triples

The set of primitive Pythagorean triples is unbounded. Figure 1 shows all PPTs whose hypotenuse c is such that c ≤ 2 12 . We can observe that they cover a wide range of angles over ]0, π/2[, and that this set is finite.

In this article, we use the Barning-Hall tree [START_REF] Barning | On pythagorean and quasi-pythagorean triangles and a generation process with the help of unimodular matrices[END_REF] to build sets of PPTs that exploit its ternary-tree structure [START_REF] Price | The pythagorean tree: A new species[END_REF]. From a given PPT (a, b, c), represented as a column vector, three new PPTs are computed by multiplying the former with the following three constants matrices:

  1 -2 2 2 -1 2 2 -2 3   ,   -1 2 2 2 -1 2 -2 2 3   , and   1 2 2 2 1 2 2 2 3   .
It is proven that all PPTs can be generated from the root (3, 4, 5) with increasing hypotenuse lengths. For every generated PPT (a, b, c), we also consider its symmetric PPT (b, a, c), as it corresponds to a different angle in ]0, π/2[. For instance, the three children of the triple (3, 4, 5) using the Barning-Hall tree are:

(5, 12, 13), [START_REF] Stehlé | Gal's accurate tables method revisited[END_REF][START_REF] Muller | Elementary functions -algorithms and implementation[END_REF][START_REF] Wang | k)-friendly points: A table-based method to evaluate trigonometric function[END_REF], and (21, 20, 29),

and their symmetric counterparts are [START_REF] Defour | Cache-optimised methods for the evaluation of elementary functions[END_REF][START_REF] Defour | A new scheme for tablebased evaluation of functions[END_REF][START_REF]Computing elementary functions: a new approach for achieving high accuracy and good performance[END_REF], [START_REF] Muller | Elementary functions -algorithms and implementation[END_REF][START_REF] Stehlé | Gal's accurate tables method revisited[END_REF][START_REF] Wang | k)-friendly points: A table-based method to evaluate trigonometric function[END_REF], and (20, 21, 29).

C. Selection of Primitive Pythagorean Triples

Recall that the table is addressed by x * h , the p leading bits of the reduced argument x * . Then for each entry i of the table, we want to select only one primitive Pythagorean triple such that its corresponding angle θ is as close as possible to x * h , where: θ = arcsin(a/c) and x * h = i • 2 -p . This means that for each table entry, we have to select the triple that minimizes the corrective term corr defined by:

corr = x * h -arcsin(a/c).
Then, once the terms a, b, c, and corr are computed, a naive solution would consist in storing exactly a, b, c in double precision, and an approximation of corr on as many bits as x * l and computed in multiple precision. As presented in Section III, instead of doing this, the solution we propose consists in storing two integers A and B of the form

A = a c • k and B = b c • k
where k ∈ N * is the same for all table entries. In this article, in order to minimize the error coming from the multiplication, we are looking for table entries such that A and B are the smallest. This consists in looking for a value k that is small or even better the smallest.

Depending on the number of bits p for x * h , recall that we have a table with π/4 × 2 p entries. This corresponds to 101 entries for p = 7. For each considered entry, we need to have at least one PPT in order to find a set of A's and B's. We have reported in Figure 2(a) the number of PPTs per entry for p = 7, such that the generated hypotenuses were less than 2 18 . This corresponds to the minimum number of bits for the hypotenuse so that we have at least one PPT per entry. On this figure, we observe that we have a mean number of 413 PPTs per entry, and a standard deviation of 33. We are looking for the lowest common multiple (LCM) k of one PPT hypotenuse per table entry. This means that we have to compute the LCM for the set made of 101 elements which can take 413 different values. This corresponds to approximately 413 101 combinations. This number is definitely too large to be tested.

We can observe on Figure 2(a) that the number of PPTs for the table entry 0, that is, near the angle 0, is equal to 106 whereas the mean is around 413. We have reported in Figure 2(b), the number of PPTs per entry for p = 7, such that the generated hypotenuses were less than 2 14 . This corresponds to the minimum number of bits for the hypotenuse so that we have at least one PPT per entry, but when the first entry from the table is excluded from the search space. On this figure, we observe that we have a mean number of 26 PPTs per entry and a standard deviation of 4.6. Therefore, if the first entry is excluded, the search space is reduced to 26 100 combinations to test. There are still too many combinations to compute an LCM in a reasonable time. To further reduce the space search, the LCM will not be computed but rather be looked for among generated hypotenuses. Nevertheless, it is still possible to exclude the first entry by selecting the degenerated PPT (0, 1, 1) for this entry. This PPT corresponds to the angle θ = 0 • exactly with a corrective term corr = 0. The advantage of selecting this PPT is twofold: (1) First, the search space is reduced. (2) Second, its hypotenuse is equal to 1 which will not impact the search of LCM and the corrective term is exact.

V. IMPLEMENTATION AND NUMERIC RESULTS

In this section we present how we implemented the proposed method to generate tables of exact points that fulfill the properties listed in Section III. We have designed two solutions to look for a small common multiple k. The first solution is based on an exhaustive approach and allows us to build tables indexed by up to 7 bits in a reasonable amount of time. The second solution uses a heuristic approach which reduces memory consumption during the execution and the search time.

A. Exhaustive Search 1) Algorithm: As seen in Section IV-C, we need to look for a small common multiple k of a combination of one PPT per table entry. We designed a C++ program that takes as input the number p of bits used to index the table and look for the minimal number n of bits used to represent the hypotenuse that corresponds to k • c. The search is exhaustive amongst the generated hypotenuses, which warranties that if a common multiple is found, it is the lowest. We initialize n to 4 and we build the program around the following three steps:

1) Generate all PPTs (a, b, c) such that c ≤ 2 n . 2) Search for the LCM k among all generated hypotenuses c.

3) If such a k is found, build tabulated values (A, B, corr) for every entry. Otherwise, set n = n + 1 and go back to step 1.

Values (A, B, corr) are computed using the PPTs (a, b, c) such that k is a multiple of c. In case several PPTs per entry fulfill this property, the one for which arcsin(a/c) is the closest to x * h is selected as this minimizes the error on x * l . 2) Numeric Results: Table I shows the results we obtained on an Intel(R) Xeon(R) CPU E5-2650 v2 @ 2.6 GHz (32 cores) with 125 GB of RAM running on GNU/Linux. It describes for the number p of bits that is targeted, the value of k min that was found followed by the number n of bits used to represent k min , the time necessary to find this value, and the numbers of considered PPTs and hypotenuses. We observe that it is possible to find k min and to build tables indexed by up to p = 7 bits in a reasonable amount of time. However, during our test we noticed that the heap memory was heavily solicited and that it was not possible to go beyond 8 bits.

Table II describes the table of exact points for p = 5, where k = 160 225 and the absolute value of the corrective term is at most 135188896813938 × 2 -53 , that is, ≈ 0.0150090 for input index i = 11.

B. Heuristic Search

To build tables indexed by a larger number of bits, it is therefore mandatory to use another solution. In order to reduce the search space, we have developed a heuristic that selects "good" hypotenuses and rejects others during the PPT generation phase.

We collected information related to the decomposition in prime factors of each k min found using the exhaustive search. Such a decomposition is given in Table III. These factorizations show that every k in the table is a composite number divisible by relatively small primes. Furthermore all those small primes are of the form 4n + 1, better known as Pythagorean primes [START_REF] Fässler | Multiple Pythagorean number triples[END_REF].

Therefore, the heuristic we propose follows a simple rule: only store primitive Pythagorean triples whose hypotenuse is of the form:

c = i p ri i with p i ∈ P and r i ∈ {0, 1} if p i = 5 r i ∈ N * else , (2)
where P is the set of Pythagorean primes lower than or equal to 73: P = 5, 13, 17, 29, 37, 41, 53, 61, 73 .

Results, timings, and numbers of considered triples and hypotenuses for this heuristic are given in Table IV. We can observe that this algorithm considers a number of hypotenuses several order of magnitude lower than the exhaustive search solution. This reduces the amount of necessary memory and execution time. For instance, for p = 7, only 3308 triples are stored, compared to the 1 347 953 triples for the exhaustive algorithm. In this case, the execution time was reduced from 31 seconds to 0.4 seconds. With this heuristic, the bottleneck is no longer memory but the selection of PPTs during their generation. Indeed, checking if a given hypotenuse satisfies (2) requires checking if it is multiple of Pythagorean prime numbers, which is an expensive test.

VI. COMPARISONS WITH OTHER METHODS

We have presented a range reduction based on exact points and how to efficiently build those points. In order to compare this solution with the other solutions presented in Section II, we consider a two-phase evaluation scheme of the sine function that targets correct rounding in double precision. The quick and the accurate phases target a relative error less than 2 -66 and 2 -150 , respectively. We choose p = 10 which corresponds to π/4 × 2 10 = 805 entries in each considered table.

In order to ease comparisons, we are only considering the number of memory accesses required by the second range reduction and the number of floating-point operations involved in the reconstruction step. We will consider that expansion algorithms are used whenever high accuracy is required as it is the case in the correctly rounded library CR-Libm [START_REF]CR-Libm, a library of correctly rounded elementary functions in double-precision[END_REF]. An expansion of size n consists in storing a given number as the unevaluated sum of n floating-point numbers [START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF]. We will use Table V extracted from [START_REF] Lauter | Basic building blocks for a triple-double intermediate format[END_REF] as the reference cost for those algorithms when no fma (fused multiply-add) is available. The notation E n stands for expansion of size n, such that, with this notation, E 1 represents a regular floating-point number.

= E2 + E2 12 E3 = E3 + E3 27 E2 = E1 × E2 20 E3 = E1 × E3 47 E2 = E1 × E1 17 E2 = E2 × E2 26 E3 = E3 × E3 107
A. Tang's Solution

In order to reach an accuracy of 66 bits, Tang's solution requires accessing to tabulated values sin h and cos h that are stored as expansions of size 2. These values need to be multiplied by the two results of the polynomial evaluations, which can be considered stored as expansion of size 2 as well. The total cost of the quick phase becomes: 4 double precision memory accesses (MA), 2 multiplications E 2 × E 2 , and 1 addition E 2 + E 2 , that is, 64 floating-point operations (FLOP).

In case the quick phase failed to return correct rounding, the accurate phase is launched. This requires accessing to 2 extra tabulated values to represent sin h and cos h as expansions of size 3. Those values need to be multiplied by the two results of the polynomial evaluations, which can be considered stored as expansion of size 3 as well. The total cost of the accurate phase becomes: 2 extra memory accesses, 2 multiplications E 3 × E 3 , and 1 addition E 3 + E 3 . This corresponds to 6 memory accesses, 241 floating-point operations, and a computed table of 38 640 bytes.

B. Gal's Solution

Using Gal's method, the corrective term allows to achieve around 66 bits of accuracy, and Stehlé's improvement allows to reach 74 bits. In both cases, only one double precision number is required for sin h , cos h , and the corrective term. Again, these values need to be multiplied by the two results of the polynomial evaluations, which can be considered stored as expansions of size 2. Thus the quick phase requires 2 + 1 double precision memory accesses, 1 addition for the corrective term, 2 multiplications E 1 × E 2 and 1 addition E 2 + E 2 . The cost of the quick phase with this table becomes 3 memory accesses and 53 floating-point operations.

To reach the 150-bit accuracy required by the accurate phase, it is necessary to get 2 extra floating-point numbers for each tabulated values. The corrective term is then integrated in the final result using addition between expansions of size 3. The rest of the operations need to be done using size-3 expansions. The total cost for the accurate phase becomes: 6 extra memory accesses, 2 multiplications E 3 × E 3 and 2 additions E 3 + E 3 , that is, 9 memory accesses, 268 floatingpoint operations, for a 57 960-byte computed table.

C. The Exact Points Solution

With our solution, as shown in Table IV, at most 39 bits are required to store A and B, that is, only one floating-point number per entry. Hence, the cost of the quick phase is the same as Gal's approach in Section VI-B.

However, for the accurate phase, values sin h and cos h that were accessed during the quick phase are exact, and do not require any extra memory access. The corrective term is stored as an expansion of size 3 and it requires 2 extra floatingpoint numbers to reach the 150 bits of accuracy. The corrective term is integrated in the final result using an addition with expansion of size 3. Multiplications correspond to E 3 = E 1 × E 3 as the results of the polynomial evaluations are stored as expansions of size 3. The final addition is done using E 3 = E 3 + E 3 operation. That is, the total cost of this step becomes 5 memory accesses and 148 floating-point operations, for a computed table of 32 200 bytes. Second, regarding the number of operations for both the quick and accurate phases, we observe that our solution provides similar performance to Gal's solution for the quick phase. For the accurate phase, we observe an improvement in favor of our approach of ≈ 39% and 45% over Tang and Gal, respectively.

D. Comparison Results

Third, we notice that the proposed solution reduces the number of memory accesses. The quick phase requires 3 accesses, while Tang's approach uses 4 accesses, that is, an improvement of 25%. The benefit is even more significant for the accurate phase, as the number of accesses is reduced from 9 to 4 compared to Gal's approach. This is an improvement of ≈ 45%.

VII. CONCLUSIONS AND PERSPECTIVES

In this article, we have presented a new approach to address the second range reduction of the evaluation process of elementary functions based on tabulated values. It relies on Pythagorean triples, which allow to simplify and accelerate the evaluation of these functions. Compared to other solutions, the table of Exact Points eliminates the rounding error on certain tabulated values, and transfers this error in the remaining reduced argument. We have focused the method on sine and cosine functions as it is the most difficult. However, the concept remains valid for other functions. For those functions and thanks to the proposed method, it is possible to reduce up to 45% the table sizes, the number of memory accesses and the number of floating-point operations involved in the reconstruction process.

Our further research direction is threefold: First, it would be interesting to plug the computed table into a full sine and cosine implementation to observe its actual impact. This could be done within the CR-Libm project. Second, as we have seen in Section V-B, relevant hypotenuses are factors of small Pythagorean primes. Following this, it would be interesting to characterize "good" hypotenuses, so that, instead of generating a huge set of triples and then choosing the relevant ones, they could be computed directly possibly using Bresenham's algorithm. Doing this way, it would speed up the table computation process. Third, we have focused on looking for the lowest common multiple so that tabulated values would be stored on the minimal number of bits. This property is essential for hardware implementations where each bit is important. For software solution, it would be more interesting to look for tabulated values that could fit in a given format (i.e double precision floating-point number) but that would minimize the corrective terms. This may result in corrective terms that could be representable as expansions of size strictly less than n, and thus save some extra memory accesses and associated floating-point operations.

4 Figure 1 .

 41 Figure 1. Primitive Pythagorean triples whose hypotenuse c is less than 2 12 .

  with a, b, c ∈ N * . It follows from the Pythagorean theorem that such a triple corresponds to the lengths of a right triangle edges. Sine and cosine of right triangle are defined as quotients of these lengths. Hence, given any Pythagorean triple (a, b, c), there exists an angle θ ∈]0, π/2[ such that: sin(θ) = a c and cos(θ) = b c .

  When entry 0 is omitted.

Figure 2 .

 2 Figure 2. Minimum number of PPTs per entry for p = 7.

Table I .

 I EXHAUSTIVE SEARCH RESULTS.

	p	kmin	n	time (s)	Triples	Hypotenuses
	3	425	9	1	86	66
	4	5525	13	1	1404	889
	5	160 225	18	0.2	42 328	24 228
	6	1 698 385	21	7	335 344	179 632
	7	6 569 225	23	31	1 347 953	686 701
	8	> 2 27	> 27	> 6700?	> 21 407 992	> 10 144 723

Table II .

 II VALUES (A, B, corr) COMPUTED FOR p = 5.

	Input index	A	B	corr
	0	0	160225	+0x0.0000000000000p+0
	1	5772	160121	+0x1.3966f27bfc9f0p-8
	2	11385	159820	+0x1.1a56677e409a8p-7
	3	15225	159500	+0x1.73400355ad200p-10
	4	19775	159000	-0x1.4b6e4ab5496c0p-10
	5	24505	158340	-0x1.62b596fa18d80p-9
	6	30044	157383	+0x1.27ac440de0a80p-10
	7	33820	156615	-0x1.8df1bd239dc40p-8
	8	39440	155295	-0x1.522b2a9e84900p-10
	9	44863	153816	+0x1.4d7623b1c6e80p-9
	10	50375	152100	+0x1.e0220454e0100p-8
	11	56257	150024	+0x1.ebd078b04dc80p-7
	12	58305	149240	-0x1.4eccbaa9e7000p-9
	13	63504	147103	+0x1.4f7dfcae44600p-10
	14	67860	145145	-0x1.53f734851f800p-13
	15	73593	142324	+0x1.158024b185460p-7
	16	75400	141375	-0x1.49140da6fe460p-7
	17	81345	138040	+0x1.48c11cc509200p-10
	18	85255	135660	-0x1.75ed3145ed600p-10
	19	90596	132153	+0x1.d824b6b5b8300p-8
	20	93247	130296	-0x1.f7e24ff929500p-9
	21	97760	126945	-0x1.5696f82020000p-17
	22	101500	123975	-0x1.7caa112f28800p-10
	23	105465	120620	-0x1.29d074350d800p-12
	24	110500	116025	+0x1.68ddad9df7000p-7
	25	111969	114608	-0x1.eb6d097519180p-8

Table III .

 III PRIME FACTORIZATION OF FOUND COMMON MULTIPLES.

	k	Prime factorization
	425	5 2 • 17
	5525	5 2 • 13 • 17
	160 225	5 2 • 13 • 17 • 29
	1 698 385	5 • 13 • 17 • 29 • 53
	6 569 225	5

2 • 13 • 17 • 29 • 41

Table IV .

 IV HEURISTIC SEARCH RESULTS.

	p	n	kmin	time (s)	Triples	Hypotenuses
	3	9	425	1	41	8
	4	13	5525	1	210	17
	5	18	160 225	1	995	40
	6	21	1 698 385	0.1	2171	66
	7	23	6 569 225	0.4	3452	69
	8	29	314 201 225	9.5	10 467	100
	9	34	12 882 250 225	294	20 311	109
	10	39	279 827 610 985	9393	33 056	110

Table V .

 V COST OF ADDITION AND MULTIPLICATION OF EXPANSIONS.

	Operation	Floating-point operations
	E2

Table VI .

 VI Table VI synthesizes the estimated costs for those three range reduction algorithms based on tabulated values. This table reports the number of floating-point operations for the quick COMPARISONS BETWEEN THREE TABLE-BASED RANGE REDUCTIONS, FOR p = 10. THE NUMBER OF MEMORY ACCESSES (MA) AND THE NUMBER OF FLOATING POINT OPERATIONS (FLOP) ARE First, we notice from this table that the proposed tablebased range reduction requires less memory per table entry than others tested solutions. Tang's method needs 48 bytes per entry and Gal's method, 72 bytes per entry, while 40 bytes are enough for the table with exact points. This is an improvement in memory usage of ≈ 17% and ≈ 45% over Tang's and Gal's methods, respectively.

			REPORTED.
	Solution	Quick phase	Accurate phase	Table size (bytes)
	Tang	4 MA + 64 FLOP	6 MA + 241 FLOP	38640
	Gal	3 MA + 53 FLOP	9 MA + 268 FLOP	57960
	Proposed	3 MA + 53 FLOP	5 MA + 148 FLOP	32200
	and accurate phases, together with the size in bytes of the
	computed table.

ACKNOWLEDGMENT

This work was supported by the ANR MetaLibm project (Ingénierie Numérique et Sécurité 2013, ANR-13-INSE-0007).


Razvan Gabriel Iagar
email: razvan.iagar@icmat.es

Philippe Laurençot
email: philippe.laurencot@math.univ-toulouse.fr

Large time behavior for a quasilinear diffusion equation with critical gradient absorption

Keywords: AMS Subject Classification: 35K59, 35K65, 35K92, 35B40 large time behavior, degenerate diffusion, gradient absorption, decay estimates, logarithmic scales

We study the large time behavior of non-negative solutions to the nonlinear diffusion equation with critical gradient absorption

for p ∈ (2, ∞) and q * := p -N/(N + 1). We show that the asymptotic profile of compactly supported solutions is given by a source-type self-similar solution of the p-Laplacian equation with suitable logarithmic time and space scales. In the process, we also get optimal decay rates for compactly supported solutions and optimal expansion rates for their supports that strongly improve previous results.

Introduction and results

The goal of the present paper is to study the large time behavior of non-negative solutions to the following equation which combines a nonlinear diffusion and a gradient absorption term:

∂ t u -∆ p u + |∇u| q * = 0 in (0, ∞) × R N , (1.1)
u(0) = u 0 in R N , (1.2)
where the p-Laplacian ∆ p u is given by ∆ p u := div(|∇u| p-2 ∇u) for p ∈ (2, ∞), and the critical exponent is q * := p -N/(N + 1) > 1. We consider the Cauchy problem (1.1)-(1.2) for initial data u 0 satisfying

u 0 ∈ W 1,∞ (R N ), u 0 ≥ 0, u 0 ≡ 0, supp u 0 is compact in R N . (1.
3)

The theory for the more general equation

∂ t u -∆ p u + |∇u| q = 0, in (0, ∞) × R N ,
with general values of the exponents p and q, developed very quickly in the last years, after understanding better how the competition between the two effects-nonlinear diffusion in form of a p-Laplacian term, and first order absorption-works with respect to different ranges of these exponents.

The semilinear problem, when p = 2, has been understood first, due to the possibility of using, at least partially, techniques related to the heat equation. It has been noticed that there exists two critical values of the exponent q, namely q = 1 and q = q * = (N + 2)/(N + 1), and the qualitative theory together with the large time behavior are now well understood after a series of works [START_REF] Andreucci | The Cauchy problem for degenerate parabolic equations with source and damping[END_REF][START_REF] Ben-Artzi | T he local theory for viscous Hamilton-Jacobi equations in Lebesgue spaces[END_REF][START_REF] Benachour | Global solutions to viscous Hamilton-Jacobi equations with irregular initial data[END_REF][START_REF] Bidaut-Véron | L ∞ estimates and uniqueness results for nonlinear parabolic equations with gradient absorption terms[END_REF][START_REF] Biler | Asymptotic properties of solutions of the viscous Hamilton-Jacobi equation[END_REF][START_REF] Gilding | The Cauchy problem for u t = ∆u + |∇u| q , large-time behaviour[END_REF][START_REF] Gilding | The Cauchy problem for u t = ∆u + |∇u| q[END_REF]. More precisely, for q > q * we have an asymptotic simplification with dominating diffusion, meaning that the absorption plays no role in the large time behavior, while for 1 < q < q * , there exists a special, unique self-similar solution of very singular type, which is the asymptotic profile for the evolution [START_REF] Benachour | Asymptotic profiles of solutions to viscous Hamilton-Jacobi equations[END_REF]. The critical case corresponding to q = q * requires a different treatment and has been investigated in [START_REF] Th | Asymptotic behavior for a viscous Hamilton-Jacobi equation with critical exponent[END_REF] using central manifold theory techniques. Finally, the case q = 1 is the subject of a well-known open problem (see [START_REF] Benachour | Asymptotic estimates of solutions of u t -1 2 ∆u = -|∇u| in R + × R d , d ≥ 2[END_REF] for the best estimates available), while for q ∈ (0, 1), finite time extinction occurs [START_REF] Benachour | Extinction and decay estimates for viscous Hamilton-Jacobi equations in R N[END_REF][START_REF] Benachour | Extinction and nonextinction for viscous Hamilton-Jacobi equations in R N , Asympt[END_REF][START_REF] Gilding | The Cauchy problem for u t = ∆u + |∇u| q , large-time behaviour[END_REF]. As heat equation techniques are not anymore available for p > 2, its study is more involved and it came later. In this range, the qualitative behavior is very different: on the one hand, the nonlinear slow diffusion implies a finite speed of propagation which entails that solutions emanating from compactly supported initial data stay compactly supported for all times. On the other hand, there exists a range of exponents q ∈ (1, p -1] where the dynamics is governed by the Hamilton-Jacobi part, a new fact that does not appear for p = 2 and q > 1; indeed, for q ∈ (1, p -1], the large time behavior is given by profiles in form of "sandpiles" or "regularized sandpiles", reminding of the solutions to Hamilton-Jacobi equations, see [START_REF] Iagar | Asymptotic behavior of a nonlinear parabolic equation with gradient absorption and critical exponent[END_REF][START_REF] Ph | Localized non-diffusive asymptotic patterns for nonlinear parabolic equations with gradient absorption[END_REF]. For higher values of the exponent q, asymptotic simplification with dominating diffusion is expected for q > q * = p -N/(N + 1) and the existence of very singular self-similar solutions is known for q ∈ (p -1, q * ) [START_REF] Shi | Self-similar singular solution of a p-Laplacian evolution equation with gradient absorption term[END_REF], the lack of a uniqueness result preventing a complete understanding of the large time asymptotics in that case. In the present work we complete the panorama of the slow diffusion case by studying the case q = q * , where finer estimates for the solution than the general ones in [START_REF] Bartier | Gradient estimates for a degenerate parabolic equation with gradient absorption and applications[END_REF] are needed.

We finally mention that, for p ∈ (1, 2), a qualitative theory is developed starting from the paper [START_REF] Iagar | Positivity, decay and extinction for a singular diffusion equation with gradient absorption[END_REF], where a new critical exponent q = p/2 is discovered, and the qualitative theory for any q > 0 is established. In particular, there are again critical exponents q = q * and q = p/2, limiting ranges of parameters with different behaviors: diffusive for q > q * , algebraic decay for q ∈ (q * , p/2), exponential decay for q = p/2, finite time extinction for q < p/2. Recent work by the authors helped to establish the existence of special solutions and the large time behavior for some of these cases, such as q ∈ (q * , p/2), with very singular self-similar solutions, see [START_REF] Iagar | Existence and uniqueness of very singular solutions for a fast diffusion equation with gradient absorption[END_REF][START_REF] Iagar | Asymptotic behavior for a singular diffusion equation with gradient absorption[END_REF], and q = p/2 where eternal solutions have been discovered [START_REF] Iagar | Eternal solutions to a singular diffusion equation with critical gradient absorption[END_REF]. The critical case q = q * is still open in this range.

Main results. Coming back to Eq. (1.1), our goal is to determine a profile that the solutions approach as t → ∞. As we are in a critical case which plays the role of an interface between purely diffusive behavior and mixed-type behavior, some logarithmic time scales are expected to appear. We introduce the following constant which will be important in the analysis

η := 1 N (p -2) + p = 1 p(N + 1) -2N . (1.4)
We are now in a position to state our main asymptotic result.

Theorem 1.1. Consider an initial condition u 0 satisfying (1.3), q = q * = p -N/(N + 1), and let u be the solution to the Cauchy problem (1.1)-(1.2) with initial condition u 0 . Then 1) , and A * is uniquely determined and given by:

lim t→∞ t N η (log t) p(N +1)η u(t, x) - 1 t N η (log t) p(N +1)η B A * x t η (log t) (2-p)(N +1)η = 0, where B A (y) = A -B 0 |y| p/(p-1) (p-1)/(p-2) + , B 0 = p -2 p η 1/(p-
A * :=     (N + 1) ∞ 0 1 -B 0 r p/(p-1) (p-1)/(p-2) + r N -1 dr η q/(p-1) ∞ 0 1 -B 0 r p/(p-1) q/(p-2) + r N -1+q/(p-1) dr     p(p-2)(N +1)η/(p-1)
.

Let us remark that B A is the well-known Barenblatt profile, the fundamental solution (and also the asymptotic profile) of the pure diffusion equation (without the gradient absorption term). The effect of absorption is seen in the fact that both the time-decay rate and the expansion of supports are changed with respect to the p-Laplacian equation.

In order to prove Theorem 1.1, we have noticed that the precise time decay rate of the solutions was missing from the theory. Indeed, it follows from [START_REF] Bartier | Gradient estimates for a degenerate parabolic equation with gradient absorption and applications[END_REF] that the L 1 -norm of the solutions to the Cauchy problem (1.1)-(1.2) with initial condition u 0 satisfying (1.3) decays with a logarithmic rate as t → ∞; however, as we will prove, the rate in [START_REF] Bartier | Gradient estimates for a degenerate parabolic equation with gradient absorption and applications[END_REF] is not optimal and it can be improved. Based on that, we find the exact logarithmic correction for the profile. Once optimal estimates are obtained, the final step will be an application of the stability technique of Galaktionov and Vázquez [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF][START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF] together with some analysis concerning the variation of the L 1 -norm of the solution, which is needed to establish the uniqueness (and precise parameter A = A * ) of the asymptotic profile.

A by-product of our analysis is the following expansion property of the positivity set of non-negative compactly supported solutions to (1.1).

Proposition 1.2. Consider an initial condition u 0 satisfying (1.3), q = q * = p-N/(N +1), and let u be the solution to the Cauchy problem for the equation (1.1) with initial condition u 0 . Introducing the positivity set

P(t) := {x ∈ R N : u(t, x) > 0} of u at time t ≥ 0, there are ̺ 2 > ̺ 1 > 0 such that B 0, ̺ 1 t η (log t) (2-p)(N +1)η ⊂ P(t) ⊂ B 0, ̺ 2 t η (log t) (2-p)(N +1)η , t ≥ 1 .
Organisation of the paper. We begin by establishing our optimal time decay rate for the solutions to (1.1), together with the rate of expansion of supports. As we already mentioned, they strongly improve existing results and we consider them as the main theoretical novelty of the present work. This is done in Section 2, making use of a fine analysis of the supports and functional inequalities of Hölder and Poincaré type. The precise estimates are given in Proposition 2.1 and Corollary 2.3. These sharp estimates give in particular the correct time scales, thus allowing for a rescaling to reach a new equation with solutions being uniformly bounded in

L 1 (R N ) ∩ W 1,∞ (R N ).
The scaling step is performed in Section 3, where we also show that the Barenblatt profile itself provides a suitable subsolution, which is fundamental in the sequel in avoiding the phenomenon of collapse to zero in the limit. Finally, in Section 4, we complete the proof of Theorem 1.1 by using the stability theorem for small perturbations of dynamical systems from [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF][START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF], which we briefly recall in a more abstract framework in the Appendix. This approach is by now rather standard; nevertheless, we present it in some details, as the presence of the gradient absorption term leads to some further technical complications.

Optimal L ∞ -bounds

Let u 0 be an initial condition satisfying (1.3) and denote the corresponding solution to the Cauchy problem (1.1) by u. We recall that the existence and uniqueness of a non-negative (viscosity) solution u ∈ BC([0, ∞) × R N ) to (1.1) which is also a weak solution follows from [2, Theorem 1.1], and moreover it satisfies

0 ≤ u(t, x) ≤ u 0 ∞ , ∇u(t) ∞ ≤ ∇u 0 ∞ , for any t ≥ 0. In addition, it follows from [1, Theorem 1.1 & Theorem 1.2] and [2, Theorem 1.6 & Corollary 1.7] that u(t)
is compactly supported for any t > 0 so that the maximal radius of the (compact) support of u(t) and u 0 defined by

̺(t) := inf{R > 0 : u(t, x) = 0 for |x| > R}, ̺ 0 := ̺(0),
are finite for each t ≥ 0. Furthermore, the following estimates for the support and the L 1 -norm are established in the above mentioned references: there is C 1 > 0 such that

̺(t) ≤ C 1 (1 + t) η for any t ≥ 0 (2.1) and u(t) 1 ≤ C 1 log(1 + t) -1/(q * -1) for any t ≥ 1. (2.2)
However, as we shall see below, the estimates (2.1) and (2.2) are far from being optimal and can be strongly improved. We finally recall that, due to [1, Lemma 4.1], the solution u to (1.1) cannot vanish identically in finite time, that is,

u(t) ≡ 0 for t ≥ 0 .
Proposition 2.1. There exists a positive constant C 2 , depending only on p, q, N and u 0 , such that for any t > 0,

̺(t) ≤ C 2 (1 + t) η (log(1 + t)) -η(p-2)(N +1) , (2.3
)

and u(t) 1 ≤ C 2 (log(1 + t)) -(N +1) . (2.4)
Notice that this is a real improvement in the second estimate (2.4) with respect to (2.2), since

1 q * -1 = N + 1 p(N + 1) -(2N + 1) < N + 1 2(N + 1) -(2N + 1) = N + 1, for p > 2.
Proof. Since u 0 is non-negative, continuous and compactly supported in R N , there exists a non-negative, continuous and radially symmetric, radially non-increasing function U 0 with compact support, such that 0 ≤ u 0 (x) ≤ U 0 (x) for any x ∈ R N . Let U be the corresponding solution to (1.1) with initial data U 0 ; it follows from [2, Theorem 1.2] that the function x → U (t, x) is also radially symmetric, radially non-increasing and compactly supported, for any t > 0, and we deduce from the comparison principle that 0 ≤ u(t, x) ≤ U (t, x), for any (t, x) ∈ (0, ∞) × R N . Let, for t ≥ 0,

σ(t) := sup x∈R N {|x| : U (t, x) > 0}
be the radius of the support of U . We then have ̺(t) ≤ σ(t). Furthermore, we infer from [2, Theorem 1.2 & Proposition 1.4] that

|∇U (p-2)/(p-1) (t, x)| ≤ C 3 U (s) (p-2)/p(p-1) ∞ (t -s) -1/p , for 0 ≤ s < t (2.5) and U (t) ∞ ≤ C 3 U (s) pη 1 (t -s) -N η , for 0 ≤ s < t. (2.6)
It next follows from (1.1) that, for any non-negative function y ∈ C 1 ([0, ∞)), we have

d dt {|x|≥y(t)} U (t, x) dx = {|x|≥y(t)} ∂ t U (t, x) dx -y ′ (t) {|x|=y(t)} U (t, x) dx ≤ {|x|≥y(t)} div(|∇U | p-2 ∇U )(t, x) dx -y ′ (t) {|x|=y(t)} U (t, x) dx ≤ - {|x|=y(t)} |∇U (t, x)| p-2 ∇U (t, x) • x |x| dx -y ′ (t) {|x|=y(t)} U (t, x) dx ≤ {|x|=y(t)} p -1 p -2 |∇U (p-2)/(p-1) (t, x)| p-1 -y ′ (t) U (t, x) dx.
Fix now t 0 ≥ 0. For t > t 0 we deduce from the estimates (2.5) (with s = (t + t 0 )/2) and (2.6) (with ((t + t 0 )/2, t 0 ) instead of (t, s)) that

p -1 p -2 |∇U (p-2)/(p-1) (t, x)| p-1 ≤ 2 (p-1)/p p -1 p -2 C p-1 3 U t + t 0 2 p-2)/p ∞ (t -t 0 ) -(p-1)/p ≤ (2C p 3 ) (p-1)/p p -1 p -2 C 3 U (t 0 ) pη 1 t -t 0 2 -N η (p-2)/p (t -t 0 ) -(p-1)/p ≤ ηC 4 U (t 0 ) (p-2)η 1 (t -t 0 ) η-1 . Choosing y(t) := σ(t 0 ) + C 4 U (t 0 ) (p-2)η 1 (t -t 0 ) η ,
the above inequality reads

p -1 p -2 |∇U (p-2)/(p-1) (t, x)| p-1 ≤ y ′ (t),
for any t ≥ t 0 , from which we deduce that

d dt {|x|≥y(t)} U (t, x) dx ≤ 0 for t ≥ t 0 .
Since y(t 0 ) = σ(t 0 ), we end up with

{|x|≥y(t)} U (t, x) dx ≤ {|x|≥σ(t 0 )} U (t 0 , x) dx = 0,
for any t > t 0 . Owing to the non-negativity of U , this is only possible if U (t, x) ≡ 0 for |x| ≥ y(t), which means

σ(t) ≤ y(t) = σ(t 0 ) + C 4 U (t 0 ) (p-2)η 1 (t -t 0 ) η , (2.7)
for t > t 0 . In order to proceed further, we need the following Poincaré inequality:

Lemma 2.2. Given µ ∈ [1, ∞)
and R > 0 there exists a constant K depending only in N and µ such that

w µ L 1 (B(0,R)) ≤ KR µ(N +1)-N ∇w µ L µ (B(0,R)) for all w ∈ W 1,µ 0 (B(0, R)) .
Proof of Lemma 2.2. This follows by combining Hölder and Poincaré inequalities. More precisely, setting µ ′ := µ/(µ -1), one has

w µ L 1 (B(0,R)) ≤ w µ L µ (B(0,R)) |B(0, R)| µ/µ ′ = K w µ L µ (B(0,R)) R N (µ-1) ≤ KR µ ∇w µ L µ (B(0,R)) R N µ-N = KR µ(N +1)-N ∇w µ L µ (B(0,R)) ,
where we have used that the Poincaré constant is of order O(R).

In order to continue the proof of Proposition 2.1, we first note that integrating (1.1) over R N gives d dt U (t) 1 + ∇U q * q * = 0 .

Since U (t) is supported in B(0, σ(t)) for each t ≥ 0, it belongs to W 1,q * 0 (B(0, σ(t))) and we apply Lemma 2.2 to obtain

d dt U (t) 1 + 1 K U (t) q * 1 σ(t) 1/η ≤ 0 (2.8)
for any t ≥ 0. Now, fix T ≥ 1 and introduce the notation

Σ(T ) := sup t∈[1,T ] t -η (log t) A σ(t) , A := η(p -2)(N + 1).
We infer from (2.8) that, for any t ∈ [1, T ],

d dt U (t) 1 + (log t) (p-2)(N +1) Kt U (t) q * 1 (t -η (log t) A σ(t)) 1/η ≤ 0, whence d dt U (t) 1 + (log t) (p-2)(N +1) Kt(Σ(T )) 1/η U (t) q * 1 ≤ 0 , t ∈ [1, T ] .
Integrating the above inequality over (1, t), t ∈ (1, T ), we find

U (t) 1 ≤ (K(N + 1)) 1/(q * -1) Σ(T ) 1/η(q * -1) (log t) -(N +1) , t ∈ (1, T ]. (2.9)
Next, let m ≥ 1 be an integer to be determined later and consider t ∈ (1, T ].

• Either t ≤ 2 m and it follows from (2.7) with t 0 = 1 that

t -η (log t) A σ(t) ≤ (log 2 m ) A σ(1) + C 4 U (1)
(p-2)η 1

2 mη ≤ C(m) .
• Or 2 m ≤ t ≤ T (if 2 m ≤ T ), and we infer from (2.7) (with t 0 = t/2) and (2.9) that

t -η (log t) A σ(t) ≤ t -η (log t) A σ t 2 + C 4 t -η (log t) A U t 2 (p-2)η 1 t 2 η ≤ 2 -η t 2 -η log t log(t/2) A log t 2 A σ t 2 + C 4 2 -η (log t) A Σ(T ) 1/η(q * -1) log t 2 -(N +1) (p-2)η . Since t ≥ 2 m , we obtain that log t ≤ m m -1 log t 2 ,
whence, plugging this estimate in the previous inequality, we get 1) .

t -η (log t) A σ(t) ≤ m m -1 A 2 -η Σ(T ) + C 4 Σ(T ) (p-2)/(q * -
Combining the previous two estimates and taking the supremum over t ∈ (1, T ], we obtain

Σ(T ) ≤ m m -1 A 2 -η Σ(T ) + C 5 Σ(T ) (p-2)/(q * -1) + C(m).
We now fix m large enough such that

δ m := m m -1 A 2 -η < 1 .
Taking into account that p -2 < q * -1, we deduce from Young's inequality that

Σ(T ) ≤ δ m Σ(T ) + 1 -δ m 2 Σ(T ) + C(m) ,
which readily implies that Σ(T ) ≤ C(m) for each T ≥ 1, the constant C(m) being independent of T . We have thus proved that

σ(t) ≤ Ct η (log t) -η(p-2)(N +1) , t ≥ 1,
while the uniform bound for Σ(T ) together with (2.9) imply

U (t) 1 ≤ C(log t) -(N +1) .
Recalling that u(t, x) ≤ U (t, x) for any (t, x) ∈ (0, ∞) × R N and that ̺(t) ≤ σ(t), we obtain the expected estimates (2.3) and (2.4).

Corollary 2.3.

There exists a constant C 5 > 0 depending only on p, q, N and u 0 , such that for any t > 0, we have

u(t) ∞ ≤ C 5 (1 + t) -N η (log(1 + t)) -pη(N +1) (2.10) and ∇u(t) ∞ ≤ C 5 (1 + t) -(N +1)η (log(1 + t)) -2η(N +1) .
(2.11)

Proof. We combine the estimates in [2, Proposition 1.4] with the previous estimates of Proposition 2.1. We thereby obtain

u(t) ∞ ≤ C u t 2 pη 1 t -N η 2 N η ≤ Ct -N η log 1 + t 2 -(N +1)pη ≤ Ct -N η 1 2 log(1 + t) -(N +1)pη ≤ C(1 + t) -N η (log(1 + t)) -(N +1)pη and ∇u(t) ∞ ≤ C u t 2 2η 1 t -η(N +1) ≤ Ct -η(N +1) (log(1 + t)) -2η(N +1) ,
which completes the proof.

As we shall see in the next section, these estimates are optimal and are building blocks in identifying the large time behavior of compactly supported solutions to (1.1) and thus proving Theorem 1.1.

Scaling variables

Let u 0 be an initial condition satisfying (1.3) and denote the corresponding solution to the Cauchy problem (1.1)-(1.2) by u. According to the estimates derived in Section 2, we introduce the following new variables (s, y) and function w: s = log(e + t), y = x(e + t) -η log(e + t) (p-2)(N +1)η , (

and u(t, x) = (e + t) -N η log(e + t) -pη(N +1) w(s, y).

(3.2)

By (1.1) the rescaled function w solves

∂ s w -Lw = 0 in (1, ∞) × R N , (3.3a
)

w(1) = u 0 in R N , (3.3b)
where

Lz := ηy • ∇z + ηN z + ∆ p z - 1 s [|∇z| q * -pη(N + 1)z + (p -2)η(N + 1)y • ∇z] . (3.3c)
For further use, we introduce the autonomous counterpart of (3.3a) which is

∂ s v -Lv = 0 in (1, ∞) × R N , (3.4a)
with Lz := ηy • ∇z + ηN z + ∆ p z.

(3.4b)

The boundedness of w readily follows from the previous section.

Lemma 3.1.

There is a positive constant C 6 depending only on p, q, N , and u 0 such that

w(s) 1 + w(s) ∞ + ∇w(s) ∞ ≤ C 6 ,
for any s > 1. Moreover, the support of w(s) is localized: there exists R 0 > 0 such that supp(w(s)) ⊆ B(0, R 0 ) for any s ∈ (1, ∞).

Proof. The bounds for the W

w(s) 1 = (e + t) N η log(e + t) pη(N +1) R N u t, y(e + t) η log(e + t) -(p-2)(N +1)η dy = log(e + t) N +1 u(t) 1 ≤ C 2 .
Finally, the assertion about the localization of the support follows from estimate (2.3) and the definition of the new variable y in (3.1) (with R 0 = C 2 ).

Lemma 3.1 provides a fine upper bound for w but its optimality can only be guaranteed by a lower bound of the same order. In addition, such a lower bound would prevent the possibility of collapsing to the trivial solution in the limit s → ∞. This is a consequence of the following result.

Lemma 3.2.

There exists A sub > 0 sufficiently small such that the Barenblatt profile 1) , is a subsolution to Eq. (3.3a) for A ∈ (0, A sub ).

B A (y) = A -B 0 |y| p/(p-1) (p-1)/(p-2) + , B 0 = p -2 p η 1/(p-
Proof. A simple computation shows that LB A = 0 for any A > 0, where L is the autonomous operator defined in (3.4b). Moreover,

∇B A (y) • y = - B 0 p p -2 A -B 0 |y| p/(p-1) 1/(p-2) + |y| p/(p-1) , hence -sLB A = |∇B A | q * + η(p -2)(N + 1)y • ∇B A -η(N + 1)pB A = B 0 p p -2 q *
A -B 0 |y| p/(p-1) q * /(p-2)

+ |y| q * /(p-1) -η(N + 1)pB 0 A -B 0 |y| p/(p-1) 1/(p-2) + |y| p/(p-1)
-η(N + 1)p A -B 0 |y| p/(p-1) (p-1)/(p-2) + = η q * /(p-1) A -B 0 |y| p/(p-1) q * /(p-2)

+ |y| q * /(p-1) -η(N + 1)pA A -B 0 |y| p/(p-1) 1/(p-2) + = η A -B 0 |y| p/(p-1) 1/(p-2)
+ × η (q * -p+1)/(p-1) A -B 0 |y| p/(p-1) (q * -1)/(p-2) + |y| q * /(p-1) -(N + 1)pA .

Either |y| ≥ (A/B 0 ) (p-1)/p and -sLB A (y) = 0. Or |y| ≤ (A/B 0 ) (p-1)/p and, since A -B 0 |y| p/(p-1) ≤ A, we find

-sLB A ≤ η A -B 0 |y| p/(p-1) 1/(p-2) + η (q * -p+1)/(p-1) A (q * -1)/(p-2) A B 0 q/p -(N + 1)pA ≤ ηA A -B 0 |y| p/(p-1) 1/(p-2) + η (q * -p+1)/(p-1) B -q * /p 0 A θ -(N + 1)p ,
where

θ := q * -1 p -2 + q p -1 = (N + 1)(p -1)((N + 1)p -2N ) p(p -2) > 0
for p > 2, taking into account that q * = p -N/(N + 1). Consequently, there exists

A sub = (N + 1)pB q/p 0 η -(q-1)/(p-1) 1/θ > 0, such that -sLB A ≤ 0 in (1, ∞) × R N for A ∈ (0, A sub ), ending the proof.
A first consequence of Lemma 3.2 is the optimality of the temporal decay estimates established in Lemma 3.1.

Proposition 3.3.

There are A 0 ∈ (0, A sub ) and t 0 > 0 such that

u(t) ∞ ≥ A 0 (e + t) -N η log(e + t) -pη(N +1) , t > 0 , (3.5)
and

u(t 0 , x) ≥ B A 0 (x) , x ∈ R N . (3.6)
Proof. Since u 0 ≡ 0 and is continuous by (1.3), there are x 0 ∈ R N , ̺ 0 > 0, and

ε 0 > 0 such that u 0 (x) ≥ ε 0 , x ∈ B(x 0 , ̺ 0 ) . (3.7)
Introducing ũ0 (x) := u 0 (x + x 0 ), x ∈ R N , and denoting the corresponding solutions to (1.1) and (3.3a) by ũ and w, respectively, we observe that the invariance of (1.1) with respect to translations entails that ũ(t, x) = u(t, x + x 0 ), while we infer from (3.7) that

ũ0 (x) ≥ B A (x) , x ∈ R N , (3.8)
provided

A ≤ ε (p-2)/(p-1) 0 and A ≤ B 0 ̺ p/(p-1) 0 . (3.9)
Choosing A ∈ (0, A sub ) satisfying (3.9), we deduce from Lemma 3.2, (3.8), and the comparison principle that w(s, y)

≥ B A (y) , (s, y) ∈ (1, ∞) × R N .
Coming back to u, we realize that (e + t) N η log(e + t) p(N +1)η u(t, x) ≥ B A (x -x 0 ) log(e + t) (p-2)(N +1)η (e + t) η (3.10) for (t, x) ∈ (0, ∞) × R N . A first consequence of (3.10) with x = x 0 is the lower bound (3.5). We also infer from (3.10) with x = 0 that (e + t) N η log(e + t) p(N +1)η u(t, 0) ≥ B A (x 0 ) log(e + t) (p-2)(N +1)η (e + t) η and the right-hand side of the above inequality is positive provided t is large enough. Therefore there is t 0 > 0 such that u(t 0 , 0) > 0, and we argue as in the proof of (3.8) to complete the proof of (3.6), possibly taking a lower value of A if necessary.

4 Convergence. Proof of Theorem 1.1

Thanks to the outcome of Sections 2 and 3 we are in a position to prove Theorem 1.1.

To this end, we follow the lines of the analysis developed by Galaktionov & Vázquez in [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF][START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF], the central tool being a stability theorem which is recalled in Section A for the reader's convenience.

We fix an initial condition u 0 satisfying (1.3) and denote the corresponding solution to (1.1)-(1.2) by u. By Proposition 3.3, there is t 0 > 0 and A 0 ∈ (0, A sub ) such that

u(t 0 , x) ≥ B A 0 (x) , x ∈ R N . (4.1)
We then define w by (3.1)-(3.2) with u(• + t 0 ) instead of u, that is, u(t + t 0 , x) = (e + t) -N η log(e + t) -pη(N +1) w(s, y) , the variables (s, y) being still given by (3.1). We infer from (4.1), Lemma 3.1, Lemma 3.2, and the comparison principle that, for all s ≥ 1,

w(s) 1 + w(s) ∞ + ∇w ∞ ≤ C 6 , (4.2)
and w(s, y) ≥ B A 0 (y) , y ∈ R N , and w(s, y) = 0 , |y| ≥ R 0 .

We define the set

X := z ∈ L 1 (R N ) ∩ BC(R N ) : z(y) ≥ B A 0 (y) , y ∈ R N and z 1 ≤ C 6 ,
which is a complete metric space for the distance induced by the L ∞ -norm, the parameters C 6 and A 0 being given in (4.2) and (4.3), respectively. We also set S := {w} , and deduce from (3.1), (4.2), (4.3), and the properties of u that w ∈ C([0, ∞); X). We now check that the set S enjoys the three properties (H1)-(H3) required to apply the stability result from [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF][START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF] recalled in Theorem A.1 below. In our setting, the non-autonomous operator is the operator L defined in (3.3c) and its autonomous counterpart L is defined in (3.4b), the associated evolution equations being (3.3a) and (3.4a), respectively. Clearly, w is a solution to (3.3a) and it readily follows from the W 1,∞ -bound (4.2), the uniform localization of the support (4.3), and the Arzelá-Ascoli theorem that {w(s)} s≥0 is compact in L ∞ (R N ), so that (H1) is satisfied. We next infer from the same properties (4.2) and (4.3) that, for (s, y)

∈ (0, ∞) × R N , |Lw(s, y) -Lw(s, y)| ≤ C s [ ∇w(s) q * ∞ + w(s) ∞ + |y| ∇w(s) ∞ ] ≤ C s [C q * 6 + C 6 + R 0 C 6 ] .
Therefore,

Lw(s) -Lw(s) ∞ ≤ C 7 s , s > 0 ,
from which (H2) follows.

Finally, we fix

A 1 ∈ (A 0 , ∞) such that A 1 > B 0 R p/(p-1) 0 and B A 1 1 ≥ C 6 and set Ω := {B A : A ∈ [A 0 , A 1 ]} .
Clearly Ω is a non-empty and compact subset of X. Since (H3) concerns only the autonomous operator (3.4a), which also arises from the standard p-Laplacian equation via a self-similar change of variable, the uniform stability of Ω holds true as a by-product of classical results for the p-Laplacian equation (see [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF]Section 6] and also [START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF]Section 4.6]).

We have thus checked the validity of (H1)-(H3) and may thus apply the stability theorem (Theorem A.1) to conclude that the ω-limit set (for the topology of the uniform convergence) of w is included in Ω, that is, ω(w) := w ∈ X : there is a sequence (s j ) j≥1 , s j → ∞ , such that lim j→∞ w(s j ) -w ∞ = 0 ⊂ Ω . (4.4)

Mass analysis. Uniqueness of the limit. It remains to show that the asymptotic profile is in fact a uniquely determined Barenblatt profile with the parameter A * as in Theorem 1.1. To this end, we perform a mass analysis along the lines of a similar argument in [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF]. We first observe that classical properties of the Barenblatt profiles ensure that, given s ≥ 0, there is a unique A(s) > 0 such that

Θ(s) := w(s) 1 = B A(s) 1 . (4.5)
More precisely,

Θ(s) = N N + 1 ω N I 1 A(s) (p-1)/[p(p-2)η] , (4.6)
where In addition, integrating (3.3a) with respect to space shows that, for s ≥ 0,

I 1 := (N + 1) ∞ 0 1 -B 0 r p/(p-1) (p-1)/(p-2) + r N -1 dr (4.
dΘ ds (s) = G(w(s)) s with G(z) := (N + 1) z 1 -∇z q * q * . (4.9)
Owing to (4.9), (4.6), and the regularity properties of w, we realize that

A ∈ C([0, ∞)) ∩ C 1 ((0, ∞)). We next claim that lim s→∞ w(s) -B A(s) ∞ = 0 . (4.10)
Indeed, assume for contradiction that there are an increasing sequence (s j ) j≥1 of positive real numbers, s j → ∞, and ε > 0 such that

w(s j ) -B A(s j ) ∞ ≥ ε , j ≥ 1 . (4.11)
On the one hand, we infer from (4.2), (4.3), (4.4), and the Arzelá-Ascoli theorem that there are a subsequence of (s j ) j≥1 (not relabeled) and Ā ∈ [A 0 , A 1 ] such that lim j→∞ w(s j ) -B Ā ∞ = 0 .

On the other hand, it follows from (4.8) that, after possibly extracting a further susbequence, we may assume that there is

A ∞ ∈ [A 0 , A 1 ] such that A(s j ) → A ∞ as j → ∞.
This readily implies that lim j→∞ B A(s j ) -B A∞ ∞ = 0 .

We may then let j → ∞ in (4.11) to conclude that

B Ā -B A∞ ∞ ≥ ε . (4.12) Now, B Ā 1 = lim j→∞ Θ(s j ) = lim j→∞ B A(s j ) 1 = B A∞ 1 ,
so that Ā = A ∞ , which contradicts (4.12). We have thus proved (4.10). We next infer from (4.2) and (4.3) that, introducing f (s, y) := ηy • ∇w(s, y) + ηN w(s, y) where I 1 is defined in (4.7), and

- 1 s [|∇w(s, y)| q * -pη(N + 1)w(s, y) + (p -2)η(N + 1)y • ∇w(s, y)] for (s, y) ∈ (0, ∞) × R N , Eq. (3.3a) reads ∂ s w -∆ p w = f in (1, ∞) × B(0, R 0 + 1) with w ∈ L ∞ (1, ∞; W 1,∞ (B(0, R 0 + 1))) and f ∈ L ∞ ((1, ∞) × B(0, R 0 + 1)).
I 2 := η q/(p-1)
∞ 0 1 -B 0 r p/(p-1) q/(p-2) + r ((p-1)(N -1)+q)/(p-1) dr .

We observe that g vanishes only once in (0, ∞), for a = A * , which is defined in (1.5) and reads

A * := I 1 I 2 p(p-2)(N +1)η/(p-1)
, with the notation introduced in this section. In fact, g(a) < 0 for a > A * and g(a) > 0 for a ∈ (0, A * ) . We finally note that Proposition 1.2 follows at once from (4.3) after translating these properties in terms of u.

A The stability theorem

We briefly recall here for the reader's convenience the stability theorem introduced by Galaktionov and Vázquez in [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF][START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF] and used in Section 4. As a general framework, consider a non-autonomous evolution equation

∂ s ϑ = Lϑ , (A.1)
that can be seen as a small perturbation of an autonomous evolution equation with good asymptotic properties ∂ s Φ = LΦ , (A.2) in the sense described below. We consider a set S of solutions ϑ ∈ C([0, ∞); X) to (A.1) with values in a complete metric space (X, d). We assume that:

(H1) For each ϑ ∈ S, the orbit {ϑ(t)} t>0 is relatively compact in X. Moreover, if we let ϑ τ (t) := ϑ(t + τ ) , t ≥ 0 , τ > 0, then {ϑ τ } τ >0 is relatively compact in C([0, T ]; X) for any T > 0.

(H2) Let ϑ ∈ S for which there is a sequence of positive times (t k ) k≥1 , t k → ∞, such that ϑ(• + t k ) -→ Θ in C([0, T ]; X) as k → ∞ for any T > 0. Then Θ is a solution to (A.2).

(H3) Define the ω-limit set Ω of (A.2) in X as the set of f ∈ X enjoying the following property: there are a solution Φ ∈ C([0, ∞); X) to (A.2) and a sequence of positive times (t k ) k≥1 such that t k → ∞ and Φ(t k ) -→ f in X. Then Ω is non-empty, compact and uniformly stable, that is: for any ε > 0, there exists δ > 0 such that if Φ is a solution to (A.2) with d(Φ(0), Ω) ≤ δ, then d(Φ(t), Ω) ≤ ε for any t > 0.

The stability theorem (also known as the S-theorem) then reads:

Theorem A.1. If (H1)-(H3) above are satisfied, then the ω-limit set of any solution ϑ ∈ S is contained in Ω.

For a detailed proof we refer the reader to [START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF][START_REF] Galaktionov | A Stability Technique for Evolution Partial Differential Equations. A Dynamical Systems Approach[END_REF].

7 )

 7 and ω N denotes the volume of the unit ball of R N , see[START_REF] Vázquez | Smoothing and Decay Estimates for Nonlinear Diffusion Equations[END_REF] Section 11.4.1] for instance. Since B A 1 1 ≥ C 6 ≥ Θ(s) ≥ B A 0 1 , s ≥ 0 , by (4.2), (4.3), and the choice of A 1 , we deduce from the monotonicity of A → B A 1 and (4.5) that A 0 ≤ A(s) ≤ A 1 , s ≥ 0 . (4.8)

  We then infer from[START_REF] Dibenedetto | Hölder estimates for nonlinear degenerate parabolic systems[END_REF] Theorem 1.1] that there are C 8 > 0 and ν ∈ (0, 1) such that|∇w(s 1 , y 1 ) -∇w(s 2 , y 2 )| ≤ C 8 |y 1 -y 2 | ν + |s 1 -s 2 | ν/2for all s 2 ≥ s 1 ≥ 2 and (y 1 , y 2 ) ∈ B(0, R 0 ) × B(0, R 0 ). Combining the above property with (4.2) we deduce that {∇w(s)} s≥1 is bounded and equicontinuous in C(B(0, R 0 ); R N ) and thus compact in that space by the Arzelá-Ascoli theorem. Recalling (4.3) and (4.10) we conclude thatlim s→∞ w(s) -B A(s) r + ∇w(s) -∇B A(s) r = 0 , r ∈ [1, ∞] . (4.13)An immediate consequence of (4.13) is thatlim s→∞ G(w(s)) -G(B A(s) ) = 0 ,(4.14)while the explicit formula for B A(s) givesG(B A(s) ) = g(A(s)) (4.15)with g(a) := N ω N I 1 a (p-1)/ηp(p-2) -N ω N I 2 a (N +2)(p-1)/(N +1)ηp(p-2) , a > 0 ,

(4. 16 )

 16 Thanks to (4.4), (4.9), (4.13), (4.14),(4.15), and (4.16), we are in a position to proceed as in the proof of[START_REF] Galaktionov | Asymptotic behaviour of nonlinear parabolic equations with critical exponents. A dynamical systems approach[END_REF] Proposition 5.2] to establish that ω(w) = {B A * }. Undoing the rescaling (3.1)-(3.2), this property readily gives Theorem 1.1.