1 Introduction
In this paper, we consider the computation of an integral
(1) 
where is a slowly decaying oscillatory function, such as
It is difficult to compute this integral by conventional methods such as the DE formula [8]. In this paper, we propose a numerical method for computing integrals of this type. In the proposed method, we consider an integral of the FourierLaplace transform type
(2) 
which is a complex analytic function in the upperhalf complex plane under some conditions on . We obtain the function in the upperhalf plane, and we compute the analytic continuation of on the real axis . Then, we obtain the desired integral by
(3) 
Our method is related to hyperfunction theory, a theory of generalized functions proposed by Sato [7]. Roughly speaking, a hyperfunction is given by the boundary values of functions which are analytic in subdomains of . From (3), the desired integral is a boundary value on of the function analytic in the upper half plane , and, in this sense, the integral can be regarded as a value of a hyperfunction.
Previous studies related to this paper are as follows. Toda and Ono proposed a method of computing (1) by the limit
which is obtained by the DE formula and the Richardson extrapolation [9]. Sugihara improved Toda and Ono’s method and proposed a method of computing (1) by the limit
which is obtained by the DE formula and the Richardson extrapolation. Ooura and Mori proposed a DEtype formula designed especially for Fourier transform type integrals using a unique DE transform
[5] and improved it for integrands with singularities in the complex plane [6]. The author proposed a method of computing Fourier transform in a way similar to that of this paper [4].The remainder of this paper is structured as follows. In Section 2, we propose a numerical method for computing the oscillatory integral (1). In Section 3, we give a brief review of hyperfunction theory and mention a relation between the proposed method and hyperfunction theory. In Section 5, we give the summary and concluding remarks of this paper and refer to problems for future study.
2 Numerical method
We consider an integral (1), where is a slowly decaying oscillatory function and propose a numerical method for computing it. As mentioned in the previous section,

we obtain a function , which is a complex analytic function in the upper half plane under some conditions on ,

and we obtain the desired integral by the analytic continuation of onto the real axis.
The details of each step are as follows.

(Computation of in the upper half plane ) We obtain the analytic function as a Taylor series
where is a complex constant given by the user such that . The Taylor coefficient is given by the integral
(4) This integral involves an exponentially decaying factor , where we remark that , and it is easy to compute it by conventional numerical integration formulas.

(Analytic continuation of ) We obtain the analytic continuation of onto the real axis using its continued fraction expansion since the continued fraction expansion of an analytic function has a larger convergence region compared with that of its Taylor series. We can convert a Taylor series
into a continued fraction
(5) using the socalled quotient difference (QD) algorithm (see Theorem 12.4c of [3]). Namely, we generate the sequences and by the recurrence relations
and, then, we obtain
Once we obtain the coefficients of the continued fraction (5), we can compute the continued fraction by
where and are the polynomials given by the recurrence relations
In the procedure of the QD algorithm, we use multiple precision arithmetic since the algorithm is numerically unstable.
Finally, we obtain the desired oscillatory integral by
(6) 
3 Relation with hyperfunction theory
In this section, we give a brief review of hyperfunction theory and its relation with the proposed method. For the detail of hyperfunction theory, see the textbook [2].
Intuitively, a hyperfunction is defined by the difference between the boundary values of a complex analytic functions in a subdomain of the upper half plane and a complex analytic function in a subdomain of the lower half plane
where the analytic functions are called defining functions of the hyperfunction . We also use the notation
where
For example, the Dirac delta function is defined as a hyperfunction by
and the Heaviside step function is defined as a hyperfunction by
where the complex logarithmic function is the principal value, that is, the branch such that it takes a real value on the positive real axis. Figure 1 shows the graphs of the defining functions of the Dirac delta function and the Heaviside step function . From these figures, we find easily that the difference between the boundary values on the real axis of the defining functions gives the hyperfunctions.
x
y
z

x
y
z

