Approximation to Elementary Transcendental Functions

1 Introduction

The second programming language I ever learned was the C language. About the first half of the ’80s, I got The C Programming Language book (Kernighan and Ritchie 1978). Well written, easy to understand, I love it. At that time, a friend of mine in college gave me a copy of “Lattice C,” a C compiler, to run it in the DOS command line, of course. At that time internet didn’t exist; to browse some books, you had to go to a library. With the C compiler, I practiced the exercises in the book.

However, something was missing. The compiler came without a mathematical library to compute sine, cosine, exponential function, etc. Those are called elementary transcendental functions. I was studying Mechanical Engineering at that time, and those functions were vital to solving college problems. Those days I had a part-time job teaching the Basic language to the owners of the ZX81 (Times Sinclair) home computer. The folks that hired me for the classes gave me the ZX81 manual. I asked them if they had something about the computing of those transcendental functions. Fortunately, they have it. They lent me a book describing the ZX81 ROM where there was an explanation about the source of the implementation of their transcendental functions.

I found that those transcendental functions, calculated by the ZX81 ROM, were using formulas from the book “Handbook of Mathematical Functions by Abramowitz and Stegun” (Abramowitz and Stegun 1964). Which eventually, I found it in the School of Mathematics Library at my college. Finally! I had the formulas referenced by the ROM book!

Figure 1.1: All the approximations in the Handbook of Mathematical Functions has 10 or less digits of precision.

Unfortunately, those formulas in the book have just ten digits of precision, and in C language, I had double-precision arithmetic with 16 precisión digits. I had to find a way to compute the coefficients in those formulas, but this time with 16 precision digits. Nonetheless, that was not an easy task. The maths behind the coefficients’ computation were complicated and valuable enough that the handbook gives the authorship attribution to the Los Alamos Scientific Laboratory and has to ask permission to publish the Chebyshev approximation from Clenshaw (Clenshaw 1954). Another hint that they were not easy is that they had numerical form only, and they were from the 1950’.

Long story short, I had to evaluate growingly complicated integrals over the so-called Chebyshev Polynomials (Gil, Segura, and Temme 2007) to get the coefficients. In the end, for some transcendental functions, I could find analytical solutions for these integrals and get the desired approximation.

The beauty of the analytical form is that you can have arbitrary precision. You are limited only by the accuracy of the arithmetic you use. In the code below, we shall test the \(atan\) approximation using 30 precision digits.

For \(atan\), the inverse of the tangent function, I found the following approximation:

\[ \begin{align} atan(x) &= \sum_{k=1}^\infty b_k·P_k(x) ~, ~~~\text{for}~ \lvert x \rvert \leq 1 \\ b_k &= (-1)^{(k-1)}\, (\frac{2}{2\,k-1}) \, (\sqrt{2}-1)^{(2k-1)} \\[5pt] P_{0}(x) &= x \\ P_{1}(x) &= x \\ P_{k}(x) &= (4\,x^2-2)·P_{k-1} - P_{k-2} \end{align} \]

We only need the squared root of \(2\) and arithmetic to compute the \(atan\) function with these formulas. The Handbook has the values of some irrational numbers with a lot of digits of precision, including \(\sqrt{2}\).

from mpmath import mp, mpf
mp.dps = 30 # Number of precision digits

def atan_poly_cheb(x):
    atan_poly_cheb.p2 = 4*x*x-2
    atan_poly_cheb.p_nm1 = x
    atan_poly_cheb.p_nm2 = x
    yield x
    while True:
        p_np1 = atan_poly_cheb.p2*atan_poly_cheb.p_nm1 \
          - atan_poly_cheb.p_nm2
        atan_poly_cheb.p_nm2 = atan_poly_cheb.p_nm1
        atan_poly_cheb.p_nm1 = p_np1
        yield p_np1

We will use Python code to test this approximation works as intended. The polynomials \(P_k(x)\) are computed with the atan_poly_cheb(x) function, and the \(b_k\) coefficients with atan_coef(). Both are Python generators, so you can call them as many times as you want and get as many values you need. We use the mpmath module to compute \(atan\) with 30 precision digits.

