<figure><img src="image_1.png" /></figure>

HAL
open science

# Approximation methods in Numerical Analysis

Charbel Tannous

► To cite this version:

Charbel Tannous. Approximation methods in Numerical Analysis. Master. Analyse Numérique, UBO Brest, France. 2016, pp.23. hal-05011363v2

HAL Id: hal-05011363

https://hal.science/hal-05011363v2

Submitted on 25 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Approximation methods in Numerical Analysis

C. Tannous

<figure><img src="image_3.png" /></figure>

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

Fast and accurate approximation methods are described for number estimation and function evaluation. Several approaches are reviewed: rational, continued fraction estimation and Complex integration techniques for number estimation. Regarding function evaluation, Polynomial estimation of functions, rational, continued fraction, Padé approximation of functions, least-square fitting, Levenberg-Marquardt method and Neural Network based methods.

PACS numbers: 02.60.-x, 02.60.Gf, 02.60.Ed, 84.35.+i

Keywords: Numerical methods (mathematics), Algorithms for Functional approximation, Least square approximation, Neural networks

Version: April 25, 2025

## Contents

I. Introduction 1

II. Estimation of numbers 2

A. Vieta and $\pi$ estimation 2

B. Continued fraction representation 3

C. Complex integration for number estimation 5

III. Polynomial representation of functions 7

A. Taylor expansion 7

B. Series expansion 8

C. Euler-MacLaurin Formula 9

D. Chebyshev Polynomial 9

IV. Rational Padé approximants 10

V. Continued fraction representation of functions 11

VI. Fast evaluation of functions 13

A. John Napier and the invention of logarithms 13

B. CORDIC algorithms 14

C. GPU architecture, AI, SFU and TPU 16

VII. Approximation of functions by regression 17

A. Least Squares Estimates 18

B. Bias, Variance and Standard Error of Parameter Estimates 18

C. Maximum Likelihood vs. Least Squares 19

D. Non-linear Least Squares: Levenberg-Marquardt method 20

References 23

## I. INTRODUCTION

Enabling a digital system to perform sophisticated calculations is based on developing a number representation, arithmetic and logic operations, functions and various operations on functions such as integration, derivation, series expansion... in order to tackle solutions of Ordinary Differential Equations (ODE), Partial Differential Equations (PDE), Finite Element Methods (FEM)...

In this work, fast and accurate numerical approximation is restricted to the framework of number estimation and evaluation of functions with several approaches such as rational estimation of numbers, Polynomial estimation of
functions, rational, continued fraction, Padé approximation of functions, linear least-square or non-linear least-square fitting without or with the use of Neural Networks.

## II. ESTIMATION OF NUMBERS

The best way to learn number estimation is to consult books or articles related to Fermi questions.

Fermi [1] had an amazing ability to answer with reasonable accuracy any question addressed to him, questions that would seem impossible to answer by an ordinary person.

Fermi used also to ask specific questions to any student who would approach him to do graduate work under his supervision.

The classic example of such questions that is attributed to him is to ask for the number of piano tuners in the city of Chicago.

Actually one might infer that this question can be split into a series of simpler questions which admit approximate answers leading eventually to an approximate answer to the original question:

1. Size of Chicago population

2. Number of families living in Chicago

3. Fraction of families possessing pianos at home

4. Average number of piano tunings per year and per family

5. Average number of piano tunings per year that a tuner can accomplish

While Fermi questions are enlightening and provide depth and fine training to students for doing high quality research, these notes are rather focused on the numerical aspect involved in estimation of mathematical constants starting with $\pi$ that triggered a lot of interest in the Greek world [2] and was responsible for the famous circle quadrature paradox.

Actually $\pi$ holds a special place in mathematics. It is an irrational number that continues infinitely without repetition.

People celebrate "Pi day" on March 14, because this date is written 3/14 in the USA, which is evocative of the decimal representation of $\pi = 3.14...$.

There are several other approximations: $\pi \approx 22/7$. as well as 355/113... and the multiple precision representation of $\pi$ leading to trillion of digits as was done by Emma Haruka Iwao, a Google employee, who spent four months working on a project in which she calculated $\pi$ to 31.4 trillion digits (i.e. 31,415,926,535,897 digits) and announced by Google on 'Pi day' March 14, 2019!

### A. Vieta and $\pi$ estimation

Vieta (François Viète) a French 16-th century mathematician decided to improve Archimedes who made the first evaluation of $\pi$ by estimating it from the area of a circle or perimeter of its circumference a problem known in Antiquity as "Squaring the circle". Ancient Greek mathematicians developed rules to make geometric construction of numbers that lead either to rational or special types of algebraic numbers that are roots of integer polynomials.

Given a polynomial:

$$a_n x^n + a_{n-1} x^{n-1} + a_{n-2} x^{n-2} + \dots + a_1 x + a_0 \quad (1)$$

if the $a_i$ coefficients are integers or rational, any root of the polynomial is considered as an algebraic number of degree $n$.

$\pi$ is a transcendental number i.e. is not root of any integer polynomial of any degree. This is why it took about 2000 years to understand the nature of $\pi$.

In fact, the value of $\pi$ has been known with remarkable accuracy for a long time [2, 3]. An Egyptian text dating to 1650 B.C. and known as the Rhind Papyrus (named after the Scottish Egyptologist Henry Rhind, who purchased it) contains the statement that a circle has the same area as a square whose side is 8/9 the diameter of the circle. Denoting the diameter by $d$, the statement translates into the equation $\pi(d/2)^2 = [(8/9)d]^2$, from which one gets, $\pi/4 = 64/81$, or $\pi = 256/81 = 3.16049$. This result comes within 0.6 percent of the true value of $\pi$ (3.14159, rounded to
<table><caption>**Table 1:** Iterations for the polygon areas inscribed (2<sup>i</sup>b<sub>i</sub>) within and inscribing (2<sup>i+1</sup>c<sub>i</sub>) a unit circle as depicted in Fig. 1 in double precision. At iteration 10 the values 3.1415914215112033 and 3.1415932178497314 are valid bounds for π.</caption><thead><tr><th>Iteration i</th><th>2<sup>i</sup>b<sub>i</sub></th><th>2<sup>i+1</sup>c<sub>i</sub></th></tr></thead><tbody><tr><td>0</td><td>2.00000000000000</td><td>4.00000000000000</td></tr><tr><td>1</td><td>2.82842712474619</td><td>3.313708543777466</td></tr><tr><td>2</td><td>3.0614674589207174</td><td>3.1825978755950928</td></tr><tr><td>3</td><td>3.121445152258051</td><td>3.1517248153686523</td></tr><tr><td>4</td><td>3.1365484905459375</td><td>3.144118309020996</td></tr><tr><td>5</td><td>3.1403311569547507</td><td>3.142223596572876</td></tr><tr><td>6</td><td>3.14127725093277</td><td>3.1417503356933594</td></tr><tr><td>7</td><td>3.1415138011442982</td><td>3.141632080078125</td></tr><tr><td>8</td><td>3.1415729403670887</td><td>3.1416025161743164</td></tr><tr><td>9</td><td>3.141587725277156</td><td>3.1415951251983643</td></tr><tr><td>10</td><td>3.1415914215112033</td><td>3.1415932178497314</td></tr></tbody></table>

five decimal places). This result is remarkably accurate for a text written almost four thousand years ago improving the round Biblical value $\pi \approx 3$ [2].

Archimedes [4] inscribed 96 triangles within the unit circle (cf. Fig. 1) to derive:

$$3\frac{1137}{8069} < \pi < 3\frac{1335}{9347} \quad \text{or} \quad 3.1409 < \pi < 3.1428 \qquad (2)$$

In 1593 Vieta considered two types of polygons, those inscribed into a circle and those containing a circle inscribed into them in order to bound $\pi$ rigorously (cf. Fig. 1).

<figure><img src="image_4.png" /><figcaption>**Fig. 1:** (left) Archimedes used the circumference of the triangles inscribed within the unit circle (in red). The approach can be extended to the circle surface as described by Rice [4]. (Center) Vieta estimated $\pi$ using polygons inscribed in a unit circle (in red) à la Archimedes and the reverse situation (Right) where the unit circle (in red) is inscribed within polygons [5]. The transitions from square AAAA to octagon ABAB... are shown in both cases.</figcaption></figure>

The iteration sequences giving the polygon areas inscribed within a unit circle [5, 6] are:

$$b_n = \sqrt{2} \sqrt{1 - \sqrt{1 - (b_{n-1}/2)^2}} \qquad (3)$$

whereas The iteration sequences giving the polygon areas inscribed within a unit circle [5, 6] are

$$c_n = \frac{4(\sqrt{1 + (c_{n-1}/2)^2} - 1)}{c_{n-1}} \qquad (4)$$

with $b_0 = c_0 = 2$ and $\lim_{n \to \infty} 2^n b_n \le \pi \le \lim_{n \to \infty} 2^{n+1} c_n$.

Details of convergence of both sequences are displayed in Table 1.

## B. Continued fraction representation

Every real number $x$ can be represented as a continued fraction where the $a_i$ are positive integers. With all of fractions having 1 in the numerator, the continued fraction can be compactly represented by specifying only the
integers:

$$x = [a_0; a_1, a_2, a_3, ...] \equiv a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + \frac{1}{a_4 + ...}}}} \quad (5)$$

Every rational number is represented by a finite sequence; every irrational number requires an infinite sequence. In fact, the decimal representation of a number does not have this property: 1/3 is rational but its decimal representation is infinite.

* The continued fraction (CF) expansion of Euler $e \approx 2.718281828459$ shows a regular pattern of 1, 1, 2n with $n > 1$ integer:

$$e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, ...] \quad (6)$$

A very interesting fact of CF is that if one wants to evaluate $[a_0; a_1, a_2, a_3, ...]$ it is possible to evaluate successively the ith value $[a_0; a_1, a_2, a_3, ...a_i]$ that converges to $p_i/q_i$ that can be evaluated recursively from

$$p_n = a_n p_{n-1} + p_{n-2}, \quad q_n = a_n q_{n-1} + q_{n-2} \quad (7)$$