More precisely, hyperfunctions are defined as follows. Let be an open interval on the real axis and be a complex neighborhood of , that is, a complex domain including as a closed subset. We denote the set of complex analytic functions in by , which can be regarded as a complex linear space by the usual summation and multiplication by complex numbers, and the set of complex analytic functions in by , which can be regarded as a linear subspace of . Then, we can consider the quotient space . We call an element of , that is, an equivalence class a hyperfunction on and denote it by
where we call a defining function of the hyperfunction . We also use the notation
where
In other words, a hyperfunction on is a complex analytic function in , where another complex analytic function in is identified with if is analytic in .
We define the value of a hyperfunction at a point by
if the limit on the right hand side exists. We do not define the value if the limit does not exist.
4 Numerical examples
In this section, we discuss show some numerical examples that illustrate the effectiveness of the proposed method. All the computations were performed using programs coded in C++ with 100 decimal digit precision, where we use the multiple precision arithmetic library exflib [1].
We computed the oscillatory integrals
by the proposed method. We took the center of the Taylor series as , and we computed the Taylor coefficients given by the integral (4), which were computed using the DE formula. Table 1 show the relative error and the number of functional evaluations for in (1) of the proposed method. Table 1 show that we can compute oscillatory integrals with high accuracy by the proposed method.
integral  (1)  (2)  (3)  (4) 
relative error  5.4E26  6.2E35  3.8E34  1.4E36 
number of functional evaluations  917  964  957  927 
integral  (5)  (6)  (7)  (8) 
relative error  1.3E35  3.8E36  1.1E33  2.1E37 
number of functional evaluations  954  958  927  947 
Here, we point out an interesting fact. Since we took , the integral in (4) which gives the Taylor coefficients is given by
which does not involve any oscillatory function. It means that we can obtain oscillatory integrals without computing oscillatory integrals. In other words, we can obtain Fourier transform type integrals by computing Laplace transform type integrals.
As a comparison with the results by a previous method, we computed some of the integrals by the Euler transform for computing alternating series. Namely, we express the desired oscillatory integral as
(7) 
where and
We computed each integral on the right hand side of (7) by the 100point GaussLegendre formula, and we obtained the desired integral by applying the Euler transform to the alternating series (7). Table 2 show the relative errors of the method using the Euler transform applied to the integrals (2,3,4,5). The number of functional evaluations for in (1) are 5000 in all the computations of the integrals. Comparing the results in Table 1 and 2, the proposed method is far superior to the method using the Euler transform.
integral  (2)  (3)  (4)  (5) 
relative error  6.3E25  1.9E25  3.1E25  1.4E25 
5 Summary
In this paper, we proposed a numerical method of computing an integral involving a slowly decaying oscillatory function. In the proposed method, we introduce a complex analytic function in the upper half plane and obtain the desired integral by the analytic continuation of this analytic function onto the real axis. The analytic function is obtained as a Taylor series by computing its coefficients, which are given by integrals involving exponentially decaying functions and easily evaluated by the conventional numerical integration formula. The analytic continuation is performed by converting the Taylor series into a continued fraction using the quotient difference (QD) algorithm. Then, the desired integral is obtained as the value of the analytic function at the origin. We also remarked that the proposed method is closely related to hyperfunction theory in the sense that the desired oscillatory integral can be regarded as the boundary value of the hyperfunction whose defining function is the above complex analytic function.
As presented in this paper, hyperfunctions, that is, the boundary values of analytic functions appear in many studies of science and engineering. In the evaluation of hyperfunctions, the techniques of analytic continuation are necessary. Therefore, the author believes that the study and the computation of hyperfunctions together with analytic continuation will be important and useful in applied mathematics.
Problems for future study related to this paper are, for example, theoretical error estimate of the proposed method and the reduction of the computational cost of analytic continuation, which needs multiple precision arithmetic in the QD algorithm applied to the conversion of a Taylor series into a continued fraction.
References
 [1] Exflib information. Note: http://wwwan.acs.i.kyotou.ac.jp/~fujiwara/exflib/ Cited by: §4.
 [2] (2010) Introduction to hyperfunctions and their integral transforms — an applied and computational approach. Birkha̋user, Basel. Cited by: §3.
 [3] (1977) Applied and computational complex analysis. Vol. 2, John Wiley & Sons, New York. Cited by: item 2.
 [4] (2018) A numerical method of fourier transform based on hyperfunction theory. Note: arXiv:1808.03460v1 [math.NA] Cited by: §1.
 [5] (1991) The double exponential formula for oscilattory functions over the half infinite interval. J. Comput. Appl. Math. 38, pp. 353–360. Cited by: §1.
 [6] (1999) A robust double exponential formula for fourier type integrals. J. Comput. Appl. Math. 112, pp. 229–241. Cited by: §1.
 [7] (1959) THeory of hyperfunctions. J. Fac. Sci. Univ. Tokyo, Sect. 1A Math 8, pp. 139–193. Cited by: §1.
 [8] (1978) Double exponential formulas for numerical integration. Publ. Res. Inst. Math. Sci., Kyoto Univ. 339, pp. 721–741. Cited by: §1.
 [9] (1978) Some remarks for efficient usage of the double exponential formulas (in japanese). Kokyuroku, Res. Inst. Math. Sci. Kyoto Univ. 339, pp. 74–109. Cited by: §1.