def atan_coef():
    SQRT2_M1 = mp.sqrt(2) - mpf(1)
    SQRT2_M1_2 = SQRT2_M1*SQRT2_M1
    atan_coef.k = mpf(0)
    atan_coef.rho_factor = mpf(-2)/SQRT2_M1

    while True:
        atan_coef.k += mpf(1)
        atan_coef.rho_factor *= - SQRT2_M1_2
        coef = atan_coef.rho_factor/(2*atan_coef.k-1)
        yield coef

The function mpf(x), in the mpmath module, creates a floating point number with multi-precision, 30 digits in our case. Then the arithmetic with these numbers maintains that precision. The computation of \(b_k\) is made in a iterative way:

\[ \begin{align} rho\_factor_0 &= -2/(\sqrt{2}-1) \\ rho\_factor_k &= - rho\_factor_{k-1}\times(\sqrt{2}-1)^2 \\[3pt] b_k &= \frac{rho\_factor_k}{2\,k-1} \end{align} \] The atan_approx(x) function calls the coefficient and Chebyshev polynomial generators, and with 37 values from each of them, compute the atan function with 30 digits of precision.

def atan_approx(x):
    result = 0
    cheby = atan_poly_cheb(x)
    coef = atan_coef()

    for i in range(37):
        result += next(coef) * next(cheby)
        
    return result

We evaluate the approximation for \(x=1\) where the \(atan(1)\) is \(\pi/4\). We multiply the result by \(4\) and compare it with the ground truth value of \(\pi\) provided by the mpmath module. Also we evaluate \(atan\) for a random number in \([-1,1]\) and compare the result, again with the value according to mpmath.

def main():
    pi_approx = atan_approx(1)*4
    print("*** Approximation to pi ****")
    print("Pi Approx: ", pi_approx)
    print("Pi Exact : ",mp.pi)
    print("Pi Approx Error: ", pi_approx-mp.pi)
    
    from random import random
    x = mpf(random()*2 - 1)
    approx = atan_approx(x)
    exact = mp.atan(x)
    
    print("\n*** atan of a random number ****")
    print("x: ", x)
    print("Approx atan: ", approx)
    print("Exact atan : ", mp.atan(x))
    print("Atan Approx Error: ", approx-exact)

*** Approximation to pi ****
Pi Approx:  3.14159265358979323846264338328
Pi Exact :  3.14159265358979323846264338328
Pi Approx Error:  3.54987407349455312435277854377e-30

*** atan of a random number ****
x:  -0.817895132505307209669354051584
Approx atan:  -0.685557759217410550765244311608
Exact atan :  -0.685557759217410550765244311608
Atan Approx Error:  0.0

At the end we achieve the 30 digits precision, showing the approximation works.

In this post we will introduce the Chebyshev Approximation, and show you how did I found the one for the \(atan\) function. I hope you like it.

2 The Chebyshev Approximation

You can approximate a function with the following formulas:

\[ \begin{align} f(x) &\approx \frac{a_0}{2} + a_1·T_1(x)+a_2·T_2(x)+ \dots a_p·T_p(x) \tag{2.1} \\[3pt] a_k(f) &= \frac{2}{\pi} \int_{-1}^1 \frac {f(x) · T_k(x) }{\sqrt{1-x^2}} dx \tag{2.2} \end{align} \]

Where \(T_k(x)\) are the Chebyshev polynomials of the first kind:

\[ \begin{align} T_0(x) &= 1 \\ T_1(x) &= x \\ T_2(x) &= 2x^2-1 \\ T_{n+1}(x) &= 2x·T_{n}(x)-T_{n-1}(x) \\ T_{n+2}(x) &= (4x^2-2)·T_{n}(x)-T_{n-2}(x) \end{align} \]

The trick with these polynomials is that if \(x=cos(\theta)\) then \(T_k(x)=cos(k\, \theta)\). In other words, you can express \(cos(k\,\theta)\) as a polynomial of \(cos(\theta)\). For example, the trignometric identity \(cos(2\,\theta)=2·cos^2(x)-1\) means that if we do the change of variable \(x=cos(\theta)\), then \(T_2(x)=cos(2\,\theta)=2x^2-1\).

If in (2.2) we make the change of variable \(x=cos(\theta)\), and substitute \(T_k(x)=cos(k\, \theta)\), we get:

\[ \begin{align} a_k(f) &= \frac{2}{\pi} \int_0^\pi {f(cos(\theta))·cos(k~\theta) } \,d\theta \tag{2.3} \end{align} \]

An Odd function