starting from the initial conditions: $p_0 = a_0, q_0 = 1, p_1 = a_0a_1 + 1$ and $q_1 = a_1$.

In the Euler constant $e$ one takes: $a_{3i+1} = 2i, a_{3i} = a_{3i+2} = 1$ such that $\lim_{i \to \infty} p_i/q_i = e$.

* The CF expansion of $\sqrt{2} \approx 1.414213562373095$ contains a repeating sequence of 2:

$$\sqrt{2} = [1; 2, 2, 2, 2, ...] \quad (8)$$

This can be derived from Newton-Raphson recursion for finding a square root:

$$\lim_{n \to \infty} x_{n+1} = \lim_{n \to \infty} \frac{1}{2}(x_n + A/x_n) = \sqrt{A} \text{ implying that } \sqrt{2} = \lim_{n \to \infty} x_{n+1} = \lim_{n \to \infty} \frac{1}{2}(x_n + 2/x_n).$$

* The CF expansion of the golden mean $\phi \approx 1.61803398874989$ is a repeating sequence of 1:

$$\phi = [1; 1, 1, 1, 1, 1, ...] \quad (9)$$

This can be interpreted from the recursion: $x_{n+1} = (1 + 1/x_n)$ whose fixed point is found by taking $\lim_{n \to \infty}$ that yields $x^* = (1 + 1/x^*)$ or $[x^*]^2 = (1 + x^*)$. The solutions of this quadratic equation are $\phi, -1/\phi$.

* The continued fraction expansion of $\pi$ given that the first fifteen digits are

$$\pi \approx [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2, 2, 6, 3, 5, 1, 1, 6, 8, 1, 7, 1, 2, 3, 7, 1, 2, 1, 1, 12, 1, 1, 1, 3, 1, 1, 8, 1, 1, 2, 1, 6, 1, 1, 5, 2, 2, 3, 1, 2, 4, 4, 16, 1, 161, 45, 1, 22, 1, 2, 2, 1, ...]$$

This is interpreted from the fact the sequence [3; 7] corresponds to the fraction $3 + 1/7 = 22/7$. similarly, the sequence [3; 7, 15, 1] corresponds to the fraction $3 + 1/(7 + 1/(15 + 1/1)) = 355/113$.

However, if the decimal expansion is pushed farther, it is possible to build a series of different CF as given by Table 2.

As an example the following CF:

$$\pi = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2, 2, 6, 3, 5, 1, 1, 6, 8, 1, 7, 1, 2, 3, 7, 1, 2, 1, 1, 12, 1, 1, 1, 3, 1, 1, 8, 1, 1, 2, 1, 6, 1, 1, 5, 2, 2, 3, 1, 2, 4, 4, 16, 1, 161, 45, 1, 22, 1, 2, 2, 1, ...]$$

allows to yield a multiple precision of $\pi$ [7]:

$$\pi = 3.14159265358979323846264338327950288419716939937510582097494459230781...$$

This result can be derived, in fact, with the Arithmetic-Geometric mean algorithm [8, 9]:

$$\begin{align}
a_k &= \frac{a_{k-1} + b_{k-1}}{2} \quad , \quad b_k = \sqrt{a_{k-1} b_{k-1}}, \nonumber \\
c_k &= a_k^2 - b_k^2, \quad s_k = s_{k-1} - 2^k c_k, \quad p_k = \frac{2a_k^2}{s_k} \tag{11}
\end{align}$$
<table><caption>**Table 2:** Representation of $\pi$ as a truncated Continued fraction, Rational value, Decimal value and corresponding error.</caption><thead><tr><th>Truncated CF</th><th>Rational value</th><th>Decimal value</th><th>Error</th></tr></thead><tbody><tr><td>[3]</td><td>3</td><td>3.000000000000000</td><td>-1.41 ×10<sup>-1</sup></td></tr><tr><td>[3; 7]</td><td>22/7</td><td>3.14285714285714</td><td>1.23 ×10<sup>-3</sup></td></tr><tr><td>[3; 7, 15]</td><td>333/106</td><td>3.14150943396226</td><td>-8.32 ×10<sup>-5</sup></td></tr><tr><td>[3; 7, 15, 1]</td><td>355/113</td><td>3.14159292035398</td><td>2.66 ×10<sup>-7</sup></td></tr><tr><td>[3; 7, 15, 1, 292]</td><td>103993/33102</td><td>3.14159265301190</td><td>-5.77 ×10<sup>-10</sup></td></tr><tr><td>[3; 7, 15, 1, 292, 1]</td><td>104348/33215</td><td>3.14159265392142</td><td>3.31 ×10<sup>-10</sup></td></tr><tr><td>[3; 7, 15, 1, 292, 1, 1]</td><td>208341/66317</td><td>3.14159265346744</td><td>-1.22 ×10<sup>-10</sup></td></tr><tr><td>[3; 7, 15, 1, 292, 1, 1, 1]</td><td>312689/99532</td><td>3.14159265361894</td><td>2.91 ×10<sup>-11</sup></td></tr><tr><td>⋮</td><td>⋮</td><td>⋮</td><td>⋮</td></tr><tr><td>[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]</td><td>80143857/25510582</td><td>3.14159265358979</td><td>-4.44 ×10<sup>-16</sup></td></tr></tbody></table>

$\pi$ is obtained by taking $\lim_{k \to \infty} p_k$.

In fact, there is another CF representation given by two sets of coefficients $a_i, i = 0...\infty$ and $b_i, i = 1...\infty$:

$$a_0 + \frac{b_1}{a_1 + \frac{b_2}{a_2 + \frac{b_3}{a_3 + \frac{b_4}{a_4 + \dots}}}} \quad (12)$$

such that $\pi$ is represented with $a_0 = 3, a_i = 6 \quad \forall i > 0$, whereas $b_i = (2i - 1)^2 \quad \forall i$:

$$\pi = 3 + \frac{1}{6 + \frac{9}{6 + \frac{25}{6 + \frac{49}{6 + \dots}}}} \quad (13)$$

Using CF based on two sets of coefficients shows more regularity than with a single set. For example Euler $e$ is given by $a_0 = 2, a_i = i \quad \forall i > 0$, whereas $b_1 = 1, b_i = i - 1 \quad \forall i > 1$:

$$e = 2 + \frac{1}{1 + \frac{1}{2 + \frac{2}{3 + \frac{3}{4 + \dots}}}} \quad (14)$$

### C. Complex integration for number estimation

Analytic continuation allows to perform quadrature by complex integration [10] leading to another type of approximation.

Consider the integral:

$$I_n = \int_C (z - z_0)^n dz, \quad (15)$$

where $C$ is a closed contour surrounding $z = z_0$ and $n$ is an integer. Since $(z - z_0)^n$ is analytic, except possibly at $z_0$, the contour may be deformed into a circle centered at $z_0$ without affecting $I_n$. Then one may write $z - z_0 = Re^{i\theta}, dz = iRe^{i\theta}d\theta$ to write:

$$I_n = iR^{n+1} \int_0^{2\pi} e^{i(n+1)\theta} d\theta = \begin{cases} 0, & n \neq -1 \\ 2i\pi, & n = -1 \end{cases} \quad (16)$$

Now consider a function, $f(z)$, which is analytic in a region except at the point $z_0$ in the region. Then performing Laurent expansion [10] about $z_0$ gives:

$$f(z) = \sum_{n=-\infty}^{\infty} A_n(z - z_0)^n. \quad (17)$$
Take a contour C around $z_0$:

$$ \frac{1}{2\pi i} \int_C f(z) dz = \frac{1}{2\pi i} \int_C \sum_{n=-\infty}^{\infty} A_n(z-z_0)^n dz = A_{-1} \quad (18) $$

Thus the coefficient of $1/(z-z_0)$ in the Laurent series is given by $A_{-1} = \frac{1}{2\pi i} \int_C f(z) dz$. This coefficient is called the residue of $f(z)$ at $z_0$. Notice that the residue is zero if $f(z)$ is analytic at $z_0$, or if the coefficient $A_{-1}$ is zero (even if $z_0$ is a pole or isolated essential singularity).

An important and useful notion is the residue theorem [10]. Consider contour C in a region where $f(z)$ is analytic except at isolated singularities (poles or essential singularities) as in Fig. 2.

<figure><img src="image_5.png" /><figcaption>Fig.2: Contours to illustrate residue theorem. Singularities are at a, b, and c.</figcaption></figure>

In order to determine $\int_C f(z) dz$, write it in terms of a sum of integrals around each singularity:

$$ \int_C f(z) dz = \int_{C_a} f(z) dz + \int_{C_b} f(z) dz + \dots = 2i\pi(a_{-1} + b_{-1} + c_{-1} + \dots) = 2i\pi \sum_{a,b,c,\dots} R_e, \quad (19) $$

where $a_{-1}, b_{-1}, c_{-1}, \dots$ are the residues at a, b, c, ..., respectively, and $\sum R_e$ is the sum of the residues of $f(z)$ interior to the contour C.

The computation of residues is thus often an important part of evaluating integrals. At a simple pole, the Laurent series is

$$ f(z) = \frac{A_{-1}}{z-z_0} + \sum_{n=0}^{\infty} A_n(z-z_0)^n, \quad A_{-1} = \lim_{z \to z_0} [(z-z_0)f(z)] $$

For a pole of order $m$,

$$ f(z) = \frac{A_{-m}}{(z-z_0)^m} + \frac{A_{-m+1}}{(z-z_0)^{m-1}} + \dots + \frac{A_{-1}}{z-z_0} + \sum_{n=0}^{\infty} A_n(z-z_0)^n. $$

Multiplying both sides by $(z-z_0)^m$, yields:

$$ (z-z_0)^m f(z) = A_{-m} + A_{-m+1}(z-z_0) + \dots + A_{-1}(z-z_0)^{m-1} + \sum_{n=0}^{\infty} A_n(z-z_0)^{n+m}. $$

Now differentiate $m-1$ times and evaluate at $z = z_0$:

$$ \left. \frac{d^{m-1}[(z-z_0)^m f(z)]}{dz^{m-1}} \right|_{z=z_0} = (m-1)! A_{-1}, \quad \text{and} \quad A_{-1} = \frac{1}{(m-1)!} \left. \frac{d^{m-1}[(z-z_0)^m f(z)]}{dz^{m-1}} \right|_{z=z_0}. \quad (20) $$

However there are fast methods for residue evaluation as described in Ref. [10].

### Fast calculation of residues

For a simple pole $z_0$, the residue [10] is $\lim_{z \to z_0} [(z-z_0)f(z)]$.
<figure><img src="image_6.png" /><figcaption>**Fig.3:** Circle of radius 1 as contour C for complex integration of $(1/2i) \oint_C dz \frac{z^{6+1}}{z^3(2z-1)(z-2)}$. Inside C there are two real poles: a pole of order 3 at $z=0$ and a simple pole at $z=1/2$.</figcaption></figure>

For an order $m > 1$ pole, supposing the integrand $f(z) = N(z)/D(z)$, then the residue [10] is $N(z_0)/D'(z_0)$, provided $D'(z_0) \neq 0$ without differentiating $f(z)$, $m-1$ times or performing Laurent expansion around the $m > 1$ order pole. In the case $D'(z_0) = 0$, the residue is found by Laurent expansion [10].

As an example, consider evaluating the trigonometric integral [10]:

$$ \int_0^{2\pi} d\theta \frac{\cos 3\theta}{5 - 4 \cos \theta} \qquad (21) $$

Using $z = e^{i\theta}$ this is rewritten with $\cos \theta = (z+1/z)/2$, $\cos 3\theta = (z^3+1/z^3)/2$ and $dz = izd\theta$, as $\oint_C (dz/iz) \frac{(z^3+1/z^3)/2}{5-4z(z+1/z)/2}$ where C is the contour shown in Fig. 3.

The integral is transformed into $(-1/2i) \oint_C dz f(z)$ with $f(z) = \frac{z^{6+1}}{z^3(2z-1)(z-2)}$. Inside the contour, the integrand has a pole of order 3 at $z=0$ and a simple pole at $z=1/2$.

The integrand $f(z) = \frac{z^{6+1}}{z^3(2z-1)(z-2)}$ has a denominator $D(z) = z^3(2z-1)(z-2)$ thus $D'(z=0) = 0$ requiring a Laurent expansion of around $z=0$ that yields:

$$ \frac{z^6 + 1}{z^3(2z-1)(z-2)} = \frac{z^6 + 1}{z^3} \frac{1}{(2z-1)(z-2)} \approx \frac{1}{z^3} \left( 1 + \frac{5z}{4} + \frac{21z^2}{8} + \frac{85z^3}{16} + \dots \right) \qquad (22) $$

The residue is 21/8, the 1/z coefficient of Laurent expansion. The residue at the simple pole $z=1/2$ is -65/24. Thus the integral is evaluated as:

$$ \int_0^{2\pi} d\theta \frac{\cos 3\theta}{5 - 4 \cos \theta} = (-1/2i)2i\pi(21/8 - 65/24) = \pi/12 \qquad (23) $$

Thus another definition of $\pi$ is obtained from the result: $\pi = 12 \int_0^{2\pi} d\theta \frac{\cos 3\theta}{5-4 \cos \theta}$.

### III. POLYNOMIAL REPRESENTATION OF FUNCTIONS

Functions may be represented around a given point $x_0$ by polynomials following Taylor, Euler-MacLaurin, perturbation or even asymptotic expansions when $x_0 \to \infty$ [11, 12]. However, there is also need to approximate a function $f(x)$ from its series representation valid $\forall x$ as described in next section.

#### A. Taylor expansion

Suppose that $f: \mathbf{R} \to \mathbf{R}$ is a sufficiently differentiable function. The Taylor expansion of the function $f$ at a point $c \in \mathbf{R}$ is given by

$$ f(x) = p_n(x) + r_n(x), $$
where

$$p_n(x) = f(c) + f'(c)(x-c) + \frac{f'(c)}{2!}(x-c)^2 + \dots + \frac{f^{(n)}(c)}{n!}(x-c)^n, \quad r_n(x) = \frac{f^{(n+1)}(\xi(x,c))}{(n+1)!}(x-c)^{n+1} \quad (24)$$

for some $\xi(x, c)$ between $x$ and $c$. If $c \in [a, b]$ and $|f^{(n+1)}(x)| < M$ for $[a \le x \le b$, to yield:

$$\sup_{a \le x \le b} |f(x) - p_n(x)| \le \frac{M}{(n-1)!} (b-a)^{n+1}.$$

This is an uniform approximation of $f$ by a polynomial $p_n$. If $M$ is small and $b-a \le 1$, then $p$ can be a good uniform approximation of $f$. So, instead of evaluating $f(x)$, it may be simpler and sufficiently accurate to evaluate $p_n(x)$. Unfortunately, in practice, it may be very hard to compute the derivatives of $f$ and to find an upper bound on $|f^{n+1}|$. Moreover, the interval $[a : b]$ may be of length greater than 1 and so $(b-a)^{n+1} \to \infty$ as $n \to \infty$.

The fundamental Stone-Weierstrass theorem is the following:

**Theorem 1.** Given a continuous function $f : [a : b] \to \mathcal{C}$, there exists a sequence of polynomials $\{p_n\}_{n=0}^{\infty}$ over $\mathcal{C}$ such that

$$\max_{a \le x \le b} |f(x) - p_n(x)| \to 0 \quad \text{as} \quad n \to \infty.$$

If $f : [a : b] \to \mathcal{R}$, the polynomials can be assumed to be over $\mathcal{R}$.

Basically, the theorem states that for any continuous function $f : [a : b] \to \mathcal{C}$, it is possible to find a sequence of polynomials converging to $f$ uniformly over $[a : b]$.

## B. Series expansion

The difference between local expansions around some chosen point $x_0$ and series is that the latter is valid $\forall x$. This is how Euler while searching for a function [3] whose derivative and integral are same, wrote the exponential series as:

$$S(x) = \sum_{k=0}^{\infty} \frac{x^k}{k!}, \quad \forall x \quad (25)$$

The equivalence of the series with $e^x$ was done by evaluating $S(x = 1)$ that yielded $e \approx 2.718281828459$ and allowed to make the identification as $S(x) \equiv e^x$ after checking the validity of the product rule $S(x)S(y) \equiv e^x e^y = e^{x+y}$.

Actually Newton [3] was the first to derive in 1665:

$$1 + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + \frac{1}{4!} + \frac{1}{5!} + \frac{1}{6!} + \frac{1}{7!} + \frac{1}{8!} + \frac{1}{9!} + \frac{1}{10!} + \dots \approx 2.718281828459 \equiv e \quad (26)$$

It can be obtained from the binomial expansion of $(1 + 1/n)^n$ by taking $n \to \infty$ to obtain $\lim_{n \to \infty} (1 + 1/n)^n = e$. It converges very quickly, due to fast increasing values of the denominator factorials. For example, the sum of the first eleven terms (ending with $1/10!$) is 2.718281801. The true value, rounded to nine decimal places, is 2.718281828.

This result leads to the trigonometric series [3]:

$$\cos(x) = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k)!} x^{2k}, \quad \forall x, \quad \sin(x) = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)!} x^{2k+1}, \quad \forall x \quad (27)$$

and to the hyperbolic functions [3]:

$$\cosh(x) = \sum_{k=0}^{\infty} \frac{1}{(2k)!} x^{2k}, \quad \forall x, \quad \sinh(x) = \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} x^{2k+1}, \quad \forall x \quad (28)$$

which implied that $e^x = \cosh(x) + \sinh(x)$.

Analytic continuation [10] leads to Euler identity: $e^{ix} = \cos(x) + i \sin(x)$, a transition from real to complex $\cos(x) = \cosh(ix), i \sin(x) = \sinh(ix)$ and the most beautiful formula of Mathematics: $e^{i\pi} = -1$.
### C. Euler-MacLaurin Formula

Euler-MacLaurin expansion [13] is the Taylor expansion of a function eq. 24 at point $c = 0$.

In contrast, Euler-MacLaurin Formula is the most important method for solving finite and infinite series representing functions. It gives a relationship between a sum and its corresponding integral.

Let the first $2n$ derivatives of $f(x)$ be continuous over an interval $[a : b]$ divided into equal parts $h = \frac{(b-a)}{n}$.

Then Euler-MacLaurin Formula [13] is:

$$ \sum_{k=0}^{n-1} f(a + kh) = \frac{1}{h} \int_{a}^{b} f(x) dx - \frac{h}{2} (f(b) - f(a)) - \sum_{k=1}^{m-1} B_{2k} h^{2k} \frac{[f^{(2k-1)}(b) - f^{(2k-1)}(a)]}{(2k)!} \\ + \frac{h^{2m+1}}{(2m)!} \int_{0}^{1} dx [B_{2m}(x) - B_{2m}] \sum_{k=0}^{n-1} f^{(2m)}(a + hx + hk) \quad (29) $$

where $m$ is an arbitrary number, $B_{2k}$ are Bernoulli numbers whereas $B_{2m}(x)$ are Bernoulli polynomials defined below.

The $B_m$ are Bernoulli numbers given by the expansion:

$$ \frac{x}{e^x - 1} = \sum_{m=0}^{\infty} \frac{B_m x^m}{m!} = 1 - \frac{x}{2} + \frac{x^2}{12} - \frac{x^4}{720} + \frac{x^6}{30240} + \dots \quad (30) $$

such that: $B_0 = 1, B_1 = 0, B_2 = \frac{1}{6}, B_3 = 0, B_4 = -\frac{1}{30}, B_5 = 0, B_6 = \frac{1}{42}, \dots$

The $B_m(x)$ are Bernoulli polynomials given by the expansion:

$$ \frac{te^{xt}}{e^t - 1} = \sum_{n=0}^{\infty} \frac{B_n(x)t^n}{n!} \quad (31) $$

For an illustration, consider the sum: $S(x) = \sum_{k=1}^{x} k^3$, Euler-MacLaurin Formula yields:

$$ S(x) = \int_{1}^{x} du \, u^{3} + \frac{1}{2}(x^{3} + 1) + \frac{1}{2}B_{2}(3x^{2} - 3) + \frac{1}{4!}B_{4}(6 - 6) \quad (32) $$

Then $S(x) = \frac{1}{4}x^2(x+1)^2 = \frac{1}{4}(x^4 + 2x^3 + x^2)$.

### D. Chebyshev Polynomial

Chebyshev polynomial $T_n(x)$ are orthogonal over interval $[-1, 1]$ considering a weight function $w(x) = 1/\sqrt{1-x^2}$. There is an orthogonal set of functions $\{\phi_0, \phi_1, \dots, \phi_n\}$, over $[a : b]$ with respect to the weight function $w$ such that:

$$ \int_{b}^{a} w(x) \phi_{j}(x) \phi_{k}(x) dx = \begin{cases} 0, & j \neq k \\ \alpha_{k} > 0, & j = k \end{cases} $$

If in addition $\alpha = 1$ for each $k = 0, \dots, n$, the set is considered orthonormal.

Chebyshev polynomial $T_n(x)$ are defined as:

$$ T_0(x) = 1, \quad T_1(x) = x, \quad T_2(x) = 2x^2 - 1, \quad T_3(x) = 4x^3 - 3x, \quad T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x) \quad (33) $$

$T_n(x)$ can also be defined as

$$ T_n(x) = \cos[n(\cos^{-1}(x))] $$

Substituting $\theta = (\cos^{-1}(x))$ changes the expression to

$$ T_n(\theta(x)) = \cos n\theta \quad \theta \in [0, \pi] $$

The recursion relation is derived by noting that:

$$ T_{n+1}(x) = \cos(n\theta) \cos\theta - \sin(n\theta) \sin\theta, \quad T_{n-1}(x) = \cos(n\theta) \cos\theta + \sin(n\theta) \sin\theta \quad (34) $$
so that $T_{n+1}(x) = 2 \cos(n\theta) \cos\theta - T_{n-1}(x)$.

Returning to the variable $x$ gives $T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x)$.

$$ \int_{-1}^{1} \frac{T_n(x)T_m(x)}{\sqrt{1-x^2}} dx = \begin{cases} 0, & m \neq n \\ \pi, & n = m = 0 \\ \frac{\pi}{2}, & j = k \neq 0 \end{cases} $$

It is possible to express $x$ in terms of $T_n$:

$$ 1 = T_0, \quad x = T_1, x^2 = \frac{1}{2}(T_0 + T_2), \quad x^2 = \frac{1}{4}(3T_1 + T_3)\dots \quad (35) $$

If $P_n(x)$ is the interpolation polynomial of at most degree $n$, with nodes that are roots of $T_{n+1}(x)$ such that $P_n(x) = \sum_{i=0}^{n} c_i T_i(x)$ where $c_i$ are evaluated as follows:

$$ c_i = \frac{\int_{-1}^{1} dx f(x) T_i(x) w(x)}{\int_{-1}^{1} dx [T_i(x)]^2 w(x)} \quad (36) $$

The Chebyshev approximation is estimated as:

$$ \max_{x \in [-1,1]} |f(x) - P_n(x)| \le \frac{1}{2^n(n+1)!} \max_{x \in [-1,1]} |f^{n+1}(x)| \quad f \in C^{n+1}[-1,1] $$

As an example, consider the $\tan(x)$ function that is Taylor expanded to order 10 about $x = 0$:

$$ \tan(x) \approx x + \frac{1}{3}x^3 + \frac{2}{15}x^5 + \frac{7}{315}x^7 + \frac{62}{2835}x^9 + O(x^{10}) \quad (37) $$

Its Chebyshev approximation is given by [7]:

$$ \tan(x) \approx 1.0040x + 0.325610x^3 + 0.173344x^5 - 0.026645x^7 + 0.0846426x^9 + O(x^{10}) \quad (38) $$

In the next section, comparison is made between Taylor, Chebyshev expansions, Padé rational approximation and the exact function $\tan(x)$.

## IV. RATIONAL PADÉ APPROXIMANTS

A rational representation of a function is expected to be superior to a polynomial expansion. Let $f$ be a function holomorphic in a neighborhood of 0 and let $f(z) = \sum_{v=0}^{\infty} a_v z^v$ its Taylor series. A Padé approximant $[p/q]_f$ of $f$, $p, q \in \{0, 1, 2, \dots\}$, is a rational function of the form

$$ \frac{\sum_{v=0}^{p} n_v z^v}{\sum_{v=0}^{q} d_v z^v}, \quad d_0 = 1. $$

such that its Taylor series $\sum_{v=0}^{\infty} b_v z^v$ coincides with $\sum_{v=0}^{\infty} a_v z^v$ up to the first $p+q+1$ terms; that is $b_v = a_v$ for $v = 0, \dots, p+q$.

In case of $q = 0$ there exists always a unique Padé approximant [14, 15] of $f$ and $[p/q]_f(z) = S_p(z)$, where $S_p(z) = \sum_{v=0}^{p} a_v z^v$. For $q \ge 1$ it is true that there exists a unique Padé approximant of $f$, if and only if the following determinant is not zero:

$$ \det \begin{vmatrix} a_{p-q+1} & a_{p-q+2} & \cdots & a_p \\ a_{p-q+2} & a_{p-q+3} & \cdots & a_{p+1} \\ \vdots & \vdots & & \vdots \\ a_p & a_{p+1} & \cdots & a_{p+q-1} \end{vmatrix} \neq 0, \quad a_i = 0, \text{ for } i < 0. \quad (*) $$
Then $f \in D_{p,q}$ is obtained.

If $f \in D_{p,q}$, then $[p/q]_f (q \ge 1)$ is given by the explicit formula:

$$ [p/q]_f = \frac{\det \begin{vmatrix} z^q S_{p-q}(z) & z^{q-1} S_{p-q+1} & \cdots & S_p(z) \\ a_{p-q+1} & a_{p-q+2} & \cdots & a_{p+1} \\ \vdots & \vdots & & \vdots \\ a_p & a_{p+1} & \cdots & a_{p+q} \end{vmatrix}}{\det \begin{vmatrix} z^q & z^{q-1} & \cdots & 1 \\ a_{p-q+1} & a_{p-q+2} & \cdots & a_{p+1} \\ \vdots & \vdots & & \vdots \\ a_p & a_{p+1} & \cdots & a_{p+q} \end{vmatrix}}, $$

$$ \text{with } S_k(z) = \begin{cases} \sum_{v=0}^{k} a_v z^v, & k \ge 0 \\ 0, & k < 0. \end{cases} $$

For example, the exponential $\exp(x)$ function when Taylor expanded around $x = 0$ to degree 8:

$$ e^x \approx \sum_{k=0}^{8} \frac{x^k}{k!} = 1 + x + \frac{1}{2}x^2 + \frac{1}{6}x^3 + \frac{1}{24}x^4 + \frac{1}{120}x^5 + \frac{1}{720}x^6 + \frac{1}{5040}x^7 + \frac{1}{40320}x^8 + \dots \quad (39) $$

and one makes it coincide with a [4/4] Padé approximant [16] given by:

$$ e^x \approx \frac{1680 + 840x + 180x^2 + 20x^3 + x^4}{1680 - 840x + 180x^2 - 20x^3 + x^4} = \frac{1 + \frac{1}{2}x + \frac{9}{84}x^2 + \frac{1}{84}x^3 + \frac{1}{1680}x^4}{1 - \frac{1}{2}x + \frac{9}{84}x^2 - \frac{1}{84}x^3 + \frac{1}{1680}x^4} \quad (40) $$

Expanding $\exp(x)$ function to degree 10 around $x = 0$:

$$ e^x \approx \sum_{k=0}^{10} \frac{x^k}{k!} \quad (41) $$

yields [5/5] Padé approximant [16] given by:

$$ e^x \approx \frac{1 + \frac{1}{2}x + \frac{1}{9}x^2 + \frac{1}{72}x^3 + \frac{1}{1008}x^4 + \frac{1}{30240}x^5}{1 - \frac{1}{2}x + \frac{1}{9}x^2 - \frac{1}{72}x^3 + \frac{1}{1008}x^4 - \frac{1}{30240}x^5} \quad (42) $$

then it is possible to show in Fig. 4, a comparison between [4/4], [5/5] approximants and the exact exponential over the interval [0:6]. It is remarkable that [5/5] Padé approximant [7] stays close to the exact $\exp(x)$ despite the fact [0:6] range is far from the origin $x = 0$. In contrast, the [4/4] approximant starts moving away from $\exp(x)$ above $x \sim 4.5$.

Moving on to the $\tan(x)$ approximation, Taylor expansion about $x = 0$ to order 10 provides Padé [5/5] rational approximant as:

$$ \tan(x) \approx \frac{x^5 - 103x^3 + 945x}{15x^4 - 420x^2 + 945} \quad (43) $$

Note that Padé [5/4] rational approximant is same as Padé [5/5]. Comparison between exact $\tan(x)$, Taylor and Chebyshev expansions and Padé [5/5] approximant is shown in Fig. 5 displaying the superior approximation of Padé [5/5] representation over the interval [0:1.5].

## V. CONTINUED FRACTION REPRESENTATION OF FUNCTIONS

A function represented by a continued fraction [17, 18] can be defined from:
<figure><img src="image_7.png" /><figcaption>Fig.4: Comparison of [4/4], [5/5] Padé approximants and the exact exponential over the interval [0:6]. Discrepancies start showing above x ~ 4.5 for [4/4] approximant whereas [5/5] approximant stays close to the exact exp(x) function.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig.5: Exact tan(x), Taylor, Chebyshev polynomial and Padé [5/5] approximations revealing the best approximation as Padé [5/5] over the interval [0:1.5]. Exact tan(x) and Padé [5/5] are indistinguishable over [0:1.5].</figcaption></figure>

Definition 1. A sequence {($a_j, b_j$), $f_j$}, where $b_0, a_j, b_j, f_j \in \mathbb{R} \forall j \in \mathbb{N}$ is called a continued fraction if

