A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
arXiv:2607.23196v1 [math.CA] 25 Jul 2026
FRANÇOIS ALOUGES, GIOVANNI DI FRATTA, ALBERTO FIORENZA, AND RENATO FIORENZA Abstract. We present a unified fixed-point construction of the elementary transcendental functions, encompassing the real exponential, the complex exponential (sine and cosine), and the natural logarithm. Each function is characterized as the unique solution of a duplication identity established through the Banach contraction principle. These foundational identities are e(2x) = e2 (x) for the exponentials, and log(x2 ) = 2 log x for the logarithm. Because a direct iteration of these identities is numerically unstable, owing to local expansiveness at the target, the central idea transfers the analysis to a residual function, on which the operator becomes a strict contraction with an explicit convergence rate. Beyond its theoretical economy, which dispenses with differential equations and power series, this characterization translates into efficient algorithms for the machine evaluation of elementary functions: the underlying framework yields floating-point kernels whose accuracy and iteration depth are governed by the theoretical contraction rate. We also present a computational study showing that, in a throughput-bound vectorized regime, these kernels are competitive with standard production libraries, and in favorable configurations exceed them, with the sine–cosine kernel faster at every tested iteration depth. These implementations operate without lookup tables or memory traffic, an architectural advantage for modern high-performance and energy-efficient computing.
1. Introduction The elementary transcendental functions, namely the exponential, the logarithm, and the trigonometric functions, are among the most frequently evaluated primitives in scientific computing. Large-scale numerical simulations routinely require billions of evaluations and the routines that serve those calls, refined over decades in libraries such as the system libm, Intel’s SVML, and the vectorised library SLEEF [SP20], almost invariably rest on the same architectural choice: a precomputed lookup table brings the argument into a narrow interval, after which a short polynomial finishes the evaluation. Lookup tables are efficient in scalar code but poorly suited to modern Single Instruction Multiple Data (SIMD) architectures. SIMD execution achieves high throughput by evaluating multiple independent arguments simultaneously, but different vector lanes generally require different table entries. The resulting irregular memory accesses are implemented through gather operations, which incur substantially higher latency than contiguous loads and increase memory traffic. As a result, arithmetic pipelines are often stalled while waiting for data, reducing the throughput that SIMD parallelism is intended to deliver. In this work, we develop a table-free construction of the elementary transcendental functions based on a structural characterisation rather than on polynomial approximation. Our starting point is the observation that each function is uniquely determined by a duplication identity, namely a functional equation relating the value at an argument to the value at either half or twice that argument, together with a natural normalisation condition. For the real and complex exponentials, this characterisation is captured by the following axioms: (A1 ) Duplication identity: e(2x) = e(x)2
for all x ∈ R,
(1)
2020 Mathematics Subject Classification. 39B22, 47H10, 65D15, 65D20, 33B10, 26A09, 68W40. Key words and phrases. Elementary transcendental functions, exponential function, logarithm, fixed-point iteration, Banach contraction principle, functional equations, duplication identity, floating-point arithmetic, table-free evaluation, special function computation, high-performance computing, SIMD vectorization.
1
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
2
(A2 ) Initial conditions: e(0) = 1,
e′ (0) = γ with γ ∈ {1, i}.
(2)
The parameter γ explicitly selects either the real exponential ex or the imaginary exponential eix , and through the latter, the sine and cosine functions. For the natural logarithm, the analogous structural definition imposes: (B1 ) Functional equation: f(x2 ) = 2f(x)
for x > 0,
(3)
f(1) = 0,
f′ (1) = 1.
(4)
(B2 ) Cauchy initial conditions:
The common principle unifying these definitions is that the target function is characterised as the unique fixed point of an operator induced directly by the respective duplication identity. Once this operator is shown to be contractive, the Banach–Caccioppoli theorem immediately yields existence, uniqueness, and a convergent iterative construction on arbitrary compact intervals. The route to a contraction is the one subtle point, and it is the same for all three functions. Applied naively, the duplication map is unstable: the linearisation of the squaring z 7→ z 2 at the initial value e(0) = 1 has derivative 2, so the operator is locally expansive and the fixed point repelling. We overcome this by transferring the analysis to a residual function. Writing e(x) = 1 + γx + x2 q(x) isolates the second-order remainder q, and the duplication identity, rewritten for q, becomes a genuine contraction once the ambient space is equipped with a suitable weighted norm. The logarithm admits an analogous reformulation based on repeated square-root reduction rather than repeated squaring. Although the two constructions are mathematically dual, Section 6 shows that this distinction has important consequences in floating-point arithmetic. The main contribution of this work is that the above characterization is not merely theoretical, but intrinsically constructive. The contraction argument underlying the existence theorem directly yields an algorithmic procedure: its two parameters, namely the degree of the polynomial seed and the number of iteration steps, are exchangeable at the explicit rate determined by the contraction. The argument reduction relies exclusively on arithmetic operations, including exponent extraction and square roots, rather than lookup tables. Consequently, the resulting computational kernels operate without memory accesses beyond register-resident data. Section 6 substantiates these claims through a detailed floating-point analysis. On throughput-oriented vector hardware, where arithmetic latency is amortized across SIMD lanes and the elimination of table accesses becomes advantageous, the fixed-point kernels are competitive with a production vector library, and in favorable configurations outperform it, most consistently for the sine and cosine. Moreover, the iterative formulation naturally exposes an accuracy-performance trade-off: the target precision can be deliberately relaxed to reduce the computational cost and achieve higher throughput when full machine accuracy is not required. By systematically replacing memory traffic with arithmetic operations, the proposed kernels also offer the potential for improved energy efficiency, a direction left for future experimental investigation. Organization of the paper. The remainder of the paper is organized as follows. Section 2 contextualizes our approach within the classical methodologies for defining elementary functions, highlighting the conceptual shift from binary functional equations to unary duplication identities. Section 3 analyzes the fixed point operator governing the real and complex exponentials, establishing the main existence and uniqueness results derived from the Banach–Caccioppoli contraction principle. Section 4 supplies the rigorous proofs for the foundational lemmas of this operator by introducing the concept of admissible weight sequences to strictly control its nonlinear behavior. Section 5 carries out the parallel construction for the logarithm, defining its associated linear operator and proving its global strict contraction within an appropriately constructed complete metric space. Finally, Section 6 details the practical realization of these theoretical operators, translating them into table-free floating-point kernels and benchmarking their performance in a vectorized regime against the production library SLEEF.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
3
2. Ways of defining the elementary functions The definition of an elementary function through a functional relation, rather than an explicit formula, is a recurring theme in analysis. Before introducing our own characterization, we briefly review the classical alternatives, both to place our approach in context and because each relies on a prerequisite that our method avoids. From a computational perspective, the numerical realization of the approaches reviewed here has been extensively studied and is now well understood; see, for example, Muller et al. [MBdD+ 18]. Accordingly, we take this opportunity to adopt a more pedagogical viewpoint, focusing on the mathematical structure underlying the different characterizations before presenting the computational framework developed later in this work. Power series: One may define ez , sin x, and cos x directly as the sums of their power series [Rud76, Lan97, Ahl79]. Although algebraically clean, this route demands a comprehensive understanding of infinite series merely to grasp the definitions. By contrast, characterizing these functions through natural duplication identities offers a distinct pedagogical advantage. The defining statements are elegant and entirely comprehensible without any prior knowledge of sequences or series. From a didactical perspective, one can introduce the functions through the natural properties they satisfy and assert their unique existence early on, deferring the proof of the underlying fixed point theorem to a later stage. Consequently, the foundational definitions remain highly accessible and require only minimal prerequisites in analysis. Integral functions and ODEs: Transcendental functions may be defined either as integral functions [Spi08, Har08] or as the unique solutions to initial value problems. For instance, one commonly defines Z x Z x dt dt √ log x = , x ∈ (0, +∞), arcsin x = , x ∈ [−1, 1], t 1 − t2 1 0 with the exponential and trigonometric functions subsequently recovered through functional inversion. Alternatively, the sine and cosine functions are frequently introduced as the unique solutions to the second order linear ordinary differential equation y ′′ +y = 0, subject to the respective initial conditions y(0) = 0, y ′ (0) = 1 and y(0) = 1, y ′ (0) = 0. Similar pedagogical remarks apply to these methods. Axiomatic approaches: To permit an early development of calculus, several authors introduce the elementary functions by postulating their defining properties and deriving the calculus from these, deferring the constructive existence proof to a later chapter. Apostol, in his classic text [Apo67, §2.5], characterises the sine and cosine not by differential equations but by their fundamental properties. He postulates functions cos and sin on R with cos 0 = 1, cos π = −1, and sin(π/2) = 1, satisfying the difference formula cos(y − x) = cos y cos x + sin y sin x
(x, y ∈ R),
together with the local inequality 0 < cos x < sinx x < cos1 x for 0 < x < π/2. From these, he derives their derivatives immediately, while postponing the proof that such functions exist to the chapter on infinite series. Shilov proceeds in the same spirit [Shi96, Ch. 5], defining sin and cos as functions satisfying the Pythagorean identity, the addition formulas sin(x + y) = sin x cos y + cos x sin y,
cos(x + y) = cos x cos y − sin x sin y
(x, y ∈ R),
sin x and the fundamental inequality 0 < sin x < x < cos x for sufficiently small x > 0, again deferring existence to the later theory of power series. For the remaining two functions, the two authors part Rx ways. Apostol defines the logarithm by integral inversion, log x = 1 dt/t. Shilov, by contrast, characterises the logarithm axiomatically as the unique increasing function on (0, +∞) satisfying the homomorphism law ℓ(xy) = ℓ(x) + ℓ(y) together with a normalisation ℓ(a) = 1 for a given a > 1, and subsequently obtains the exponential as its inverse function. The axiomatic characterisations just recalled share a structural feature that, on reflection, one might have thought unavoidable: each defines its function through a relation between its values at a pair of independent points. The addition and difference formulas for the trigonometric functions relate
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
4
cos(x ± y) to the values at x and at y. The exponential, in the standard treatment, is fixed by the Cauchy relation e(x + y) = e(x)e(y), expressing that it is a homomorphism of the additive group into the multiplicative one; the logarithm, dually, is fixed by ℓ(xy) = ℓ(x) + ℓ(y). Which of the last two is taken as primitive and which as the inverse of the other is a matter of expository taste. Shilov starts from the logarithm, whereas texts such as Prodi’s [Pro77] start from the exponential. However, in either case, the defining equation is binary: it constrains the function on R × R, and it is exactly this binary, homomorphism-theoretic form that has long been taken to be the essence of an axiomatic definition. It is also what obstructs a fixed-point reading, because a relation in two independent arguments is not an operator acting on functions of one variable, leaving the contraction principle with nothing to act upon. The starting point of the present paper is the observation that the binary relation is not necessary. The duplication identities e(2x) = e(x)2
and
log(x2 ) = 2 log x
(5)
are unary: each relates the value of the function at a single point to its value at the doubled or squared argument, imposing no constraint that involves a second, independent variable. Supplemented only by an elementary normalisation at the origin and a minimal regularity requirement, they suffice to determine the function uniquely. Consequently, the binary homomorphism law is recovered as a consequence, rather than assumed as a hypothesis. This reduction from a binary to a unary relation is precisely what makes the fixed-point formulation possible. A unary duplication identity naturally defines an operator acting on a space of single-variable functions, and it is exactly upon this operator that the Banach contraction principle can be applied. The remainder of the paper develops this observation and, in Section 6, details its computational consequences. Two further remarks locate our construction more sharply. First, the geometric definition of the trigonometric functions through arc length on the unit circle conceals a well-known circularity, noted by Richman [Ric84] and revisited recently [BCFG26]: arc length is used to define sine and cosine, yet those functions are later used to define arc length by integration [Ste16, SC17]. A fixed-point characterisation is free of this circularity because it never invokes arc length. Second, the duplication structure we exploit is closely related to the “scaling and squaring” method for the matrix exponential [AMH09, Hig05, Hig08, ML03], which is one of the most effective algorithms in numerical linear algebra. Furthermore, the study of the underlying functional equations is itself classical [Acz66, Kuc09, GD03, Zei86, Tre19]. Our contribution is to make the contraction global and explicit, achieving uniformity on arbitrary compact intervals through a weighted norm, and thereby turning the characterisation into a numerically stable iteration. 3. Analysis of the Fixed-Point Operator In this section, we prove the existence and uniqueness of a continuous function e : R → C satisfying the duplication and initial-condition axioms (A1 , A2 ), defined in (1) and (2). The key ingredient of the proof is the classical Banach–Caccioppoli contraction principle, which we recall for completeness (see [GD03, Zei86] for further details). Theorem 1 (Banach–Caccioppoli). Let (X, d) be a complete metric space, and let T : X → X be a contraction mapping; that is, there exists a constant κ ∈ [0, 1) such that d(T (x), T (y)) ⩽ κd(x, y)
for all x, y ∈ X.
(6)
Then T admits a unique fixed point: there exists a unique x∗ ∈ X such that T (x∗ ) = x∗ . Furthermore, for any arbitrary initial point x0 ∈ X, the iterative sequence defined by xn+1 = T (xn ) converges to x∗ in the metric d, and the following error estimates hold for all n ⩾ 1: (i)
A priori estimate: d(xn , x∗ ) ⩽
κn d(x1 , x0 ). 1−κ
(7)
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
(ii)
5
A posteriori estimate: d(xn , x∗ ) ⩽
κ d(xn , xn−1 ). 1−κ
(8)
Remark 1. The error estimates of Theorem 1 supply both theoretical and practical guarantees for the convergence behaviour of iterative schemes. The a priori estimate establishes a theoretical convergence rate: the error decays at least geometrically, with a factor κn , and thus quantifies how rapidly iterates approach the fixed point in the asymptotic regime. The a posteriori estimate, by contrast, is computable from the iterates and furnishes an on-line measure of the current error. In numerical implementations this distinction is fundamental: the a priori bound informs algorithm design and expectations about worst-case performance, while the a posteriori bound is indispensable for reliable error control, stopping criteria, and the detection of spurious stagnation. In particular, the latter guarantees that if the inter-iterate distance d(xn , xn−1 ) is small, then the iterate xn is provably close to the true solution x∗ , so that apparent stagnation cannot occur at a significant distance from the fixed point. Our argument relies on the following two lemmas, whose proofs (which utilize Theorem 1 operating on suitable weighted norms) are deferred to the next section. Lemma 1. Let γ ∈ {1, i}. There exists a unique continuous function q∗ ∈ C 0 (R, C) which solves the functional equation γ2 1 x x x2 x q(x) = + +γ q + q2 (x ∈ R) . (9) 4 2 4 2 16 2 Moreover, q∗ is the unique solution to (9) on any interval Ia := [−a, a], a > 0, within the class B(Ia , C) of bounded (not necessarily continuous) functions. Lemma 2. Let γ ∈ {1, i}. Let e : R → C be a continuous function satisfying the duplication and initial condition axioms (A1 , A2 ). Then, there exists a unique continuous function q : R → C, which we refer to as the residue, such that e(x) = 1 + γx + x2 q(x)
(x ∈ R).
(10)
Moreover, this function q coincides with the unique continuous solution of the functional equation (9). In particular, q(0) = γ 2 /2. Remark 2. We emphasize that the uniqueness of q must be understood within the class of continuous functions. Without this restriction, the value of q(0) could be chosen arbitrarily, since the factor x2 in (10) makes the representation insensitive to the value at the origin. Continuity therefore ensures that the residue is uniquely determined. Theorem 2. There exists a unique continuous function e : R → C satisfying the duplication and initial condition axioms (A1 , A2 ). Specifically, e(x) = 1 + γx + x2 q(x) with q being the unique continuous solution of the functional equation (9). Proof. We treat uniqueness and existence separately. Uniqueness. Suppose that e1 and e2 are two continuous functions satisfying both axioms (A1 , A2 ). By Lemma 2, there exist continuous residue functions q1 and q2 such that ej (x) = 1 + γx + x2 qj (x) for j ∈ {1, 2}, and both residues obey the functional equation (9). By Lemma 1, this functional equation has a unique continuous solution; hence q1 = q2 . It follows that e1 (x) = e2 (x) for all x ∈ R, proving uniqueness. Existence. Let q∗ denote the unique continuous solution of (9) provided by Lemma 1, and define e(x) := 1 + γx + x2 q∗ (x). We verify that e satisfies both axioms.
(11)
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
6
To check the duplication identity (A1 ), we compute the square of e(x): e2 (x) = (1 + γx + x2 q∗ (x))2 = 1 + γ 2 x2 + x4 q2∗ (x) + 2γx + 2x2 q∗ (x) + 2γx3 q∗ (x) 2 (2x)2 2 1 2x 2 γ = 1 + γ(2x) + (2x) q∗ (x) + + +γ q (x) . 4 2 4 16 ∗
(12)
Since q∗ satisfies the functional equation (9), the expression inside the square brackets is precisely the definition of q∗ (2x). Therefore, the equation simplifies to: e2 (x) = 1 + γ(2x) + (2x)2 q∗ (2x) = e(2x),
(13)
which proves that axiom (A1 ) holds. Finally, we verify the initial conditions (A2 ). Direct evaluation yields e(0) = 1. For the derivative at the origin, we compute the limit of the difference quotient: e(x) − e(0) 1 + γx + x2 q∗ (x) − 1 = lim = lim (γ + xq∗ (x)). (14) x→0 x→0 x→0 x x Since q∗ is bounded (and continuous) at 0, we have limx→0 xq∗ (x) = 0; thus e′ (0) = γ. This confirms (A2 ) and completes the proof. □ e′ (0) = lim
4. Proofs of Lemma 1 and 2 The existence and uniqueness of the fixed point on an arbitrary compact interval Ia := [−a, a] with a > 0, are established through the construction of a weighted Banach space. This setting enables effective control of the contraction induced by the linear component of the operator, while simultaneously suppressing the quadratic nonlinearity and ensuring a uniform bound on the parameter γ. 4.1. Admissible weight sequences. A central difficulty in proving existence and uniqueness lies in handling the nonlinearity of the operator T . When C 0 is equipped with the standard uniform norm, the quadratic term prevents T from being a global contraction. To address this issue, we introduce the notion of admissible weight sequences. By working within a suitably weighted Banach space, one can control the linear contraction of the domain and suppress the quadratic growth by allowing a tunable parameter λ to tend to infinity. In this subsection, we first formulate a set of abstract axioms that characterize these admissible weight sequences. Next, we demonstrate the practical applicability of this framework by explicitly verifying the conditions for two concrete families of functions: Gaussian weights and polynomial power weights. Definition 1 (Admissible weight sequences). Let Ω = (ωλ )λ∈N be a sequence of continuous functions ωλ : R → R. We say Ω is an admissible weight sequence if, for any compact interval Ia , it satisfies the following properties eventually in λ, in the sense that for every compact interval Ia there exists λ∗ = λ∗ (a) ∈ N such that for all λ ⩾ λ∗ the following conditions hold: (i)
Strict positivity and lower bound: inf ωλ (x) ⩾ 1.
x∈Ia
(ii)
(15)
Linear suppression: Defining the ratio ωλ (x/2) , ωλ (x)
(16)
lim sup [|x|R1 (x, λ)] = 0.
(17)
R1 (x, λ) := we have sup R1 (x, λ) ⩽ 1 x∈Ia
and
λ→∞ x∈Ia
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
7
(iii) Quadratic suppression: Defining the quadratic ratio R2 (x, λ) :=
ωλ2 (x/2) , ωλ (x)
(18)
we require lim sup [x2 R2 (x, λ)] = 0.
λ→∞ x∈Ia
(19)
We emphasize that all estimates are uniform on Ia . However, the threshold λ∗ is allowed to depend on the interval size a. Note that all estimates are uniform on Ia , but the threshold λ∗ and the associated constants may depend on a. Let Xa = C 0 (Ia , C), ∥· ∥λ be the normed space of continuous functions equipped with the weighted norm (note that Xa also depends on λ, though we suppress this in the notation for brevity): |q(x)| ∥q∥λ := sup . (20) x∈Ia ωλ (x) Since ωλ is continuous and strictly positive on the compact set Ia , there exist constants 0 < mλ ⩽ Mλ < ∞ such that mλ ⩽ ωλ (x) ⩽ Mλ for all x ∈ Ia . It follows that 1 1 ∥f ∥∞ ⩽ ∥f ∥λ ⩽ ∥f ∥∞ , (21) Mλ mλ so the weighted norm ∥ · ∥λ is equivalent to the standard uniform norm on C 0 (Ia , C), and hence Xa is a Banach space. Let Bλ,1 = {q ∈ Xa : ∥q∥λ ⩽ 1} denote the closed unit ball. Endowed with the metric induced by ∥ · ∥λ , the set Bλ,1 is a complete metric space, being a closed subset of a Banach space. Our goal is to show that, for sufficiently large λ, the operator T maps Bλ,1 into itself and is a strict contraction on Bλ,1 . The Banach–Caccioppoli fixed-point theorem (Theorem 1) will then guarantee the existence and uniqueness of a fixed point of T in Bλ,1 . Before that, to illustrate the effectiveness of this framework, we now verify the admissibility conditions for two concrete classes of weights: Gaussian weights and polynomial power weights. In both cases, the key estimates rely on the following elementary bounds: for any c > 0 and λ ∈ N, the 2 2 functions |x|e−cλx and x2 e−cλx attain their global maxima on R at x2 = 1/(2cλ) and x2 = 1/(cλ), respectively. Consequently, we have the uniform global bounds: 1 1 2 2 and sup[x2 e−cλx ] = sup[|x|e−cλx ] = √ (22) ceλ 2ceλ x∈R x∈R In particular, both quantities vanish as λ → ∞. Proposition 1 (Gaussian Weights). The sequence of Gaussian weights 2
ωλ (x) = eλx ,
λ ∈ N,
(23)
is an admissible weight sequence. Proof. Let Ia = [−a, a] be an arbitrary compact interval. We verify the three conditions of Definition 1. 2 The lower bound (15) is clear: Since λ ⩾ 0 we trivially have ωλ (x) = eλx ⩾ e0 = 1 for all x ∈ Ia . 2 For the linear suppression, we observe that R1 (x, λ) = e−(3/4)λx ⩽ 1. Hence, applying the bound from (22) with c = 3/4, we obtain: 1 λ→∞ sup [|x|R1 (x, λ)] ⩽ p −−−→ 0. (3/2)eλ x∈Ia It remains to prove the quadratic suppression condition (19). The quadratic ratio is R2 (x, λ) = 2 e−(1/2)λx . Applying (22) with c = 1/2: 2 λ→∞ sup [x2 R2 (x, λ)] ⩽ −−−→ 0. eλ x∈Ia
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
8
□
Thus, the Gaussian sequence is admissible.
Remark 3. While the family of Gaussian weights has the distinct advantage of being universally admissible independently of the interval size a, it can seem conceptually awkward to use the exponential function to prove the existence of the exponential function. Although this approach is formally and logically sound, from a didactic perspective it is highly desirable to construct admissible sequences that do not rely on the exponential function at all. This ensures that our framework can serve as a rigorous, self-contained primary definition of the exponential function in real and complex analysis. A simple and elegant alternative is the family of polynomial power weights. Note, however, that these weights are locally tailored, meaning the choice of the sequence now explicitly depends on a. Proposition 2 (Polynomial Power Weights). For any fixed a > 0, choose a constant c such that 0 < c < 8/a2 . Then the sequence of polynomial power weights ωλ (x) = (1 + cx2 )λ ,
λ ∈ N,
(24)
is an admissible weight sequence on the interval Ia = [−a, a]. Proof. Let Ia = [−a, a] be an arbitrary compact interval. We verify the three conditions of Definition 1. The lower bound (15) is trivial: since c > 0, we have ωλ (x) ⩾ 1 for all x ∈ Ia . For the linear suppression, the ratio is R1 (x, λ) = (C(x2 ))λ , where for every t ⩾ 0 the function C is given by: C(t) =
1 + ct/4 3c =1− t. 1 + ct 4(1 + ct)
Clearly C(0) = 1. For t ∈ [0, a2 ], we can bound the denominator from above by 1 + ca2 . Therefore: C(t) ⩽ 1 − αt,
where
3c > 0. 4(1 + ca2 )
α :=
Using the standard inequality 1 − s ⩽ e−s for s ⩾ 0, we obtain C(t) ⩽ e−αt , and thus R1 (x, λ) ⩽ 2 e−αλx ⩽ 1. Applying the exponential bound (22) with the constant α, the linear suppression holds: 2
sup [|x|R1 (x, λ)] ⩽ sup [|x|e−αλx ] ⩽ √ x∈Ia
x∈Ia
1 λ→∞ −−−→ 0. 2αeλ
For the quadratic suppression (19), the ratio simplifies to R2 (x, λ) = (B(x2 ))λ , with the base: B(t) =
ct(8 − ct) (1 + ct/4)2 =1− . 1 + ct 16(1 + ct)
By our strict choice of c < 8/a2 , the term (8 − ct) is strictly bounded from below by (8 − ca2 ) > 0 for all t ∈ [0, a2 ]. Bounding the denominator from above by 16(1 + ca2 ), we obtain: B(t) ⩽ 1 − βt,
where
β :=
c(8 − ca2 ) > 0. 16(1 + ca2 ) 2
Using 1 − s ⩽ e−s again, B(t) ⩽ e−βt , which gives R2 (x, λ) ⩽ e−βλx . Applying (22) with the constant β yields: 2
sup [x2 R2 (x, λ)] ⩽ sup [x2 e−βλx ] ⩽ x∈Ia
x∈Ia
Thus, the polynomial power sequence is locally admissible.
1 λ→∞ −−−→ 0. βeλ □
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
4.2. Proof of Lemma 1. We define the operator T acting on functions q : R → C by x γ2 1 x x2 x q (T q)(x) = + +γ + q2 , 4 2 4 2 16 2
9
(25)
so that the proof reduces to the existence of q satisfying (9) to the fixed point problem T q = q. To prove uniqueness and existence of the fixed point, we verify that for every a > 0 there exists a sufficiently large λ > 0 such that the hypotheses of theBanach–Caccioppoli contraction principle (Theorem 1) apply to the Banach space C 0 (Ia , C), ∥ · ∥λ , yielding the existence and uniqueness of a fixed point of T in Bλ,1 . We proceed in three steps. Step 1: Self-Mapping (T (Bλ,1 ) ⊆ Bλ,1 ). We verify that for any admissible weight sequence, choosing λ sufficiently large ensures T maps the unit ball into itself. Let q ∈ Bλ,1 , which implies |q(x)| ⩽ ωλ (x). Noting that |γ|2 = 1 because γ ∈ {1, i}, we apply the triangle inequality to (25) to estimate: 1 1 |x| x x 2 x 2 |(T q)(x)| ⩽ + ωλ + |γ| + ωλ . 4 2 4 2 16 2 Dividing by ωλ (x) to switch to the weighted norm, we obtain: 1 1 |x| x2 |(T q)(x)| ⩽ + + |γ| R1 (x, λ) + R2 (x, λ) ωλ (x) 4ωλ (x) 2 4 16 with R1 , R2 defined by (16), (18). Taking the supremum over Ia and applying the admissibility properties (15) and (17) (i.e., ωλ ⩾ 1 and R1 ⩽ 1): ∥T q∥λ ⩽
1 1 |γ| 1 + + sup [|x|R1 (x, λ)] + sup [x2 R2 (x, λ)] . 4 2 4 x∈Ia 16 x∈Ia
As λ → ∞, the last two terms vanish by the properties of the admissible weights, bounding the norm by 3/4 + o(1) < 1 for λ sufficiently large. Thus, T (Bλ,1 ) ⊆ Bλ,1 for such values of λ. Step 2: Strict Contraction on Compact Intervals. Let q1 , q2 ∈ Bλ,1 . We evaluate the difference: x h x x i 1 x x2 x (T q1 )(x) − (T q2 )(x) = q1 + q2 q1 − q2 . +γ + 2 4 16 2 2 2 2 Since both functions belong to Bλ,1 , we bound the sum by |q1 (x/2) + q2 (x/2)| ⩽ 2ωλ (x/2). Extracting the norm of the difference, we have |q1 (x/2) − q2 (x/2)| ⩽ ωλ (x/2)∥q1 − q2 ∥λ . Taking the absolute value and dividing the entire expression by ωλ (x) yields: |x| |(T q1 )(x) − (T q2 )(x)| 1 x2 ⩽ + |γ| R1 (x, λ) + R2 (x, λ) ∥q1 − q2 ∥λ ωλ (x) 2 4 8 Taking the supremum over Ia , we obtain the Lipschitz constant κ(λ): 1 |γ| x2 κ(λ) = sup R1 (x, λ) + |x|R1 (x, λ) + R2 (x, λ) . 4 8 x∈Ia 2 Because R1 (x, λ) ⩽ 1 (see (17)), the first term is strictly bounded by 1/2. By the admissibility properties of the weight, the second and third terms vanish uniformly on Ia as λ → ∞. Consequently: lim sup κ(λ) ⩽ λ→∞
1 < 1. 2
Therefore, for sufficiently large λ, T is a strict contraction on Bλ,1 . By the Banach-Caccioppoli theorem, there exists a unique fixed point qa ∈ Bλ,1 .
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
10
Step 3: Global Consistency and Uniqueness. We have proven that for any a > 0, there exists a unique continuous solution qa on Ia . If 0 < a1 < a2 , the restriction of qa2 to Ia1 must be a solution to the fixed point problem on Ia1 . By uniqueness, this restriction must coincide with qa1 . This compatibility allows us to define a unique global function q∗ : R → C by setting q∗ (x) = qa (x) for any x ∈ Ia . Furthermore, the same contraction argument holds if we replace C 0 (Ia , C) with the Banach space of bounded functions B(Ia , C). Thus, q∗ is the unique solution even among bounded functions. □ 4.3. Proof of Lemma 2. In this subsection we prove Lemma 2 whose utility is in showing that there is a one to one correspondence between the set of quadratic residuals and the corresponding “exponential” functions e satisfying axioms (A1 , A2 ). Note that, at this stage, we still don’t know that a function e satisfying these axioms exists, nor that it is unique. Proof. For x ̸= 0, the function q is uniquely determined by inverting (10): q=
e(x) − (1 + γx) . x2
(26)
Since e is continuous, q is continuous on R \ {0}. To ensure continuity on the whole of R, we must verify that the limit exists as x → 0 and equals q(0). We proceed in three steps. Step 1: Functional Equation. We first verify that q satisfies (9) for x ̸= 0. Substituting the expression (10) into the duplication identity e(x) = e2 (x/2), we obtain: 2 x x2 x 2 . (27) 1 + γx + x q(x) = 1 + γ + q 2 4 2 Expanding the square on the right-hand side, cancelling the linear term 1 + γx, and dividing by x2 yields exactly the functional equation (9). We define q(0) = γ 2 /2 so that the equation holds at x = 0 as well. However, we must rule out the possibility that q is unbounded near the origin. Step 2: Boundedness near the origin. We show that any residue q : R → C satisfying the ansatz must be bounded around the origin. From the differentiability of e at 0 (Axiom A2 ), we deduce that xq(x) → 0 as x → 0. Consequently, there exists a constant c > 0 and an interval Iδ := (−δ, δ) with δ > 0, such that |xq(x)| ⩽ c for every x ∈ Iδ . Using the functional equation (9), we can bound |q| as follows 1 x x x 2 1 x 1 1 |q(x)| ⩽ q + |γ| + q ⩽ q + (|γ| + c)2 . 2 2 4 2 2 2 2 4 We set κ := 41 (|γ| + c)2 and iterate the previous inequality j times to yield, at least for every x ̸= 0, |q(x)| ⩽
j−1 j−1 X X 1 x 1 1 x x 1 q = q . + κ + κ 2j 2j 2i |x| 2j 2j 2i i=0
i=0
Taking the limit as j → ∞, the first term vanishes (because aj q(aj ) → 0 for any sequence (aj )j∈N that converges to zero), and the geometric series converges to 2κ. Thus |q(x)| ⩽ 2κ for every x ∈ Iδ , proving that q is uniformly bounded near 0. Step 3: Uniqueness of the residue. So far, we have shown that q is a global solution of the functional equation (9) that is continuous on R \ {0} and bounded on Iδ . By Lemma 1, there exists a unique bounded solution q∗ of the functional equation(9), which is known to be continuous everywhere. Therefore q = q∗ on Iδ . This implies that q is continuous at the origin as well. □ 5. Fixed point characterization of the Natural Logarithm In this section, we prove that the natural logarithm is the unique fixed point of a suitably defined linear operator acting on a carefully constructed complete metric space.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
11
5.1. Definition of the Operator and the Functional Setup. Let I = [1/A, A] be a closed interval with A > 1. We define the linear functional operator T as: √ T [f ](x) = 2f ( x) . (28) √ √ √ √ The transformation x 7→ x maps the interval I strictly into itself, since x ∈ [1/ A, A] ⊂ I for all x ∈ I. As a consequence, all iterates of T remain within the domain I, and the operator is well-defined on any subspace of Cb (I, R), the Banach space of bounded continuous functions on I. Also, the functional equation T [f ] = f naturally admits a one-parameter family of solutions, f (x) = c log x for c ∈ R. To isolate the natural logarithm, corresponding to the normalization c = 1, we impose the initial conditions f (1) = 0 and f ′ (1) = 1 (cf. axioms (B1 , B2 ) in (3) and (4)). Within our topological framework, these conditions are incorporated by introducing the affine reference function f0 (x) = x−1 and by working in a suitably weighted function space. We now introduce the functional framework in which the fixed point analysis will be carried out. Let B (I, R) denote the subspace of Cb (I, R) consisting of all continuous functions h : I → R such that h(1) = 0 and |h(x)| · x < +∞ . (29) ∥h∥B := sup 2 x∈I\{1} (x − 1) This weighted norm captures the quadratic vanishing of functions at x = 1 and will play a crucial role in controlling the action of the operator T . We also need the space M (I, R) as the affine translation of B (I, R) by the reference function f0 (x) = x − 1: M (I, R) := {f : I → R : f − f0 ∈ B} . (30) Lemma 3. The space (B (I, R) , ∥ · ∥B ) is a Banach space. Therefore its affine translation M (I, R) = f0 + B (I, R) is a closed subset of the ambient functional space and is therefore a complete metric space under the induced metric d(f1 , f2 ) = ∥f1 − f2 ∥B . Proof. It is clear that B (I, R) is a normed space. It remains to prove completeness. Let (hn )n∈N be a Cauchy sequence in B (I, R). We need to show the existence of an element h ∈ B (I, R) such that hn → h in B (I, R). x For each n ∈ N, define the auxiliary function Hn (x) := hn (x) (x−1) 2 on the punctured compact interval I \ {1}. By the definition of the norm on B (I, R), we have for any n, m ∈ N: ∥hn − hm ∥B = sup |Hn (x) − Hm (x)| = ∥Hn − Hm ∥∞ .
(31)
x∈I\{1}
This isometry implies that (Hn )n∈N is a Cauchy sequence in the Banach space Cb (I\{1}, R) of bounded, continuous functions on I \ {1}, endowed with the supremum norm. Consequently, Hn converges uniformly to a bounded and continuous function H (in particular, its supremum norm is finite: ∥H∥∞ < ∞). We now define a candidate limit function h : I → R by (x − 1)2 H(x) , x ̸= 1, (32) h(x) := x 0, x = 1. Since x ∈ I = [1/A, A] implies x > 0, this definition is meaningful. Moreover, h is continuous on I \ {1} because H is continuous there. To check continuity at x = 1, note that (x − 1)2 (x − 1)2 ⩽ ∥H∥∞ −→ 0 as x → 1. (33) x x Thus, h ∈ Cb (I, R). Next, we verify that h ∈ B (I, R). We have already established h(1) = 0. Furthermore, evaluating its norm yields: |h(x) − h(1)| = |h(x)| = |H(x)|
|h(x)| · x = sup |H(x)| = ∥H∥∞ < ∞ . 2 x∈I\{1} x∈I\{1} (x − 1)
∥h∥B = sup
(34)
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
12
So h ∈ B (I, R). Moreover, ∥hn − h∥B = ∥Hn − H∥∞ −→ 0 because Hn → H uniformly. This proves that every Cauchy sequence in B (I, R) converges, and thus (B (I, R) , ∥ · ∥B ) is a Banach space. □ Theorem 3 (The natural logarithm). Let I = [1/A,√A] be a closed interval with an arbitrary constant A > 1. Then, the linear operator T [f ](x) = 2f ( x) is a strict global contraction on the complete metric space M (I, R) with contraction constant κ = 1/2. Consequently, T admits a unique fixed point in M (I, R), and this fixed point is the natural logarithm. Remark √ 4. The fact that the natural logarithm log x is a fixed point is immediate: T [log](x) = 2 log( x) = log x. It remains only to check that log ∈ M (I, R). This amounts showing that h∗ (x) := log x − (x − 1) belongs to B (I, R). Since h∗ (1) = 0, it remains to show that |h∗ (x)| · x < +∞. 2 x∈I\{1} (x − 1)
∥h∗ ∥B = sup
(35)
Using the Taylor expansion of log x at x = 1, we get h∗ (x) = − 21 (x − 1)2 + o((x − 1)2 ). Therefore |h∗ (x)| · x (x − 1)2
(36)
extends continuously to x = 1 with finite limit 1/2. Being continuous on the compact interval I, the supremum is finite, confirming h∗ ∈ B (I, R) and log ∈ M (I, R). Because the Banach fixed-point theorem yields a unique fixed point of T in M (I, R), and log is one such fixed point, it follows that log x is the unique fixed point of T in M (I, R). Proof. We divide the proof in two steps Step 1: The operator T maps M into itself (T (M) ⊂ M). Let f ∈ M (I, R). Then f = f0 + h for some h ∈ B (I, R). We compute √ √ √ √ T [f ](x) − f0 (x) = 2( x − 1 + h( x)) − (x − 1) = −( x − 1)2 + 2h( x) (37) Thus it is enough to √ show that both terms on the right-hand side belong to B (I, R). For the first term, define g(x) = −( x − 1)2 . Then g(1) = 0, and for x ̸= 1, x ∥g∥B = sup √ . 2 x̸=1 ( x + 1) This function is continuous and clearly bounded√on I (approaching 1/4 as x → 1), hence g ∈ B (I,√R). For the second term, let h ∈ B (I, R). Since x ∈ I whenever x ∈ I, √ the composite x 7→ h( x) is well defined on I. In Step 2 below we show that the map h → 7 2h( ·) sends B (I, R) into itself. √ Hence 2h( x) ∈ B (I, R). Therefore, T [f ] − f0 ∈ B and T [f ] ∈ M. Step 2: Proof of Strict Contraction. Let f1 , f2 ∈ M (I, R) and set h = f1 − f2 ∈ B (I, R). By linearity, d(T [f1 ], T [f2 ]) = ∥T [h]∥B . We explicitly compute this norm: √ 2|h( x)| · x ∥T [h]∥B = sup . 2 x∈I\{1} (x − 1) √ √ √ Applying the substitution y = x and setting J := [1/ A, A] ⊂ I, by the elementary inequality (y + 1)2 ⩾ 4y, we obtain 2y 1 |h(y)| · y 1 |h(y)| · y ∥T [h]∥B = sup · ⩽ sup = ∥h∥B . 2 2 2 (y + 1) 2 y∈I\{1} (y − 1) 2 y∈J\{1} (y − 1) Therefore, T is a strict global contraction on B (I, R) with contraction constant κ = 1/2. Finally, since B (I, R) is Banach, the affine space M (I, R) = f0 + B (I, R) is complete with respect to the metric d(f, g) = ∥f − g∥B . The Banach fixed-point theorem therefore applies and gives a unique fixed point of T in M.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
13
As observed in the remark, log x ∈ M and satisfies T [log] = log. Hence this unique fixed point is precisely the natural logarithm. □ 6. Algorithms and Numerical Results This section evaluates the numerical implementation of the methods developed in the preceding sections. We refer to the computational routines implementing the real exponential, complex exponential, and natural logarithm as kernels. 6.1. Overview. All three kernels follow a fixed-point construction pattern comprising three stages: argument reduction to a small interval, polynomial approximation (seed) of the residual function on that interval, and reconstruction of the final value. The design space is governed by two parameters: the degree K of the polynomial seed and the number N of fixed-point iterations. Because each iteration halves the interval, a higher degree can be traded for more iterations at roughly equal accuracy. On general-purpose CPUs, high-degree/low-iteration configurations generally win, whereas architectures with limited multiplier width or storage favor short-seed/iterated configurations. This fixed point approach replaces traditional table based argument reduction with purely arithmetic transformations, using repeated halvings for the exponential and repeated square roots for the logarithm. It relies exclusively on IEEE 754 arithmetic primitives, including fused multiply add (FMA) and direct exponent manipulation, without memory lookups, making it particularly well suited to vectorization. The choice of seed differs structurally between functions. The residuals of the exponential and sine–cosine functions are entire; their Taylor series converge quickly, so simple truncated series are used. Conversely, the logarithm’s residual converges slowly, requiring a fitted minimax polynomial to save degrees. Accuracy is reported in ulp (units in the last place) against an 80-bit extended-precision reference. Speed is reported as a ratio to the SLEEF vector library (grades u35 and u10). All kernels are AVX2, table-free, and run in plain double precision unless marked double-double (dd). Note that these are research prototypes assuming finite, normal, in-range inputs, whereas SLEEF handles all IEEE special values; the comparison is therefore generous to our prototypes by construction, and its outcomes are read with this boundary in mind. 6.2. Stability and Error Propagation. The numerical stability of a fixed-point kernel depends heavily on the coordinate system used to represent its discrete orbit. While the underlying functional operators are contractive analytically, finite-precision arithmetic injects rounding errors that can either be amplified or attenuated during the reconstruction stage. 6.2.1. The Exponential: Contractive Reconstruction. The exponential kernel leverages the homomorphism ea+b = ea eb to reduce the argument to a small value t0 , evaluates a polynomial seed, and then ascends back to the original argument using N successive doublings (t 7→ 2t). Crucially, the kernel computes this entirely within the residual coordinate q(t) = (et − 1 − t)/t2 which removes the constant and linear contributions and isolates the nonlinear remainder. In this coordinate, the doubling map is represented by the rational recurrence 1 1 t t2 Dq (q) = + + q + q2. (38) 4 2 2 4 The reconstruction process is computationally contractive. Indeed, differentiating the recurrence with respect to q and evaluating at the exact residual coordinate q(t) gives ∂q Dq |q=q(t) = et /2. Since the argument reduction guarantees that the reduced argument remains √ in the interval |t| ⩽ log 2/2 < 1/2, the contraction factor is uniformly bounded by elog 2/2 /2 = 2/2 < 1. Therefore, perturbations introduced by floating point rounding during an ascent step are geometrically attenuated by subsequent doubling operations. The resulting error accumulation remains bounded, and the attainable accuracy is essentially independent of the number of reconstruction steps N . In contrast, the naive approach of storing y = et and repeatedly squaring it (y 7→ y 2 ) has a fundamental drawback. Each squaring step doubles the relative error, causing perturbations to grow
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
N
√ logarithm, T [f ] = 2f ( ·)
0 1 2 3 4 5 6
1× 2× 4× 8× 16× 32× 64×
14
exponential, naive y 7→ y 2 exponential, residual 1× 2× 4× 8× 16× 32× 64×
1.00× 0.59× 0.32× 0.17× 0.09× 0.04× 0.02×
Table 1. Error propagation through the reconstruction. A seed perturbation ε0 = 10−12 propagated through N reconstruction steps. The logarithm’s ascent multiplies it by 2N ; the naive exponential, which stores the value y = et and ascends by squaring, amplifies the relative error by exactly the same factor 2N ; the contractive residual coordinate attenuates it instead. The middle column shows that the amplification is a property of the coordinate, not of the function: the same exponential, in the value coordinate, is precisely as unstable as the logarithm.
by a factor of 2N after N iterations. The residual coordinate replaces this simple multiplication with a slightly more expensive rational recurrence, but in return it ensures that the reconstruction phase attenuates rather than amplifies the errors introduced at each step (see Table 1). Remark 5 (Relation to scaling and squaring). The reconstruction phase is an instance of the classical scaling and squaring paradigm underlying modern algorithms for the matrix exponential [Hig05, AMH09]. The argument is first scaled by a power of two, an approximation is constructed on the reduced interval, and the original value is recovered by repeated squaring. Standard implementations of this paradigm, such as MATLAB’s expm, execute this upward phase directly on the macroscopic matrix value (X 7→ X 2 ). As for the scalar case (see Table 1), this macroscopic reconstruction amplifies injected rounding errors by a factor of two at each step. Consequently, current matrix algorithms must restrict the scaling depth and rely on high-degree Padé approximants to control the backward error. Formulating the matrix ascent in a contractive residual coordinate (analogous to the scalar recurrence) could suppress this error amplification; we leave this as an open direction. 6.2.2. The Logarithm: Stabilization via the Product Form. The logarithm is evaluated by first partitioning the input using its IEEE 754 representation. Any x > 0 is √ expressed as x = 2E m with m ∈ [1, 2). A conditional scaling, m 7→ m/2 and E 7→ E + 1 for m ⩾ 2, recenters the mantissa to the symmetric interval m ∈ [2−1/2 , 21/2 ). Folding one bit of the exponent into a half-integer count eh further restricts the domain to m ∈ [2−1/4 , 21/4 ). These operations modify only the exponent field and therefore incur no rounding error. The logarithm then decomposes additively as log x = eh log2 2 + log m,
(39)
reducing the problem to evaluating log m on a bounded interval near 1, followed by the addition of the exponent term eh log2 2 . Unlike the exponential kernel, where the recurrence attenuates rounding N errors at each step, the logarithm iteration relies on the identity log m = 2N log(m1/2 ). Because the reconstruction of log m involves multiplying the polynomial approximation by the scalar 2N , numerical errors introduced during the approximation phase are amplified by this factor. To mitigate this error growth, the iterative sequence is reformulated. The state variable is defined as the deviation from unity, u0 = m − 1, rather than evaluating the argument m directly. Because the recentered mantissa m lies within [0.5, 2.0], Sterbenz’s lemma guarantees that this initial floating-point subtraction is computed exactly.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
N
naive
product
2 3 4 5 6
10.6 18.4 31.3 58.2 115.3
5.9 5.9 5.9 5.9 5.9
15
Table 2. Naive versus product logarithmic descent. Maximum absolute error (×10−16 ) of the scalar logarithm over [10−3 , 103 ] with a degree-10 seed. The naive form grows exponentially, whereas the product form maintains a flat error floor.
k
The fixed-point iteration computes a sequence of nested square roots, sk = m1/2 , which converges to 1. Tracking sk directly leads to numerical instability, as the final logarithm polynomial requires evaluating the difference sk − 1, which is subject to subtractive cancellation. To prevent this, the coordinate uk := sk − 1 is tracked at each step. Computing uN = sN − 1 directly at the end of the iteration would still incur subtractive cancellation, and the subsequent reconstruction would amplify the resulting error by 2N . This is resolved by applying the difference of squares, s2k − 1 = (sk − 1)(sk + 1), yielding the recurrence uk = uk−1 /(sk + 1). This leads to: u0 uN = QN (40) k=1 (sk + 1) As sk → 1, each denominator factor approaches 2. This sequence of divisions reduces the exact initial deviation u0 by a factor of approximately 2N . The product formulation replaces subtractions with additions and divisions. Consequently, when the final reconstruction step scales this attenuated value by 2N , it offsets the prior division, preventing the error amplification that characterizes the direct evaluation (see Table 2). 6.2.3. Functional Formulations and Pseudocode. Algorithms 1, 2, and 3 detail the pointwise execution paths. All three share the fixed-point architecture outlined in the overview (§6.1): argument reduction (Cody–Waite for the exponential family, quarter-octave for the logarithm), evaluation of a degree-K polynomial seed, and N reconstruction iterations, carried out in the residual coordinate for the exponential family and in the stabilized product form for the logarithm. The complex exponential additionally applies a final quadrant selection, mapping the result on the reduced argument r ∈ [−π/4, π/4] back to the full circle. The implementations rely exclusively on two standard, table-free arithmetic primitives: • Cody–Waite Reduction: Computes the reduced argument r = x − kC (where C is an irrational period such as log 2 or π/2) without catastrophic cancellation. To avoid floating-point precision loss, we use standard techniques consisting of representing the constant in simulated extended precision as an unevaluated sum, C ≈ Chigh + Clow . This allows the dominant subtraction x − k · Chigh to be computed exactly, preserving the significant digits of the mathematical remainder before safely applying the small correction k · Clow . • Exponent Scaling (scalbn): Restores the 2k factor isolated during the reduction step. Instead of explicitly computing 2k and performing a floating-point multiplication, which introduces additional rounding errors and may cause intermediate overflow, the integer k is directly added to the IEEE 754 binary exponent field of the operand. This results in an exact scaling operation with O(1) computational cost. By utilizing these exact operations, the reduction and scaling phases avoid introducing additional rounding error, preserving the numerical accuracy of the core approximations detailed in the algorithms.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
16
Algorithm 1 Fixed-point real exponential (residual form) Require: x ∈ R; parameters N, K ∈ N Ensure: y ≈ ex 1: k ← round(x/ log 2) 2: r ← fma(k, − log 2high , x) − k log 2low ▷ Cody–Waite two-limb reduction −N 3: t ← r · 2 ▷ Exact bit-shift decrement of the floating-point exponent PK j a t ▷ Truncated Taylor series of q(t) = (et − 1 − t)/t2 , aj = 1/(j + 2)! 4: q ← j=0 j 5: loop (N times) 2 6: a ← 21 + 2t , b ← t4 7: q ← 14 + a q + b q 2 ▷ Contractive residual recurrence Dq , Eq. (38) 8: t ← 2t ▷ Ascend one level toward r 9: end loop 10: return scalbn(fma(r 2 , q, 1 + r), k) ▷ Reconstruction 2k (1 + r + r2 q) Algorithm 2 Fixed-point complex exponential (sine and cosine together) Require: x ∈ R; parameters N, K ∈ N Ensure: (cos x, sin x) 1: κ ← round(x · 2/π) mod 4 ▷ Quadrant index 2: r ← CodyWaiteReduce π/2 (x) ▷ Three-limb reduction yielding r ∈ [−π/4, π/4] 3: t ← r · 2−N ▷ Exact scaling of the reduced argument P⌊K/2⌋ 2j 4: qr ← ▷ Truncated Taylor seed: even part qr = (cos t − 1)/t2 j=0 b2j t P⌊(K−1)/2⌋ 5: qi ← b2j+1 t2j+1 ▷ odd part qi = (sin t − t)/t2 j=0 6: loop (N times) ▷ Pointwise evaluation of the complex operator equation 7: hr ← −0.25 + 0.5qr − 0.5tqi + 0.25t2 (qr2 − qi2 ) 8: hi ← 0.5qi + 0.5tqr + 0.5t2 qr qi 9: qr ← hr , qi ← hi , t ← 2t 10: end loop 11: c ← 1 + r 2 qr , s ← r + r 2 qi ▷ (cos r, sin r) on the reduced argument 12: return QuadrantSelect(κ, c, s) ▷ Recover (cos x, sin x) from (cos r, sin r) and the quadrant κ Algorithm 3 Fixed-point logarithm (product form) Require: x ∈ R+ ; parameters N, K ∈ N Ensure: y ≈ log x 1: (m, eh ) ← QuarterOctaveReduce(x) ▷ Decomposition mapping m ∈ [2−1/4 , 21/4 ) 2: u0 ← m − 1.0 ▷ Exact evaluation via Sterbenz’s lemma 3: s ← m, D ← 1.0 4: loop (N times) Q √ 5: s ← s, D ← D · (s + 1.0) ▷ Accumulate the telescoping denominator D = (sk + 1) 6: end loop 7: u ← u0 /D ▷ A single division executing range-reduction to the origin PK j 8: P ← c u ▷ Minimax polynomial seed of log(1 + u)/u j=0 j log 2 N 9: return fma(eh , 2 , 2 uP ) ▷ Reconstruction; log2 2 scales eh , carried in half-integer units 6.3. The Vector Regime: Throughput Characteristics. Under SIMD vectorization, where a single instruction operates concurrently across multiple register lanes (e.g., four double-precision operands), eliminating table lookups avoids irregular memory traffic. Rather than paying the latency of gather operations at data-dependent addresses, table-free kernels run entirely within registers.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
17
function
prec.
N
K
max ulp
vs sleef
exp exp exp exp
plain plain plain plain
0 1 2 3
12 10 8 7
0.67 (u10) 0.65 (u10) 0.66 (u10) 0.66 (u10)
1.08× u10† 0.91× u10† 0.81× u10† 0.77× u10†
sin, cos sin, cos sin, cos sin, cos
plain plain plain plain
0 1 2 3
8 6 5 4
1.42 (u10∗ ) 1.42 (u10∗ ) 1.44 (u10∗ ) 2.21 (u35)
2.38× u35, 3.85×u10 1.88× u35, 3.05×u10 1.49× u35, 2.41×u10 1.24× u35, 2.00×u10
log log log log log
plain plain plain plain dd
0 1 2 3 0
15 12 10 10 14
1.99 (u35) 3.21 (u35) 3.54 (u35) 3.32 (u35) 0.59 (u10)
1.64× u35, 2.21×u10 1.27× u35, 1.72×u10 1.16× u35, 1.57×u10 0.85× u35, 1.14×u10 1.16× u35, 1.57×u10
Table 3. Fixed-point kernels versus SLEEF (AVX2). Accuracy is reported as maximum ulp error, categorized by grade: u10(⩽ 1), u35 (⩽ 3.5), and u10∗ (just above u10). Speed ratios are measured against the respective SLEEF grade († the exp kernel compares to u10, as SLEEF lacks a u35 routine). Full bordeaux indicates configurations that are faster and meet the target grade; light bordeaux indicates faster configurations that exceed the error threshold; plain black denotes slower performance. Seeds are Taylor polynomials for exp and sin, cos, and minimax for log. As established in §6.2, the rising error floor of the log kernel reflects its amplifying macroscopic reconstruction.
To evaluate this approach, the AVX2-vectorized prototypes were benchmarked against the SLEEF library. As detailed in §6.1, accuracy is reported in units in the last place (ulp) against an 80-bit reference, and speed is measured as a relative throughput ratio against SLEEF’s u10 (⩽ 1 ulp) and u35 (⩽ 3.5 ulp) baselines. The polynomial seeds approximate the analytic residuals of each coordinate system: q(t) = t (e − 1 − t)/t2 for the real exponential, qr (y) = (cos y − 1)/y 2 and qi (y) = (sin y − y)/y 2 for the complex exponential, and a minimax approximation of log(1 + u)/u for the logarithm. The logarithmic residual is evaluated on the reduced mantissa m = 1 + u, where the deviation u = m − 1 is computed exactly. To preserve relative accuracy near m = 1 (where u → 0), we employ the standard technique of factoring out the root. Instead of approximating the logarithm directly, we approximate the well-conditioned quotient log(1 + u)/u with a polynomial. Multiplying this polynomial by u recovers the logarithm and guarantees bounded relative error, whereas a direct polynomial approximation of log(1 + u) would leave a tiny non-zero absolute error at u = 0, causing the relative error to diverge. Table 3 indicates that the performance and accuracy of these portable, table-free prototypes are comparable to those of a standard production library. Their performance and accuracy profiles depend on the algebraic structure of each function’s reconstruction stage: • Real Exponential (No Cancellation): Reconstructed via er = 1+r+r2 q, the calculation avoids subtractive cancellation. This preserves relative precision, allowing standard double precision to achieve u10 accuracy (≈ 0.65 ulp). Because the reconstruction is contractive, the error floor remains constant across iteration depths N ∈ {0, 1, 2, 3}. This decoupling permits the selection of a design point according to hardware constraints: N = 0 maximizes throughput on general-purpose CPUs (1.08× SLEEF), while higher N values reduce coefficient storage requirements on constrained architectures. • Sine and Cosine (Local Cancellation): Reconstructed via cos x = 1+r2 qr and sin x = r+r2 qi , local cancellation near the zeros limits the global precision to approximately 1.42 ulp (denoted as u10∗ ). Away from these zeros, the same seed yields an error of 0.68 ulp, indicating
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
18
that the error bound is determined by the algebraic structure rather than the polynomial approximation. This kernel yields higher throughput than SLEEF across all evaluated iteration depths, showing a 2× to 4× increase compared to the strict u10 baseline, with an accuracy reduction of approximately 0.5 ulp. • Logarithm (Macroscopic Shift and Amplification): Reconstructed via log x = eh log2 2 + log m, the overall precision is limited by the addition of the exponent term eh log2 2 . In standard double precision, this operation limits accuracy to approximately 2 ulp at N = 0. As the iteration depth N increases, this reconstruction structure amplifies intermediate rounding errors by a factor of 2N . Consequently, the error increases (1.99 → 3.54 ulp) and exceeds the u35 threshold at N = 2. Achieving u10 accuracy (0.59 ulp) requires evaluating the addition in double-double (dd) precision at N = 0; this configuration maintains a higher throughput than the corresponding SLEEF baselines. • Logarithm (Macroscopic Shift and Amplification): Reconstructed via log x = E log 2 + log m, the overall precision is limited by the addition of the exponent term E log 2. In standard double precision, this operation limits accuracy to approximately 2 ulp at N = 0. As the iteration depth N increases, this reconstruction structure amplifies intermediate rounding errors by a factor of 2N . Consequently, the error increases (1.99 → 3.54 ulp) and exceeds the u35 threshold at N = 2. Achieving u10 accuracy (0.59 ulp) requires evaluating the addition in double-double (dd) precision at N = 0; this configuration maintains a higher throughput than the corresponding SLEEF baselines. In summary, the vector implementation characteristics follow a consistent pattern within the (K, N ) design space. The iteration depth N and the polynomial degree K can be adjusted inversely while maintaining approximately constant accuracy. This allows the kernels to be adapted to specific hardware constraints, such as balancing multiplier utilization against register usage. However, regardless of the chosen (K, N ) configuration, the algebraic structure of the final reconstruction formula determines the minimum achievable error in finite precision. 6.3.1. A Reduced-Accuracy, High-Throughput Regime. While standard production libraries define accuracy targets for general numerical computing, specific applications such as computer graphics, signal processing, and machine learning can operate with error tolerances near 10−6 (comparable to single precision). Standard double-precision libraries typically do not provide performance optimizations for this reduced-accuracy regime. The fixed-point methods can target this regime by reducing the degree of the polynomial seed, K. As shown in Table 4, increasing the iteration depth N with these lower-degree seeds reduces the approximation error by a factor of approximately 2−(K+2) per step for the contractive kernels (exponential and sine/cosine). For the logarithm, the error bounds near 10−9 when the polynomial approximation error intersects the error amplification inherent to its reconstruction phase. Adjusting the (K, N ) parameters allows the kernels to traverse the accuracy-throughput design space. 6.3.2. The (K, N ) Trade-Off as a Hardware Portability Axis. Inspection of Table 3 reveals an accuracy preserving frontier in the (K, N ) parameter space. For the real exponential kernel, configurations ranging from (N =0, K=12) to (N =3, K=7) achieve approximately 0.65 ulp accuracy. This suggests that the degree of the polynomial seed approximation can be reduced by increasing the number of reconstruction iterations, without affecting the final accuracy. This observation can also be justified by the formal error estimates developed above, although we do not pursue this analysis further here. The existence of this frontier is a direct consequence of the contractive residual coordinate. Since the reconstruction phase attenuates rather than amplifies rounding errors, increasing the number of iterations affects the computational cost but does not significantly alter the attainable accuracy. Consequently, configurations along the frontier are numerically equivalent but have different computational characteristics, allowing the implementation to select the most suitable configuration for the target architecture.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
function
prec.
N
K
max error
vs sleef
3.0 × 10−6 (rel.)
exp exp exp exp exp exp exp exp
plain plain plain plain plain plain plain plain
1 2 3 4 1 2 3 4
2 2 2 2 4 4 4 4
1.7 × 10−7 (rel.) 1.1 × 10−8 (rel.) 6.5 × 10−10 (rel.) 2.2 × 10−9 (rel.) 3.1 × 10−11 (rel.) 4.7 × 10−13 (rel.) 7.3 × 10−15 (rel.)
2.01× u10† 1.49× u10† 1.18× u10† 0.98× u10† 1.80× u10† 1.39× u10† 1.10× u10† 0.90× u10†
sin, cos sin, cos sin, cos sin, cos sin, cos sin, cos sin, cos sin, cos
plain plain plain plain plain plain plain plain
1 2 3 4 1 2 3 4
2 2 2 2 4 4 4 4
9.2 × 10−6 (abs.) 2.6 × 10−7 (abs.) 7.6 × 10−9 (abs.) 2.3 × 10−10 (abs.) 4.4 × 10−11 (abs.) 7.7 × 10−14 (abs.) 2.3 × 10−16 (abs.) 1.3 × 10−16 (abs.)
2.42× u35 1.83× u35 1.39× u35 1.11× u35 2.15× u35 1.61× u35 1.24× u35 1.01× u35
log log log log log log log log
plain plain plain plain plain plain plain plain
1 2 3 4 1 2 3 4
4 4 4 4 5 5 5 5
2.5 × 10−7 (abs.) 1.6 × 10−8 (abs.) 1.9 × 10−9 (abs.) 1.0 × 10−9 (abs.) 9.6 × 10−9 (abs.) 1.2 × 10−9 (abs.) 9.6 × 10−10 (abs.) 9.5 × 10−10 (abs.)
1.96× u35 1.25× u35 0.92× u35 0.72× u35 1.86× u35 1.25× u35 0.91× u35 0.72× u35
19
Table 4. Short-seed configurations for the 10−6 -accuracy regime. Metrics denote relative error for exp (whose output spans many binades) and absolute error for sin, cos, and log (whose outputs are of moderate size and cross zero). Two seed degrees (K) are evaluated across iteration depths N = 1, . . . , 4. For the contractive kernels, the error falls geometrically by ≈ 2−(K+2) per step, whereas the log kernel hits an amplification floor near 10−9 (§6.2). Ratios compare against SLEEF’s lowest-accuracy tier († u10 for exp). Full bordeaux marks configurations that are faster and meet the 10−6 target; light bordeaux marks faster configurations whose error exceeds the target.
• General purpose CPUs (high K, low N ): On architectures where polynomial coefficients can be stored in registers and evaluated with low latency, it is preferable to move along the frontier toward higher K and fewer iterations, since the extra coefficients are essentially free while each iteration adds latency (see the throughput ratios in Tables 3 and 4). • FPGAs and embedded platforms (low K, high N ): On architectures with limited coefficient storage, register availability, or multiplier resources, it is preferable to move toward lower K, compensating with a higher iteration count N to stay on the frontier. This trades coefficient storage for repeated arithmetic, reducing the memory footprint at the cost of additional iterations (see the throughput ratios in Tables 3 and 4). This trade-off is conceptually similar to CORDIC, exchanging a complex initial approximation for additional iterative steps to maintain accuracy. The underlying mechanisms differ fundamentally. CORDIC relies on precomputed look-up tables, while the present approach reconstructs the function value via a contractive rational recurrence in the residual coordinate, requiring no tables. This contractive property attenuates rounding errors from one step to the next, making the attainable
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
coordinate
N
K
max ulp
vs u35‡
vs u10‡
residual ρ residual ρ residual ρ residual ρ residual ρ
0 1 2 3 4
8 8 8 8 8
4.9 × 106 1.5 × 103 1.50 1.50 1.50
1.55× 0.90× 0.50× 0.34× 0.25×
2.41× 1.40× 0.78× 0.52× 0.39×
residual ρ residual ρ residual ρ residual ρ residual ρ
0 1 2 3 4
12 12 12 12 12
4.6 × 103 1.50 1.50 1.50 1.50
1.27× 0.80× 0.49× 0.33× 0.25×
1.97× 1.24× 0.75× 0.51× 0.39×
20
Table 5. Preliminary simulation of the logarithm in the contractive residual coordinate ρ. Unlike the value-coordinate formulation, the error stabilizes at a flat floor (≈ 1.5 ulp) as the iteration depth N increases, demonstrating damped error propagation. The rows far from the floor (N = 0; N = 1 at K = 8) reflect that the residual approach relies on the iteration: without descent, a short seed cannot cover the reduced interval. ‡ Throughput ratios are unoptimized (evaluating standard square roots in plain form) and are reported only to indicate order of magnitude relative to SLEEF.
accuracy essentially independent of the iteration count N . In contrast, a direct reconstruction would amplify rounding errors, causing accuracy to degrade with N (see Table 1). 6.4. Discussion. The preceding sections analyzed the numerical accuracy and SIMD performance of the proposed fixed-point algorithms. We conclude by summarizing the principal algorithmic trade-offs, discussing the current limitations, and outlining directions for future work. 6.4.1. Performance Trade-offs. As discussed in §6.3.2, the parameters (K, N ) define a family of implementations that exchange polynomial complexity for iteration depth while maintaining essentially the same accuracy. For the real exponential, configurations ranging from (K, N ) = (12, 0) to (7, 3) all achieve a maximum error of approximately 0.65–0.67 ulp. The choice of (K, N ) may therefore be adapted to the characteristics of the target architecture without sacrificing numerical accuracy. On the CPUs considered here, the throughput optimum lies at the low-N end of the frontier (N = 0, reaching 1.08× the SLEEF u10 throughput). Configurations with smaller polynomial seeds and deeper iteration become attractive on architectures where polynomial evaluation is comparatively expensive, such as extended-precision arithmetic or memory-constrained hardware. 6.4.2. Future Work: A Residual Formulation for the Logarithm. Natural engineering extensions include implementations for additional SIMD instruction sets, such as AVX-512 and ARM Neon, together with branchless handling of IEEE-754 exceptional values. A more significant algorithmic direction concerns the logarithm. Throughout this work, the logarithm is reconstructed from its macroscopic value. As discussed in §6.2, this reconstruction amplifies accumulated rounding errors by a factor of 2N . This amplification is not an intrinsic property of the logarithm, but rather a consequence of the chosen coordinate. As for the exponential, one may instead introduce the residual coordinate log(1 + u) − u + 12 u2 , (41) u3 √ so that log(1 + u) = u − 12 u2 + u3 ρ(u). Setting v := 1 + u − 1, one gets log(1 + u) = 2 log(1 + v), which yields 2ρ(v) + 2 + 12 v ρ(u) = . (42) (2 + v)3 ρ(u) =
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
21
The Taylor expansion of log(1 + u) around zero gives ρ(u) = 13 + O(u), showing that the residual remains bounded as u → 0 and admits the continuous extension ρ(0) = 1/3. The important stability property, however, is the propagation of perturbations through the fixed-point iteration. Since the 2 1 contribution of an error in ρ(v) to the updated value of ρ(u) is multiplied by (2+v) 3 = 4 + O(v), small perturbations are reduced by approximately a factor of four when u is close to zero. Thus, the residual formulation damps rounding errors during the iteration rather than amplifying them, in contrast with the macroscopic formulation, where reconstruction errors grow proportionally to 2N . A rigorous convergence analysis of this operator, for example in an appropriate weighted Banach space, remains the subject of future work. Nevertheless, preliminary numerical experiments exhibit the expected behavior. As shown in Table 5, increasing the iteration depth rapidly reduces the approximation error, after which the maximum error stabilizes at approximately 1.5 ulp, with no evidence of the exponential error growth observed in the value-coordinate formulation. These preliminary results indicate that the residual formulation restores the algorithmic symmetry between the exponential and logarithm. Achieving full u10 accuracy will likely require carrying the final reconstruction, E log 2 + log m, in extended precision. 6.4.3. Concluding Remarks. The numerical stability and computational efficiency of the proposed algorithms follow directly from the interaction between their fixed-point structure and finite-precision arithmetic. SIMD architectures are particularly well suited to these methods because they replace table lookups with arithmetic operations that can be efficiently vectorized. The choice of coordinate determines the numerical behavior. For the exponential, performing the iteration in a contractive residual coordinate suppresses the propagation of rounding errors. For the logarithm, reconstructing the function from its macroscopic value requires explicit stabilization, whereas performing the iteration directly in the residual coordinate recovers the same contractive mechanism as for the exponential, providing a unified framework for both elementary functions. The main contribution of this work is the identification of the residual coordinate as the fundamental mechanism that eliminates the 2N amplification of rounding errors. This is a structural property of the underlying algorithm rather than of a particular implementation, and we expect it to remain relevant across different hardware platforms and software environments. Fully exploiting its practical potential will require the extensive architecture-specific optimization, engineering effort, and benchmarking that have gone into the development of state-of-the-art mathematical libraries. We hope that the present work provides a sound theoretical foundation and a practical point of departure for researchers and developers working on high-performance mathematical software, both in academia and in industry, and that it encourages further interaction between the mathematical analysis, numerical analysis, and high-performance computing communities in the development of faster, more reliable, and more energy-efficient elementary function libraries. Declarations Author Contributions: All authors contributed equally to this work. Ethical Approval: Not applicable. This study does not involve human participants or animals. Conflict of Interest: The authors declare that they have no conflicts of interest. Data and Code Availability: No datasets were generated or analyzed during the current study. The C++ source code and benchmarking kernels evaluated in this work are currently available from the corresponding author upon reasonable request. The authors plan to clean, document, and host the complete codebase in a publicly accessible GitHub repository before final publication. Declaration of Generative AI and AI-Assisted Technologies: During the preparation of this manuscript, the authors used Grammarly and Gemini to check spelling and improve the fluency of selected portions of the text. The manuscript was initially prepared using TeXmacs and subsequently exported to LaTeX. Gemini was used to identify and correct translation or conversion errors arising during this process, to assist with LaTeX formatting, and to generate LaTeX code to improve the formatting of tables that were not converted correctly during the export.
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
22
In addition, Claude was used to assist in implementing and testing the algorithms described in the paper. No generative AI tools were used in the development of the mathematical ideas, the derivation of the theoretical results, or the proofs. The authors are fully responsible for the content of the manuscript and for the accuracy, validity, and integrity of all results presented.
Acknowledgements G.D.F. is a member of GNAMPA–INdAM. He gratefully acknowledges partial financial support from the GNAMPA Project CUP_E53C25002010001, and from the University of Naples Federico II through the FRA Project-B “VarMoCry” on Variational Analysis and Modeling of Liquid Crystals. Further support is acknowledged from the Italian Ministry of University and Research through the PRIN 2022 project Variational Analysis of Complex Systems in Material Science, Physics and Biology (No. 2022HKBF5C).
References János Aczél. Lectures on Functional Equations and Their Applications. Academic Press, New York, 1966. Lars V. Ahlfors. Complex Analysis. McGraw-Hill, New York, 3rd edition, 1979. Awad H. Al-Mohy and Nicholas J. Higham. A new scaling and squaring algorithm for the matrix exponential. SIAM Journal on Matrix Analysis and Applications, 31(3):970–989, 2009. [Apo67] Tom M. Apostol. Calculus, Vol. 1: One-Variable Calculus, with an Introduction to Linear Algebra. Wiley, New York, 2nd edition, 1967. [BCFG26] Enrico Babilio, Claudia Capone, Alberto Fiorenza, and Filomena Galizia. How to define sine and cosine as functions over reals rigorously and with minimal prerequisites. Results in Mathematics, 81(27), 2026. [GD03] Andrzej Granas and James Dugundji. Fixed Point Theory. Springer-Verlag, New York, 2003. [Har08] G. H. Hardy. A Course of Pure Mathematics. Cambridge University Press, Cambridge, 1908. [Hig05] Nicholas J. Higham. The scaling and squaring method for the matrix exponential revisited. SIAM Journal on Matrix Analysis and Applications, 26(4):1179–1193, 2005. [Hig08] Nicholas J. Higham. Functions of Matrices: Theory and Computation. SIAM, Philadelphia, 2008. [Kuc09] Marek Kuczma. An Introduction to the Theory of Functional Equations and Inequalities. Birkhäuser, Basel, 2nd edition, 2009. [Lan97] Serge Lang. Undergraduate Analysis. Springer, New York, 2nd edition, 1997. [MBdD+ 18] Jean-Michel Muller, Nicolas Brunie, Florent de Dinechin, Claude-Pierre Jeannerod, Mioara Joldes, Vincent Lefèvre, Guillaume Melquiond, Nathalie Revol, and Serge Torres. Handbook of Floating-Point Arithmetic. Birkhäuser, Cham, 2nd edition, 2018. Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of a matrix, twenty-five [ML03] years later. SIAM Review, 45(1):3–49, 2003. [Pro77] Giovanni Prodi. Analisi Matematica. Programma di matematica, fisica, elettronica. Bollati Boringhieri, Torino, 2nd edition, 1977. [Ric84] Fred Richman. A circular argument. The College Mathematics Journal, 15(4):362, 1984. [Rud76] Walter Rudin. Principles of Mathematical Analysis. McGraw-Hill, New York, 3rd edition, 1976. [SC17] Earl W. Swokowski and Jeffery A. Cole. Precalculus: Functions and Graphs. Cengage Learning, Boston, 12th edition, 2017. [Shi96] Georgi E. Shilov. Elementary Real and Complex Analysis. Dover Publications, Mineola, NY, 1996. [SP20] Naoki Shibata and Francesco Petrogalli. SLEEF: A portable vectorized library of C standard mathematical functions. IEEE Transactions on Parallel and Distributed Systems, 31(6):1316–1327, 2020. [Spi08] Michael Spivak. Calculus. Publish or Perish, Houston, 4th edition, 2008. [Ste16] James Stewart. Calculus. Cengage Learning, Boston, 8th edition, 2016. [Tre19] Lloyd N. Trefethen. Approximation Theory and Approximation Practice. SIAM, Philadelphia, extended edition, 2019. [Zei86] Eberhard Zeidler. Nonlinear Functional Analysis and Its Applications: I: Fixed-Point Theorems. SpringerVerlag, New York, 1986. [Acz66] [Ahl79] [AMH09]
A FIXED-POINT CONSTRUCTION OF THE ELEMENTARY TRANSCENDENTAL FUNCTIONS
23
François Alouges. 1 Université Paris Saclay, Université Paris Cité, ENS Paris-Saclay, CNRS, SSA, INSERM, Centre Borelli, 4 avenue des Sciences, Gif-sur-Yvette 91190, France. 2 Institut Universitaire de France. Email address: [email protected] Giovanni Di Fratta. 1 Dipartimento di Matematica e Applicazioni “R. Caccioppoli”, Università degli Studi di Napoli “Federico II”, Via Cintia, Complesso Monte S. Angelo, 80126 Naples, Italy Email address: [email protected] Alberto Fiorenza. 1 Dipartimento di Architettura, Università di Napoli, Via Monteoliveto 3, 80134 Napoli, Italy. 2 Istituto per le Applicazioni del Calcolo “Mauro Picone”, Sezione di Napoli, Consiglio Nazionale delle Ricerche, Via Pietro Castellino 111, 80131 Napoli, Italy Email address: [email protected] Renato Fiorenza. Accademia di Scienze Fisiche e Matematiche, Via Mezzocannone 8, 80134 Napoli, Italy Email address: [email protected]