In math when \(f(-x)=f(x)\) you say that \(f\) is even, and if \(f(-x)=-f(x)\) then \(f\) is odd. The \(T_k\) polynomials for k even, contain only even powers of \(x\), and for \(k\) odd only odd powers of \(x\); which causes \(T_k(x)\) is an even function when \(k\) is even, and odd when \(k\) is odd.

If we split the integration interval in (2.2), and make a change of variable to get both integrals in \([0,1]\):

\[ \begin{aligned} a_k(f) &= \frac{2}{\pi} \int_{-1}^1 \frac {f(x) · T_k(x) }{\sqrt{1-x^2}} dx \\ &= \frac{2}{\pi} \int_{-1}^0 \frac {f(x) · T_k(x) }{\sqrt{1-x^2}} dx + \frac{2}{\pi} \int_0^1 \frac {f(x) · T_k(x) }{\sqrt{1-x^2}} dx \\ &= \frac{2}{\pi} \int_0^1 \frac {f(x) · T_k(x) + f(-x) · T_k(-x) }{\sqrt{1-x^2}} dx \end{aligned} \]

Now we can see for \(f\) even and \(k\) odd \(f(-x)·T_k(-x) = -f(x)·T_k(x)\) and \(a_k(f)\) vanishes; the same happens when \(f\) is odd and \(k\) is even. In summary, we can say:

\[ \begin{aligned} \text{When } f ~\text{is even:} & \\ a_k(f) &= \begin{cases} a_k(f)\neq 0 & \text{if } k ~\text{is even,} \\ 0 & \text{if } k ~\text{is odd} \end{cases} \\ \text{When } f ~\text{is odd:} & \\ a_k(f) &= \begin{cases} 0 & \text{if } k ~\text{is even,} \\ a_k(f)\neq 0 & \text{if } k ~\text{is odd} \end{cases} \end{aligned} \]

2.1 Approximation using the Derivative of the Function

Sometimes it is easier to find the approximation using the function first derivative instead of the function itself. To find the required formula, we will use the Integration by Parts rule, which you can remember from its wikipedia article.

\[ \begin{aligned} cos(k\,\theta) &= \frac{d(sin(k\,\theta)/k)}{d\theta}, \\ a_k(f) &= \frac{2}{\pi} \int_0^\pi {f(cos(\theta))·cos(k\,\theta) } \,d\theta \\ \implies \\ &= \frac{2}{\pi} f(cos(\theta))·\frac{sin(k\,\theta)}{k} \Big|_0^\pi + \\ &\frac{2}{k\,\pi} \int_{0}^\pi f'(cos(\theta))·sin(\theta)·sin(k\,\theta) \,d\theta \\[8pt] a_k(f) &= \frac{2}{k\,\pi} \int_0^\pi f'(cos(\theta))·sin(\theta)·sin(k\,\theta) \,d\theta \end{aligned} \]

Now we will go to the polynomial form by changing \(x=cos(\theta)\):