$$f_n = b_0 + \frac{a_1}{b_1 + \frac{a_2}{b_2 + \frac{a_3}{b_3 + \dots + \frac{a_n}{b_n}}}}, \quad \forall n \in \mathbb{N}.$$

The number $f_n$ is called the n-th convergent of the continued fraction (CF). A CF {($a_j, b_j$), $f_j$} converges if there exists finite or infinite limit

$$\lim_{n \to \infty} f_n = f_*.$$
An interesting fact is the existence of inductive procedure that enables constructing the n-th convergents [17].
Namely, for a CF $\{(a_j, b_j), f_j\}$, define:

$$
\begin{align}
p_n &= b_n p_{n-1} + a_n p_{n-2}, \nonumber \\
q_n &= b_n q_{n-1} + a_n q_{n-2}, \tag{44} \\
p_0 &= b_0, \quad p_{-1} = 1, \quad q_0 = 1, \quad q_{-1} = 0. \nonumber
\end{align}
$$

Then

$$
f_n = \frac{p_n}{q_n}. \tag{45}
$$

**Examples of functions [13]:**

$$
\bullet \ln(1+z) = \frac{z}{1+z} \frac{z}{1+z} \frac{z}{1+z} \frac{4z}{3+z} \frac{4z}{5+z} \frac{9z}{6+z} \cdots
$$

$$
\bullet e^z = 1 + \frac{z}{1-z} \frac{z}{1-z} \frac{z}{1-z} \frac{z}{1-z} \frac{z}{1-z} \frac{z}{1-z} \cdots
$$

$$
\bullet \tan z = \frac{z}{1-z} \frac{z^2}{1-z} \frac{z^2}{1-z} \cdots
$$

$$
\bullet \tanh z = \frac{z}{1+z} \frac{z^2}{1+z} \frac{z^2}{1+z} \cdots
$$

VI. FAST EVALUATION OF FUNCTIONS

Arithmetic operations such as +, - are fast to perform mentally and by machines, whereas ×, / are harder and slower to perform [19].

Thus, polynomial, rational, continued fraction evaluation of functions are slow methods for function evaluation since they need ×, / operations, extensively.

Several proposals were developed to bypass ×, / operations by making a correspondence between ×, / and +, - operations as for instance with $A^m A^n = A^{m+n}$ and $A^m / A^n = A^{m-n}$. This means that a function such as the logarithm will be able to achieve such correspondence.

A. John Napier and the invention of logarithms

Development of fast procedures for evaluation of functions started in the 17th century with Bürgi and Napier.

Joost Bürgi, a Swiss clockmaker [20] discovered first, logarithms as a computational tool. He did not publish his discovery, however, until 1620. Just a few years earlier, the method of logarithms was first published in 1614, in a book titled "Mirifici Logarithmorum Canonis Descriptio", by John Napier, a Scottish Baron. This method contributed to the advance of science, and especially astronomy, by making some difficult calculations possible. Prior to the advent of calculators and computers, it was used constantly in surveying, navigation, and other branches of practical mathematics. Besides their usefulness in computation, logarithms also filled an important place in higher theoretical mathematics.

Napier made the word logarithm [20] to specify a ratio from the Greek (logos) meaning proportion, and (arithmos) meaning number. Napier realized that the difference of two logarithms determines the ratio of their corresponding numbers, so that an arithmetic series of integers (logarithms) such as 1,2,3... corresponds to a geometric series of numbers $q^1, q^2, q^3...$ where $q$ is the ratio (cf. Fig. 6).

In order to achieve a specific range of values, it is required [20] to make the ratio $q$ of the geometric series close to 1. Napier chose $q = 1 - 10^{-7} = 0.999999$, and Bürgi chose $q = 1 + 10^{-4} = 1.0001$.

Napier evaluated $(1 - 10^{-7})^\ell$ for $\ell$ ranging from 1 to 100. The result for $\ell=100$ is approximately $0.99999 = 1 - 10^{-5}$ covering the interval $[1 - 10^{-7} : 10^7]$ since in his time, this interval was the most coveted for performing evaluation of mathematical functions [20]...

His efforts involving evaluation of products such as $10^7(1 - 10^{-7})^\ell$ for $\ell = 1..50$, spanned 20 years and allowed him to find, for any number $N \in [5 - 10^7]$, the exponent $\ell$ of the expression [20]:

$$
N = 10^7 (1 - 10^{-7})^\ell \quad (46)
$$
<figure><img src="image_9.png" /><figcaption>**Fig.6:** Correspondence between logarithmic scale generated by a geometric series $q^1, q^2, q^3...$ ($q < 1$) and a uniform scale mapping powers $q^1, q^2, q^3...$ to integers 1,2,3,...</figcaption></figure>

Napier called $\ell$ "logarithm" [20].

In present day mathematics, the relation can be transformed as:

$$ \ell = \frac{\ln\left(\frac{N}{10^7}\right)}{\ln\left(1 - 10^{-7}\right)} \qquad (47) $$

Using Euler result: $\lim_{n \to \infty} (1 - 1/n)^n = 1/e$ with $n = 10^7$ leads to $(1 - 10^{-7})^{10^7} \approx \frac{1}{e}$.

This can be rewritten as: $10^7 \ln(1 - 10^{-7}) \approx \ln(\frac{1}{e})$ or $\frac{1}{\ln(1 - 10^{-7})} = 10^7 \frac{1}{\ln(\frac{1}{e})}$.

Thus:

$$ \ell \approx 10^7 \frac{\ln\left(\frac{N}{10^7}\right)}{\ln\left(\frac{1}{e}\right)} = 10^7 \log_{1/e}\left(\frac{N}{10^7}\right) = -10^7 \log_e\left(\frac{N}{10^7}\right) = -10^7 \ln\left(\frac{N}{10^7}\right) \qquad (48) $$

This result means that Napier was using $1/e$ as a logarithm basis instead of $e$ proper to the function $\ln$.

Using logarithms meant consulting tables of logarithms, which had to be created manually.

In 1617, Henry Briggs published the first version of his own table of common logarithms, containing the logarithms of all integers below 1000 to eight decimal places. In 1624, he published Arithmetica Logarithmica table, containing the logarithms of all integers from 1 to 20,000 and from 90,000 to 100,000 to fourteen places of decimals, with a theory about the nature and use of logarithms.

Tables of logarithms were replaced, much later, by slide rules that engraved the logarithmic scale (cf. Fig. 6) on rulers, allowing to make products and divisions by simply moving the slider forward and backward.

Calculators and computers were developed afterwards to perform a lot more sophisticated calculations.

## B. CORDIC algorithms

CORDIC [21, 22] is based on the fact simple Boolean shift-add operations can perform several computing tasks such as evaluation of trigonometric, hyperbolic, logarithmic functions, real and complex multiplications, division, square-root, solution of linear systems, eigenvalue estimation, singular value decomposition (SVD) [21], QR matrix factorization etc... Consequently CORDIC has been applied to diverse areas such as signal and image processing, communication systems, robotics and 3D graphics in addition to general scientific and technical computation.

The CORDIC algorithm [21, 22] was introduced in 1959 by Jack E. Volder who developed at the aero-electronics department of Convair to replace the analog resolver in the B-58 bomber navigation computer. The main challenge was the real-time determination of current position over a spherical earth. It has been used in applications like 8087 co-processor, pocket calculators, radar signal processors, robotics etc. CORDIC is nothing but a computer which contains a special serial arithmetic unit consisting of shift and addition units, and special interconnections. By using a prescribed sequence of conditional add/subtract operations, the CORDIC arithmetic unit can be controlled and thus used to solve any function [21].

The Hewlett-Packard 9100 desktop calculator, built in 1968, used CORDIC methods for the trigonometric functions. In 1971, Walther [21] generalized this algorithm to compute logarithms, exponentials, and square roots. CORDIC is not the fastest way to perform multiplications or to compute logarithms and exponentials but, since the same algorithm allows the computation of most mathematical functions using very simple basic operations, it is attractive for hardware implementations.
CORDIC is so successful that it is used in calculators, CPU, GPU (Graphic Processing Unit) chips as well in DSP [23] (Digital Signal Processors) and FPGA (Field Programmable Gate Array) circuits. Thus CORDIC processors are used for signal processing applications (FFT, DFT), discrete Hartley transform, filtering...

The input and output vectors are related by Givens Rotation and whenever vector is normally rotated, the coordinates it gives are the accurate ones. But when a real rotation is replaced by a pseudo-rotation it gives erroneous coordinate values. This is because, whenever a vector is rotated, the length of final vector is larger than the actual required one.

The following represents a Givens rotation in x - y plane by angle α:

$$
\mathbf{R}(\alpha) = \begin{bmatrix} \cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{bmatrix} \qquad (49)
$$

For an arbitrary angle θ its sine and cosine are evaluated by dividing θ into smaller angles αᵢ such that θ = Σnᵢ=₀ αᵢ and Givens rotation is performed as follows:

$$
\mathbf{R}(\theta) = \mathbf{R} \left( \sum_{i=0}^{n} \alpha_i \right) = \prod_{i=0}^{n} \mathbf{R}(\alpha_i) = \prod_{i=0}^{n} M_i \qquad (50)
$$

with:

$$
M_i = \begin{bmatrix} \cos \alpha_i & -\sin \alpha_i \\ \sin \alpha_i & \cos \alpha_i \end{bmatrix} = \cos(\alpha_i) \begin{bmatrix} 1 & -\tan \alpha_i \\ \tan \alpha_i & 1 \end{bmatrix} \quad (51)
$$

Considering angles αᵢ such that tan αᵢ = 2⁻¹, matrix Mᵢ becomes:

$$
M_i = \begin{bmatrix} 1 & -1/2^i \\ 1/2^i & 1 \end{bmatrix} \qquad (52)
$$

Thus:

$$
\mathbf{R}(\theta) = \prod_{i=0}^{n} \cos(\alpha_i) \prod_{i=0}^{n} M_i \qquad (53)
$$

implying all matrix products [23] are done simply and quickly with right shift and add operations.

