This article is mainly concerned how the double exponential formula for numerical integration was discovered and how it has been developed thereafter. For evalu1 ation of I = ∫1−1 f(x)dx H. Takahasi and M. Mori took advantage of the optimality of the trapezoidal formula over (−∞, ∞) with an equal mesh size and proposed in 1974 x = φ(t) = tanh(π/2 sinh t) as an optimal choice of variable transformation which 2 ∞ transforms the original integral to I = ∫∞−∞ f(φ(t))φ' (t) dt. If the trapezoidal formula is applied to the transformed integral and the resulting inﬁnite summation is properly truncated a quadrature formula I ≈ h ∑_n__k_=−_n_ f(φ(kh))φ'(kh) is obtained. Its error is expressed as O(exp(−_CN_/ log N )) as a function of the number N ( = 2_n_ + 1) of function evaluations. Since the integrand decays double exponentially after the transformation it is called the double exponential (DE) formula. It is also shown that the formula is optimal in the sense that there is no quadrature formula obtained by variable transformation whose error decays faster than O(exp(−_CN_/ log N )) as N becomes large. Since the paper by Takahasi and Mori was published the DE formula has gradually come to be used widely in various ﬁelds of science and engineering. In fact we can ﬁnd papers in which the DE formula is successfully used in the ﬁelds of molecular physics, ﬂuid dynamics, statistics, civil engineering, ﬁnancial engineering, in particular in the ﬁeld of the boundary element method, and so on. The DE transformation has turned out to be also useful for evaluation of indeﬁnite integrals, for solution of integral equations and for solution of ordinary diﬀerential equations, so that the scope of its applications is expected to spread also in the future.