\[ \begin{align} a_k(f) &= \frac{2}{k\,\pi} \int_{-1}^1 \frac{f'(x)·sin(k\,\theta)·sin(\theta)}{\sqrt{1-x^2}}\,dx \\ &= \frac{2}{k\,\pi} \int_{-1}^1 \frac{sin(k\,\theta)}{sin(\theta)} · f'(x)·\sqrt{1-x^2} \,dx \tag{2.4} \end{align} \]

Where the quotient \(\frac{sin(k\,\theta)}{sin(\theta)}\) in (2.4) are the Chebyshev polynomials of the second kind \(U_k\).

\[ \begin{align} U_0(x) &= 1 \\ U_1(x) &= 2x \\ U_2(x) &= 4x^2 - 1 \\ U_{n+1}(x) &= 2x·U_n(x)-U_{n-1}(x) \\ U_{n+2}(x) &= (4x^2-2)·U_n(x)-U_{n-2}(x) \end{align} \]

The quotient \(\frac{sin(k\,\theta)}{sin(\theta)}\) in (2.4), for \(k=2\) is \(2·cos(\theta)=2\,x\), which es equal to \(U_1(x)\); that means we should replace that quotient with \(U_{k-1}\).

\[ a_k(f) = \frac{2}{k\,\pi} \int_{-1}^1 U_{k-1}(x) · f'(x)·\sqrt{1-x^2} \,dx \tag{2.5} \]

3 Approximation to atan

For the inverse of the tangent function \(atan\), we can get its series expansion by using sageMath (The Sage Developers 2021) online:

var('x')
f = atan(x)
f.taylor(x, 0, 7)

\[ atan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} \tag{3.1} \]

In the graph “atan function plot” 3.1, we can see \(atan(x)\) is an odd function. Which is also something we can notice from its Taylor series. Therefore, the \(atan\) approximation will only have odd \(a_k\) terms, e.g. \(a_1, a_3, ~\dots\). Which, by the way, correspond to even polynomials \(U_0, U_2, \dots\). Another useful feature from the Taylor series is the alternating sign of the terms.

Figure 3.1: atan function plot. Here we can see \(atan(-x)=-atan(x)\); it is an odd function. Therefore \(a_k(atan)\) for even k is zero.

Using our (2.5) equation for the approximation from the function derivative, we have:

\[ \begin{align} \frac{d(atan(x))}{dx} &= \frac{1}{1+x^2} \\[3pt] \implies \\ a_k(atan) &= \frac{2}{k\,\pi} \int_{-1}^1 U_{k-1}(x)·\frac{\sqrt{1-x^2}}{1+x^2} \,dx \tag{3.2} \end{align} \]

3.1 First \(a_k\) Coefficients

The first terms of our approximation, considering only the odd terms because the others are zero, are:

\[ \begin{align} a_1(atan) &= \frac{2}{\pi} \int_{-1}^1 \frac{\sqrt{1-x^2}}{1+x^2} \,dx \\ a_3(atan) &= \frac{2}{3\,\pi} \int_{-1}^1 (4\,x^2-1)·\frac{\sqrt{1-x^2}}{1+x^2} \,dx \\ a_5(atan) &= \frac{2}{5\,\pi} \int_{-1}^1 (16\,x^4-12\,x^2+1)·\frac{\sqrt{1-x^2}}{1+x^2} \,dx \\ \end{align} \]

Fortunately, nowadays, we can count on online symbolic integrators like in sageMath to evaluate those integrals. For example, for the \(a_5(atan)\) coefficient, we use:

k = 5
x = var('x'); f(x) = (16*x^4-12*x^2+1)*sqrt(1-x^2)/(1+x^2)
a_5 = integrate(f, x, -1, 1)*2/(k*pi)
a_5

At the end we found:

\[ \begin{align} a_1(atan) &= 2·\sqrt{2} - 2 &&\approx \phantom{+}0.8284271\\ a_3(atan) &= -10/3·\sqrt{2} + 14/3 &&\approx -0.0473785\\ a_5(atan) &= 58/5·\sqrt{2} - 82/5 &&\approx \phantom{+}0.0048773\\ a_7(atan) &\approx -0.0005977 \end{align} \]

3.2 Finding a Generalized \(a_k\) Formula

At this point, manipulating the equation (3.2) with trigonometric identities –in the past– I got a generalization for \(a_k(atan)\), which is too long and tedious to put it here. Instead, we shall get another route to find it. However, somehow, I am “cheating” here because I already know it works.

Let us say, we suspect that for some \(\alpha\) and \(\beta\) constants, the \(a_{n+2}, a_n\), and \(a_{n-2}\) coefficients follow a linear relationship:

\[ (n+2)·a_{n+2} = \alpha·n·a_n + \beta·(n-2)·a_{n-2} \]

To simplify the notation, and also the following steps, let us say that \(a_k^*=\lvert k·a_{k}\rvert\). Where \(\lvert x\rvert\) is the absolute value of \(x\). Then, the previous relationship becomes:

\[ (a_{n+2}^*) = \alpha·(a_n^*) + \beta·(a_{n-2}^*) \]

To find \(\alpha\) and \(\beta\), we will use the values \(a_5,a_3,a_1\) we found. But as the relationship is for the absolute values we will invert the sign of the \(a_3^*\) value, which is negative.

\[ \begin{align} 58\,\sqrt{2} - 82 &= \alpha·(10·\sqrt{2} - 14) + \beta·(2\,\sqrt{2} - 2) \\ 58\,\sqrt{2} &= \phantom{-}10\,\alpha\,\sqrt{2} + 2\,\beta\,\sqrt{2} \\ -82 &= -14\,\alpha - 2\,\beta \\ \implies &\alpha = \phantom{-}6 \\ &\beta = -1 \\ \implies a_{n+2}^* &= 6\,a_n^* -a_{n-2}^* \\ \therefore ~a_{n+2}^* &- 6\,a_n^* +a_{n-2}^* = 0 \tag{3.3} \\ \end{align} \]

Fortunately, the equation (3.3) is a Linear Difference Equation. You can read about those equations in this Wikipedia article. However, the following shorts steps are easy to understand without reading that post.

3.2.1 Characteristic Equation and Roots

To solve the (3.3) equation, we just consider it like a polynomial of some \(\lambda\) variable like this:

\[ \lambda^4 - 6\,\lambda^2 + 1 = 0 \]

If we call \(\lambda_1, \lambda_2, \dots, \lambda_4\) to the four roots of the equation above, then we can say:

\[ \begin{align} a_k^* = c_1\,\lambda_1^k + c_2\,\lambda_2^k + c_3\,\lambda_3^k + c_4\,\lambda_4^k \end{align} \]

Where \(c_1, \dots c_4\) are some constants we will find later. Let’s do it!

\[ \begin{align} \lambda^2 &= \frac{6\pm \sqrt{(-6)^2-4}}{2}=3\pm 2\sqrt{2} \\ &=(\sqrt{2}\pm 1)^2 \\[3pt] \implies \lambda_{1,2,3,4} &= \pm(\sqrt{2}\pm 1) \end{align} \]

In the equation above, we are abusing our notation to say the four roots \(\lambda_{1,2,3,4}\) are equal to the four combinations of signs in the right-side expression.

Finding the Solution Coefficients

With the roots we just found, we can say:

\[ a_k^* = c_1\,(\sqrt{2}+1)^k + c_2\,(\sqrt{2}-1)^k + c_3\,(-(\sqrt{2}+1))^k + c_4\,(-(\sqrt{2}-1))^k \]

As in our solution \(k\) is odd, then \((-(\sqrt{2}+1))^k=-(\sqrt{2}+1)^k\). Therefore:

\[ a_k^* = c_1^*\,(\sqrt{2}+1)^k + c_2^*\,(\sqrt{2}-1)^k \]

Where \(c1^*=(c_1-c_3)\) and \(c2^*=(c_1-c_4)\). If we compare \(a_1^*\) with its value \(2\sqrt{2}-2\) we get:

\[ \begin{align} c_1^*\,(\sqrt{2}+1)^1 + c_2^*\,(\sqrt{2}-1)^1 &= 2\sqrt{2}-2 \\ \implies c_1^* &= 0 \\ c_2^* &= 2 \end{align} \]

If we do this again, but this time comparing \(a_3^*\) with its value \(10\,\sqrt{2}-14\), we get \(c_2^*=2\) again. Now we know:

\[ a_k^* = k·a_k = 2\,(\sqrt{2}-1)^k \]

3.3 Polishing the atan Approximation

One last missing thing is that, by looking at the first approximations of \(a_k\), we can see the negative sign is alternating, matching its Taylor series (3.1). The sign is \((-1)^{(k-1)/2}\). With all this, we finally have:

\[ a_k = \begin{cases} (-1)^{(k-1)/2}\, (\frac{2}{k}) \, (\sqrt{2}-1)^k & \text{if k is odd,} \\ 0 & \text{otherwise} \end{cases} \]

By using \(2\,k-1\) instead of \(k\) to get rid of the “k is odd” restriction, we get our \(atan^*\) approximation:

\[ \begin{align} atan^*(x) &= \sum_{k=1}^\infty b_k·P_k(x) ~, ~~~\text{for}~ \lvert x \rvert \leq 1 \\ b_k &= (-1)^{(k-1)}\, (\frac{2}{2\,k-1}) \, (\sqrt{2}-1)^{(2k-1)} \\[5pt] P_{0}(x) &= x \\ P_{1}(x) &= x \\ P_{k}(x) &= (4\,x^2-2)·P_{k-1} - P_{k-2} \end{align} \]

3.3.1 Extending the Approximation Domain

We can extend the domain for the \(atan\) computation by using the identity \(atan(z)+atan(1/z)=\frac{\pi}{2}sign(z)\):

\[ \begin{aligned} atan(z) &= \begin{cases} atan^*(z) & \text{if } \lvert x \rvert \leq 1, \\[4pt] \frac{\pi}{2}\,sign(z)-atan^*(1/z) & \text{otherwise} \end{cases} \end{aligned} \]

One particular case is when \(x=1\), where \(atan(1)=\pi/4\) and \(T_k=1\) for every \(k\). As \(T_k(x)=cos(k\,t)\) then \(T_k(x)\) is in \([-1,1]\) and this is a challenging case when all \(T_k(x)\) has the greatest possible value. Having \(T_k=1\), we just have to sum the approximation coefficients:

\[ \begin{aligned} \frac{\pi}{4} &= 2\rho - \frac{2}{3}\rho^3 + \frac{2}{5}\rho^5 + \dots \\ \frac{\pi}{4} &= 2\,(\rho - \frac{\rho^3}{3} + \frac{\rho^5}{5} + \dots) = 2 · atan(\rho) \\[3pt] \rho &= \sqrt{2}-1 \end{aligned} \]

At the right side we recognize the \(atan\) Taylor approximation, hence \(\frac{\pi}{4}=2·atan(\sqrt{2}-1)\). In other words, the series computes the \(atan\) of the half the angle \((\pi/8)\) and then doubles the resulting angle to get \(\pi/4\), limiting the use of the Taylor series to \(\rho=\sqrt{2}-1\approx 0.414\) where it is still a good enough estimate. This behavior is key for this approximation fast convergence. The tenth term in the series is \(\rho^{19}\) which is equal to \(2.8\times 10 ^{-9}\), and these terms are decreasing fast because each one is at least \(1/(\sqrt{2}-1)^2=5.8\) times smaller than the previous. Therefore, the approximation error is small \(atan^*(1)-\pi/4=-7.5\times 10 ^{-10}\). In contrast, in the regular Taylor approximation, the tenth term is \(1/19=0.05\), and the error is big \(-0.025\).

The graph Atan error: Chebyshev vs. Taylor. 3.2 shows the Chebyshev Cheby 12 and the Taylor Taylor 12 approximation error, both using 12 terms. The Taylor approximation, by design, is perfect for values close to zero. For arguments above 0.4, the error grows steadily and very fast. In contrast, the Chebysheb polynomial has a wavy error, which is the signature of the best approximations.

The best possible approximation is the one that minimizes the maximum error. Those are characterized by having a wavy error pattern, and the Remez algorithm is a method to find it. In the graph 3.2, the dashed vertical lines show the points where the Chebyshev approximation has a local extreme value, and the Remez algorithm works to reduce the maximum error in those points.

After one Remez iteration, we got the Remez Cheb 12.1 approximation, shown in the plot. The ratio between the maximum and minimum absolute errors in local maxima points is reduced, from 1.36 in Chebyshev to 1.00 in Remez. Which means that is the minimax approximation to \(atan\). However, the maximum error is about the same; the Remez error is smaller by a factor of just 1.16. In other words, in this case, it doesn’t seem to worth the extra precision gained with the Remez iteration, because this Chebyshev approximation is almost as good as the best possible one.

Figure 3.2: Atan error: Chebyshev vs Taylor. Chebyshev and Taylor approximations with 12 terms. The dashed lines flag the points where the Chebyshev approximation makes maximum local errors. The Chebyshev approximation is very close to the minimax one. Minimax is the polynomial approximation, which minimizes the maximum error. You can tell the Chebyshev approximation is close to the minimax when all the local top errors are about the same. The ratio between the minimum and maximum absolute errors for the points marked with dashed lines is 2.56/1.86=1.37. After one iteration with the Remez algorithm (Remez Cheb 12.1), that ratio is 1.00. However, the maximum error is just a bit smaller.

References

Abramowitz, Milton, and Irene A. Stegun. 1964. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Ninth Dover printing, tenth GPO printing. New York: Dover.

Clenshaw, C. W. 1954. “Polynomial Approximations to Elementary Functions.” Mathematics of Computation 8: 143–47.

Gil, A., J. Segura, and N. M. Temme. 2007. Numerical Methods for Special Functions. Other Titles in Applied Mathematics. Society for Industrial; Applied Mathematics. https://books.google.com.pe/books?id=c5gLFYcKHSgC.

Kernighan, Brian W., and Dennis M. Ritchie. 1978. The c Programming Language. 1st ed. Prentice Hall Professional Technical Reference.

The Sage Developers. 2021. SageMath, the Sage Mathematics Software System (Version 9.0).

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Comments