The scale factor Πnᵢ=₀ cos(αᵢ) is transformed into: Πnᵢ=₀ 1/√(1 + tan²αᵢ) leading to:

$$
K_n = \prod_{i=0}^{n} \cos(\alpha_i) = \prod_{i=0}^{n} \frac{1}{\sqrt{1 + 2^{-2i}}} \qquad (54)
$$

Finally:

$$
\mathbf{R}(\theta) = K_n \prod_{i=0}^{n} M_i \qquad (55)
$$

Scale factor K_n converges after 5 iterations toward: 0.6073518 and for large values of i toward 0.60725296.

These operations allow to calculate efficiently cos(θ), sin(θ), tan(θ), ...

After describing angles and rotations, there is need to represent ordinary numbers and elementary functions such as logarithm, square-root...

Similarly to the trigonometric approach, it is possible to perform rotations with hyperbolic functions taking angles αᵢ are such that tanh αᵢ = 2⁻¹ (the counterpart of the trigonometric relation tan αᵢ = 2⁻¹).

Rotations are accomplished by:

$$
\begin{pmatrix} c_{i+1} \\ s_{i+1} \end{pmatrix} = \begin{bmatrix} \cosh \alpha_i & \sinh \alpha_i \\ \sinh \alpha_i & \cosh \alpha_i \end{bmatrix} \begin{pmatrix} c_i \\ s_i \end{pmatrix}. \qquad (56)
$$
where a sign change in Eq. (56) leads to hyperbolic rotations.
Factoring cosh $\alpha_{n+1}$ as was done in the trigonometric case and introducing a hyperbolic scale factor

$$ \bar{K}_n = \prod_{i=1}^{n} \cosh(\alpha_i) = \prod_{i=1}^{n} \frac{1}{\sqrt{1 - 2^{-2i}}} \qquad (57) $$

all matrix products are similar to the trigonometric case and are done with shift and add operations. Scale factor $\bar{K}_n$ converges after 5 iterations toward: 1.2051361 and for large values of $i$ toward 1.2051363.

These operations allow to calculate efficiently cosh($\theta$), sinh($\theta$), tanh($\theta$), exp($\theta$), log($x$), $\sqrt{x}$...

### C. GPU architecture, AI, SFU and TPU

Graphics Processing Units (GPUs) are multiple core devices intended for performing computational graphics. They became increasingly used for high performance scientific computing, Neural Networks, Deep Learning and Artificial Intelligence (AI) [24] since their massively parallel architectures perform at least one order of magnitude faster than CPUs, and since general purpose languages adapted to GPU architecture like CUDA [25] that is similar to C and OpenCL [26] have been developed.

Google® has developed TPU (Tensor Processing Units) for dealing with deep learning problems of AI as well as cognitive problems such as Chess and Go games, Face recognition, Speaker recognition...

TPU-based processors are fast because of the massive parallelism of their architecture specially tailored for tackling stacked layers in deep neural network problems. They are faster than traditional CPU (Von Neumann architecture or with some modifications for speedup) and GPU (Graphic Processing Units) dedicated to graphic (Video) operations possessing a Harvard architecture like DSP (Digital signal processor) chips targeted for real-time applications.

Google TPU are made of parallel (systolic) arrays that contain each 65,536 (256 x 256) circuits called "Tensor Processing Elements (TPE)" that perform matrix multiplication with 8-bit based multiply-and-accumulate (MAC) operations in a single clock cycle. The TPU runs at 700 MHz, thus it can compute $65,536 \times 7 \times 10^8 = 46 \times 10^{12}$ MAC per second [27].

Performance comparisons between CPU, GPU and TPU are given in Table 3.

<table><caption>**Table 3:** Comparison of the CPU, GPU and TPU performances.</caption><thead><tr><td>Processor type</td><td>Number of Instructions per cycle</td></tr></thead><tbody><tr><td>CPU</td><td>~ 10<sup>0</sup></td></tr><tr><td>CPU (with vector extensions)</td><td>~ 10<sup>1</sup></td></tr><tr><td>GPU</td><td>~ 10<sup>4</sup></td></tr><tr><td>TPU</td><td>~ 10<sup>5</sup></td></tr></tbody></table>

From a data computing element point of view, a CPU processes (1 × 1) data units (such as product of two scalars), a GPU a (1 × N) data unit (such as the scalar product of two N dimensional vectors) whereas TPU processes (N × N) data units (such as the product of two N × N matrices). A CPU may perform tens of instructions per cycle (IPS), a GPU ~ 10⁴ IPS whereas a TPU may crunch several ~ 10⁵ IPS.

The success of GPU is built on architecture of the NVIDIA GPU initially used to test the Fermi processor architecture, GPU programming in CUDA and the divergence problems arising from the partial SIMD (Single Instruction Multiple Data) [19] execution on GPU are briefly described below.

From a hardware point of view, a GPU is composed of several Streaming Multiprocessors (SM) (Fermi C2070 contains 14), each being a SIMD unit. A SM is composed of multiple execution units (CUDA cores) (32 on Fermi) sharing the same pipeline and many registers (32768 on Fermi). GPU memory is organized in two levels: device memory, which can be accessed by any SM on the device; and shared memory, which is local to each SM [25].

From a software point of view, the developer writes in CUDA a scalar code for one function designed to be executed on the device, namely a kernel. At runtime, many threads or execution flows are created by blocks and bundled into
a grid to run the same kernel concurrently on the device. Each block is assigned to a SM. Within each block, threads are executed by groups of 32 called warps. The ratio of the number of resident warps (number of warps a SM can process at the same time) to the maximum number of resident warps per SM is named the occupancy. In order to increase the occupancy the number of blocks and their sizes have to be tuned.

GPUs since Fermi began to integrate an array of special hardware accelerators in the SMs, called Special-Functional Units (SFUs). The elementary transcendental functions include sine, cosine, division, exponential, power, logarithm, reciprocal, $\sqrt{x}$ and $1/\sqrt{x}$.

These accelerators perform complex transcendental operations and special function evaluations at the hardware level [25]. SFUs perform each operation with approximate algorithms while using IEEE 754 floating-point representation [19] in single precision. The SFU hardware module is composed of several components, the first one being a CORDIC module. This module performs estimations to obtain a valid result. All other ancillary components operate in parallel such that a single clock cycle is needed to obtain a final IEEE 754 [19] result.

## VII. APPROXIMATION OF FUNCTIONS BY REGRESSION

Approximating a function depending on a set of parameters $a$ such that $y = f(x; a)$ in order to describe experimental measurements ($x_i, y_i$) is generally a challenging problem in Science and Engineering [28]. Historically, Legendre and Gauss [28] developed in XIX century beginnings, a least-square method aiming at minimizing the quantity:

$$ \sum_{i=1}^{N} [y_i - f(x_i; a)]^2 $$

This leads to an approximate line $y = a + bx$, with coefficients given by:

$$ a = \langle y \rangle - b \langle x \rangle, \quad b = \frac{\sum (x_i - \langle x \rangle) y_i}{\sum (x_i - \langle x \rangle)^2}, $$

where $\langle x \rangle$ and $\langle y \rangle$ are mean values of $x_i$ and $y_i$.

Standard deviations of the coefficients [28] are given by:

$$ \sigma_a = \sqrt{\frac{\sum (y_i - bx_i - a)^2}{N - 2} \left[ \frac{1}{N} + \frac{\langle x \rangle^2}{\sum (x_i - \langle x \rangle)^2} \right]}, \quad \sigma_b = \sqrt{\frac{\sum (y_i - bx_i - a)^2}{(N - 2) \sum (x_i - \langle x \rangle)^2}} \quad (58) $$

Usually with same data $k$ is evaluated better than $b$ ($\sigma_k < \sigma_b$).

Consider a statistical model with two variables $X$ and $Y$, where one tries to predict $Y$ from $X$. The assumptions of the model are as follows:

1. The distribution of $X$ is arbitrary (and perhaps $X$ is even non-random).

2. If $X = x$, then $Y = \beta_0 + \beta_1 x + \epsilon$, for some parameters $\beta_0$ and $\beta_1$, and some random noise variable $\epsilon$.

3. $E[\epsilon|X = x] = 0$ (no matter what $x$ is), $Var[\epsilon|X = x] = \sigma^2$ (no matter what $x$ is).

4. $\epsilon$ is uncorrelated across observations.

When dealing with multiple data points, $(X_1, Y_1), (X_2, Y_2), \dots (X_n, Y_n)$, consider for each $i \in 1:n$,

$$ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i \quad (59) $$

where the noise variables $\epsilon_i$ all have the same expectation (0) and the same variance ($\sigma^2$), and $Cov[\epsilon_i, \epsilon_j] = 0$ (unless $i = j$).
A. Least Squares Estimates

An alternative way of estimating the simple linear regression model starts from the objective being sought, rather than from the formula for the slope. The true optimal slope and intercept [28] are the ones which minimize the Mean Squared Error $\langle \epsilon^2 \rangle$:

$$ (\beta_0, \beta_1) = \underset{(b_0, b_1)}{\operatorname{argmin}} \mathbf{E} \left[ (Y - (b_0 + b_1 X))^2 \right] \quad (60) $$

This is a function of the complete distribution, that can be approximated it from data. Data MSE is expressed as:

$$ \langle \epsilon^2 \rangle(b_0, b_1) \equiv \frac{1}{n} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_i))^2 \quad (61) $$

Notice that the MSE is a function of $b_0$ and $b_1$; it is also a function of $(x_1, y_1)$, $(x_2, y_2)$, ..., $(x_n, y_n)$ data. Minimizing $\langle \epsilon^2 \rangle(b_0, b_1)$ consists of taking the derivatives with respect to the slope and the intercept:

$$ \frac{\partial \langle \epsilon^2 \rangle}{\partial b_0} = \frac{1}{n} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_i))(-2), \quad \frac{\partial \langle \epsilon^2 \rangle}{\partial b_1} = \frac{1}{n} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_i))(-2x_i) \quad (62) $$

Set these to zero at the optimum $(\hat{\beta}_0, \hat{\beta}_1)$:

$$ \frac{1}{n} \sum_{i=1}^{n} (y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i)) = 0, \quad \frac{1}{n} \sum_{i=1}^{n} (y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i))(x_i) = 0 $$

These are often called the normal equations for least-squares estimation, or the estimating equations: a system of two equations in two unknowns, whose solution gives the estimate. Thus:

$$ \bar{y} - \hat{\beta}_0 - \hat{\beta}_1 \bar{x} = 0, \quad \bar{x}y - \hat{\beta}_0 \bar{x} - \hat{\beta}_1 \bar{x}^2 = 0 \quad (63) $$

The first equation, when rewritten, gives

$$ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} \quad (64) $$

Substituting this into the remaining equation,

$$ 0 = \bar{x}\bar{y} - \langle y \rangle \langle x \rangle + \hat{\beta}_1 \langle x \rangle \langle x \rangle - \hat{\beta}_1 \bar{x}^2, \quad \text{thus} \quad \hat{\beta}_1 = \frac{c_{XY}}{s_X^2} \quad (65) $$

The equivalence between the simple estimator and the least-squares estimator is somehow a property of linear models.

B. Bias, Variance and Standard Error of Parameter Estimates

The slope, $\hat{\beta}_1$ is given by:

$$ \hat{\beta}_1 = \frac{c_{XY}}{s_X^2} = \beta_1 + \frac{\frac{1}{n} \sum_{i=1}^{n} x_i \epsilon_i - \langle x \rangle \langle \epsilon \rangle}{s_X^2} \quad (66) $$

Since $\langle x \rangle \langle \epsilon \rangle = n^{-1} \sum_i \langle x \rangle \epsilon_i$,

$$ \hat{\beta}_1 = \beta_1 + \frac{\frac{1}{n} \sum_{i=1}^{n} (x_i - \langle x \rangle) \epsilon_i}{s_X^2} \quad (67) $$

This representation of the slope estimate shows that it is equal to the true slope ($\beta_1$) plus a contribution depending on noise characteristics ($\epsilon_i$, and their sample average $\langle \epsilon \rangle$). This is used to find the expected value [28] and the variance of the estimator $\hat{\beta}_1$.
$x_i$ are treated as non-random variables when experiments are controlled or designed specifically for that purpose. In randomized experiments or in observational studies, the $x_i$ are not necessarily fixed; however, these expressions will be correct for the conditional expectation $\mathbf{E}[\hat{\beta}_1|x_1, \dots, x_n]$ and conditional variance $\text{Var}[\hat{\beta}_1|x_1, \dots, x_n]$.

From $\mathbf{E}[\epsilon_i|X_i] = 0$, one gets:

$$ \frac{1}{n} \sum_{i=1}^{n} (x_i - \langle x \rangle) \mathbf{E}[\epsilon_i] = 0, \quad \text{thus} \quad \mathbf{E}[\hat{\beta}_1] = \beta_1 \qquad (68) $$

Since the bias of an estimator is the difference between its expected value and the truth, $\hat{\beta}_1$ is an unbiased estimator of the optimal slope.

The term bias means $\mathbf{E}[\hat{\beta}_1] - \beta_1$. An unbiased estimator could still make systematic mistakes, for instance, it could underestimate 99% of the time, provided that the 1% of the time it over-estimates, it does so by much more than it under-estimates. Moreover, unbiased estimators are not necessarily superior to biased ones: the total error depends on both the bias of the estimator and its variance.

Turning to the intercept,

$$ \mathbf{E}[\hat{\beta}_0] = \mathbf{E}[\bar{Y} - \hat{\beta}_1 \bar{X}] = \beta_0 \qquad (69) $$

thus it is also unbiased.

Using the formula for the variance of a sum, and the model assumption that all the $\epsilon_i$ are uncorrelated with each other,

$$ \text{Var}[\hat{\beta}_1] = \text{Var}[\beta_1] + \frac{\frac{1}{n} \sum_{i=1}^{n} (x_i - \langle x \rangle) \epsilon_i}{s_X^2} = \frac{\sigma^2}{ns_X^2} \qquad (70) $$

This means the variance of the slope estimate goes up as the noise around the regression line ($\sigma^2$) gets bigger, and goes down as there are more observations ($n$), which are further spread out along the horizontal axis ($s_X^2$).

The standard error of an estimator is just its standard deviation, or the square root of its variance:

$$ \text{se}(\hat{\beta}_1) = \frac{\sigma}{\sqrt{ns_X^2}} \qquad (71) $$

### C. Maximum Likelihood vs. Least Squares

The likelihood [28] of a parameter value on a data set is the probability density at the data under those parameters. Using the Gaussian-noise model, it is possible to write down a likelihood. By model assumptions, the parameters are $b_0, b_1, s^2$, then $Y|X = x \sim N(b_0 + b_1x, s^2)$, and $Y_i$ and $Y_j$ are independent given $X_i$ and $X_j$, so the over-all likelihood is

$$ \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi s^2}} e^{-\frac{(y_i - (b_0 + b_1 x_i))^2}{2s^2}} \qquad (72) $$

Using log-likelihood amounts to take logarithm of above expression to obtain:

$$ L(b_0, b_1, s^2) = -\frac{n}{2} \log 2\pi - n \log s - \frac{1}{2s^2} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_i))^2 \qquad (73) $$

In Statistics, a likelihood function is generally maximized. This means one should find a parameter value which makes the data observed as likely, as probable, as the model permits. Thus derivatives are set to zero:

$$ \frac{\partial L}{\partial b_0} = -\frac{1}{2s^2} \sum_{i=1}^{n} 2(y_i - (b_0 + b_1 x_i))(-1), \quad \frac{\partial L}{\partial b_1} = -\frac{1}{2s^2} \sum_{i=1}^{n} 2(y_i - (b_0 + b_1 x_i))(-x_i) \qquad (74) $$

Notice that when these derivatives are set to zero, all the multiplicative constants, in particular, the prefactor of $\frac{1}{2s^2}$ go away. Thus

$$ \sum_{i=1}^{n} y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i) = 0, \quad \sum_{i=1}^{n} (y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i))x_i = 0 \qquad (75) $$
These are, up to a factor of $1/n$, exactly the equations obtained from the method of least squares (Eq. 63). That means that the least squares solution is the maximum likelihood estimate under the Gaussian noise model.

The derivative with respect to $s$ is:

$$ \frac{\partial L}{\partial s} = -\frac{n}{s} + 2\frac{1}{2s^3}\sum_{i=1}^{n}(y_i - (b_0 + b_1x_i))^2 \quad (76) $$

Setting this to 0 at the optimum, including multiplying through by $\hat{\sigma}^3$, yields:

$$ \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}(y_i - (\hat{\beta}_0 + \hat{\beta}_1x_i))^2 \quad (77) $$

Notice that the right-hand side is just the data sample mean squared error.

Maximum likelihood estimates (MLE) of the regression curve coincide with least-squares estimates when the noise around the curve is additive, Gaussian, of constant variance, and both independent of $X$ and of other noise terms. If any of those assumptions fail, maximum likelihood and least squares estimates can diverge, though sometimes the MLE solves a generalized least squares problem.

## D. Non-linear Least Squares: Levenberg-Marquardt method

Levenberg-Marquardt algorithm solves non-linear least-squares problem by combining gradient descent method and Gauss-Newton method to minimize a $\chi^2(\mathbf{a})$ statistical measure [7, 29] (called also merit function) defined by:

$$ \chi^2(\mathbf{a}) = \sum_{i=1}^{N} \left[ \frac{y_i - f(x_i; \mathbf{a})}{\sigma_i} \right]^2 \quad (78) $$

where $x_i, y_i, i = 1...N$ are sampled data values and $f(x; \mathbf{a})$ is an approximate function representing the data through an ensemble of $m$ parameters represented by $\mathbf{a}$.

$\sigma_i$ represents the measurement error of the $i$-th data point.

The Levenberg-Marquardt algorithm (LMA) addresses the problem pertaining to minimize a function $\chi^2(\mathbf{a})$ according to a set of $m$ unconstrained parameters $\mathbf{a}$, as long as the second derivatives of $\chi^2(\mathbf{a})$ exist.

Given a multidimensional function $\chi^2(\mathbf{a})$ it amounts to finding zeroes $\mathbf{a}^*$ of the partial derivatives $\frac{\partial\chi^2(\mathbf{a})}{\partial a_i}$, $i = 1...m$. The latter $\mathbf{a}^*$ is maximum or minimum depending on whether the Hessian matrix $H(\chi^2)$ (or local convexity) whose components are defined by the second derivatives: $H_{ij} = (\frac{\partial^2\chi^2(\mathbf{a})}{\partial a_i \partial a_j})$, $i, j = 1...m$, is positive or negative at $\mathbf{a}^*$.

Gauss-Newton method to minimize $\chi^2(\mathbf{a})$ by minimizing the residuals $r_i(\mathbf{a})$, $i = 1...m$ expressed with respect to $\mathbf{a}$, $i = 1...m$ components.

In order to minimize $\chi^2(\mathbf{a})$, a Newton step ($Nt$) is taken, which does the following update:

$$ v_{Nt} = -\nabla^2[\chi^2(\mathbf{a})]^{-1}\nabla\chi^2(\mathbf{a}) = -(H^T H + \sum_{i=1}^{m} r_i(\mathbf{a})\nabla^2 r_i(\mathbf{a}))^{-1}H^T\chi^2(\mathbf{a}) \quad (79) $$

where $H$ is the Hessian matrix of $\chi^2(\mathbf{a})$. In the Gauss-Newton (GN) step, the second term is discarded in the inverse term and yields:

$$ v_{GN} = -(H^T H)^{-1}H^T\chi^2(\mathbf{a}) \quad (80) $$

The Gauss-Newton step always provides a descent direction when $H$ has full column rank; because $2H^T H$ is positive definite [28].

The Levenberg-Marquardt algorithm is a variant of the Gauss-Newton method add a regularization term $\lambda$ to avoid Hessian potential singularities.

$$ v_{GN} = -(H^T H + \lambda I)^{-1}H^T\chi^2(\mathbf{a}) \quad (81) $$

Here $\lambda$ defines a trust region where the quadratic approximation is a good fit. It blends both linear and quadratic approximations. When $\lambda = 0$, the update rule is completely based on the second-order approximation, when $\lambda$ is very large, the rule is completely based on the first-order approximation.
When the problem encounters a unique solution, the minimum is reached whatever the chosen initial values.

As in any iterative algorithm, convergence of LMA is achieved when convergence criteria are fulfilled. Convergence may be defined according to following criteria:

* parameters stability: $\sum_{j=1}^{m} (\mathbf{a}_{j}^{(k+1)} - \mathbf{a}_{j}^{(k)})^{2} < \epsilon_{a}$

* objective function stability: $|[\chi^{2}(\mathbf{a})]^{(k+1)} - [\chi^{2}(\mathbf{a})]^{(k)}| < \epsilon_{b}$

Actually these criteria [29, 30], do not ensure a convergence towards a saddle point instead of a true minimum. Therefore the parameter and objective function stability should account for the relative distance to minimum. Generally the Hessian matrix is required to be invertible avoiding convergence toward a saddle point, and ensuring a minimum is finally reached.

Thresholds $\epsilon_{a}$ and $\epsilon_{b}$ can be adjusted, with values around 0.001 sufficient usually to warrant good convergence.

Fig. 7 shows an application of LMA to Water ionization spectroscopy [31] requiring a multi-Gaussian Peakfit [7] such that the spectrum is written as: $\sum_{k=1}^{K} A_{k} \exp \left[-\left(\frac{x-\mu_{k}}{\sigma_{k}}\right)^{2}\right]$ where $x$ is electron energy (in eV).

<figure><img src="image_10.png" /><figcaption>Fig.7: Experimental collisional ionization spectrum [31] of gaseous Water versus ejected electron energy (in eV) and a multi-Gaussian decomposition Peakfit resulting in two Gaussians with the following parameters: Gaussian # 1 $A_{1}$= 63.52, $\mu_{1}$= 1.92 eV, $\sigma_{1}$= 1.08 eV and Gaussian # 2 $A_{2}$= 66.49, $\mu_{2}$= 3.83 eV, $\sigma_{2}$= 0.40 eV</figcaption></figure>

Note that LMA is widely used in many applications, including Neural Network training, Deep Learning and Artificial Intelligence.

As an example is Fig.8 where a Neural Network rebuilds the original sunspot sequence [32] and the fitting [33] is done by LMA.

The Neural Network procedure combined with LMA is done by considering:

$$y(n + 1) = F[y(n), y(n - 1), y(n - 2), y(n - 3), \dots, y(n - m + 1)] \quad (82)$$

where $F$ an $m$ dimensional map and $y(n)$ the value of the signal $x(t)$ (sunspot series) sampled at timestep $n$ and expressing $F$ as a sum of linear and sigmoidal (such as $g(x) = \frac{1}{2}(1 + \tanh(x))$) functions [34]. $m$ is the number of hidden layers in the Neural Network.

The linear function is simply the linear interpolation of the sequence. The unknown weights multiplying $g(x)$ are determined [34] via a non-linear least-square algorithm in order to achieve the best fit to the measured data.

Lapedes and Farber [35] built a $m = 5$ hidden layer Neural Network to reach $10^{-4}$ accuracy between the fully chaotic 1D logistic map $x(i + 1) = 4x(i)(1 - x(i))$ and the one-step prediction given by:
$$x(i+1) = (0.31x(i) - 2.04) - 0.64g(-1.11x(i) - 0.26) - 1.3g(2.22x(i) - 1.71) \\ - 2.285g(3.91x(i) + 4.82) - 3.905g(2.46x(i) - 3.05) + 5.99g(1.68x(i) + 0.60) \quad (83)$$

The first term $(0.31x(i) - 2.04)$ is a linear interpolation of the logistic sequence whereas all numerical values related to $g(x)$ are evaluated with the 5 hidden layer Neural Network.

<figure><img src="image_11.png" /><figcaption>**Fig.8:** Comparison between Non-linear one step predicted and original Sunspot sequences. Prediction is made with Neural Network rebuilt sequence [33] via non-linear least squares LMA.</figcaption></figure>

In order to make one-step non-linear prediction of the sunspot data [32], the linear function is simply a linear interpolation of the sunspot sequence as in the logistic case. The unknown weights related to $g(x)$ are determined [34] through the LMA in order to achieve the best fit to the sunspot data.

The results comparing a one-step prediction to sunspot data are displayed in Fig.8. Goodness of fit between the observed and predicted envelopes for a map of dimension $m = 5$ is confirmed by performing a study of the prediction error versus dimension $m$. A minimum is indeed observed in the prediction error for dimension $m$ equal to 5 or 6.
[1] H. C. Von Baeyer, *The Fermi Solution: essays on science*, Dover (2001).

[2] L. Hodgkin, *A History of Mathematics From Mesopotamia to Modernity*, Oxford University Press Inc., New York (2005)

[3] E. Maor, *e: The Story of a Number*, Princeton University Press (1994).

[4] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[5] K. Busawon, *Upper and lower bounds for $\pi$ based on Vieta geometrical approach*, Notes on number theory and discrete mathematics, 22, 62 (2016).

[6] T. J. Osler, *Vieta-like products of nested radicals with Fibonacci and Lucas numbers*, The Fibonacci Quarterly, 52, 202 (2007). See also *The united Vieta and Wallis product for $\pi$*, American Math. Monthly 106, 774 (1999).

[7] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, Second Edition, Cambridge University Press (New-York, 1992).

[8] E. Salamin, *Computation of Pi Using Arithmetic-Geometric Mean*, Mathematics of Computation, vol. 30, 565 (1976)

[9] R. P. Brent, *Fast Multiple-Precision Evaluation of Elementary Functions*, Journal of the ACM, vol. 23, 242 (1976)

[10] M. Ablowitz and A. Fokas, *Complex variables: Introduction and Applications*, 2nd edition, Cambridge University Press, New-York (2003).

[11] C. M. Bender, S. A. Orszag, *Advanced Mathematical Methods for Scientists and Engineers*, McGraw-Hill, New York, (1978).

[12] A. B. Migdal, *Qualitative methods in Quantum theory*, W. A. Benjamin, Inc. Advanced Book Program, Reading Mass. (1977).

[13] M. Abramowitz and I.S Stegun, *Handbook of Mathematical Tables*, Dover, New-York (1960). M. Abramowitz, *On the practical evaluation of integrals*, Journal of the Society for Industrial and Applied Mathematics (SIAM), Volume 2, (1954)

[14] C. Brezinski, *History of continued fractions and Padé approximants* Springer-Verlag, Berlin Heidelberg (1991).

[15] J. Stoer R. and Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[16] G. A. Baker Jr., *Essentials pf Padé Approximants*, Academic Press, New-York (1975).

[17] A. Ya. Khinchin, *Continued fractions*, The University of Chicago Press, London (1964).

[18] W. M. Jones, W. J. Thron, *Continued fractions. Analytic theory and applications* Addison-Wesley Publ. Co., MA (1980).

[19] J. L. Hennessy and D. A. Patterson, *Computer Architecture: A Quantitative Approach*, Fourth Edition, Morgan Kaufmann Publishers an Imprint of Elsevier, San Francisco (2007)

[20] E. W. Hobson, *John Napier and the invention of logarithms*, 1614, Cambridge: The University Press (1914)

[21] J. S. Walther, *A Unified Algorithm for Elementary Functions*. In: Proceedings of the May 18-20, 1971, Spring Joint Computer Conference, AFIPS '71, 379 (1971)

[22] P. K. Meher, J. Valls, T.B. Juang, K. Sridharan and K. Maharatna, *50 Years of CORDIC: Algorithms, Architectures, and Applications*, IEEE Trans. Circ. and Sys **IEEE CAS 56**, 1893 (2009)

[23] J. Y. Stein, *Digital Signal Processing: A computer science perspective*, John Wiley and Sons, New-York (2000).

[24] J. J. Hopfield. *Neural networks and physical systems with emergent collective computational abilities*, Proceedings of the National Academy of Sciences of the United States of America, 79, 2554 (1982).

[25] NVIDIA, *CUDA C Best Practices Guide, version 4.1*, (2012).

[26] Khronos Group, *The OpenCL Specification Version 1.2*, (2011).

[27] N. P. Jouppi et al., *In-Data-Center Performance Analysis of a Tensor Processing Unit*, International Symposium on Computer Architecture, 45(2), 1-12 (2017).

[28] P. Bevington and D K Robinson, *Data Reduction and error analysis for the Physical Sciences*, third edition, McGraw-Hill Higher Education, New-York (2003).

[29] K. Levenberg, *A Method for the Solution of Certain nonlinear Problems in Least Squares*, Quarterly of Applied Mathematics. 2 (2): 164 (1944). D. Marquardt, *An Algorithm for Least-Squares Estimation of Nonlinear Parameters*, SIAM Journal on Applied Mathematics. 11 (2): 431 (1963)

[30] R. Fletcher, *A Modified Marquardt Subroutine for Nonlinear Least Squares*, Technical Report AERE-R 6799, Theoretical Physics Division, Atomic Energy Research Establishment, Harwell, Berkshire (1971)

[31] S. Falcinelli, F. Pirani and F. Vecchiocattivi, Atmosphere, 6, 299 (2015)

[32] S. L. Marple Jr, *Digital Spectral Analysis with Applications*, Prentice-Hall, New-York (1987).

[33] M K Transtrum, B B Machta, J P Sethna, *Geometry of nonlinear least squares with applications to sloppy models and optimization*, Phys. Rev. E, PRE-Statistical, Nonlinear, and Soft Matter Physics 83, 271 (2011).

[34] H. White, *Some asymptotic results for learning in single hidden-layer feedforward network models*, J. Am. Stat Association **84**, 1003 (1989).

[35] A. Lapedes and B. Farber, *Non-linear signal processing using Neural Networks: Prediction and system modelling*, Los Alamos Report: LA-UR 87-2662 (1987).