Approximating library functions

Library of mathematical functions

https://docs.python.org/3/library/math.html

• Number-theoretic and representation functions
• Power and logarithmic functions
• Angular conversions
• Trigonometric functions
• Hyperbolic functions
• Special functions
• Constants

Uses of functions

Usage examples

• log — decision trees (machine learning)
• tanh — REINFORCE algorithm (planning)
• erf — materials discovery and design (statistical modeling)

Computing mathematical functions

Example: sin and cos, in Go.

var _cos = [...]float64{
    -1.13585365213876817300e-11, // 0xbda8fa49a0861a9b
    2.08757008419747316778e-9,   // 0x3e21ee9d7b4e3f05
    -2.75573141792967388112e-7,  // 0xbe927e4f7eac4bc6
    2.48015872888517045348e-5,   // 0x3efa01a019c844f5
    -1.38888888888730564116e-3,  // 0xbf56c16c16c14f91
    4.16666666666665929218e-2,   // 0x3fa555555555554b
}

func Cos(x float64) float64 {
    const (
        PI4A = 7.85398125648498535156e-1  // 0x3fe921fb40000000, Pi/4 split into three parts
        PI4B = 3.77489470793079817668e-8  // 0x3e64442d00000000,
        PI4C = 2.69515142907905952645e-15 // 0x3ce8469898cc5170,
    )
    // special cases
    switch {
    case IsNaN(x) || IsInf(x, 0):
        return NaN()
    }

    // make argument positive
    sign := false
    x = Abs(x)

    var j uint64
    var y, z float64
    if x >= reduceThreshold {
        j, z = trigReduce(x)
    } else {
        j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
        y = float64(j)           // integer part of x/(Pi/4), as float

        // map zeros to origin
        if j&1 == 1 {
            j++
            y++
        }
        j &= 7                               // octant modulo 2Pi radians (360 degrees)
        z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
    }

    if j > 3 {
        j -= 4
        sign = !sign
    }
    if j > 1 {
        sign = !sign
    }

    zz := z * z
    if j == 1 || j == 2 {
        y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
    } else {
        y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
    }
    if sign {
        y = -y
    }
    return y
}
                        

Approximation theory

• Polynomials
• Analytic functions
• Taylor series
• Remez algorithm

Polynomials

• Polynomial: $p(x) = \sum_{i=0}^N a_i x^i$
• Polynomials ar easy to compute:

  
                                  def poly(x, aa):
                                      xi, v = 1, 0
                                      for a in aa:
                                          v += a*xi
                                          xi *= x
                                      return v
  								

• Sum and product of polynomials are polynomials
• Ratios of polynomials are also easy to compute
• $N=\infty$ ⇒ power series: $\sum_{i=0}^\infty a_i x^i$

Analytic functions

Function $f(x)$ is analytic if for any $x_0 \in D$: $$f(x) = \sum_{i=0}^\infty a_i (x-x_0)^i$$.

• $f(x)$ is analytic if can be computed as a power series.
• Polynomials are analytic.
• We can approximate $f(x) \approx \sum_{i=0}^N a_i (x-x_0)^i$, because if $|x-x_0|<1$, $\lim_{n\to\infty} (x-x_0)^n = 0$.

Examples of analytic functions

• $\exp(x) = \sum_{i=0}^\infty \frac {x^i} {i!}$
• $\ln (x) = \sum_{i=1}^\infty (-1)^{i+1}\frac{(x-1)^i}{i}$ for $0 < x < 2$
• $\sin(x) = \sum_{i=0}^\infty \frac {(-1)^i} {(2i+1)!}x^{2i+1}$

These are Taylor series!

Taylor series

If $f(x)$ is infinitely differentiable, then $$\mathrm{Taylor}(f(x), a) = \sum_{i=0}^\infty \frac {f^{(i)}(a)} {i!} (x-a)^i$$ at point $a$.

• Taylor series of a polynomial is the polynomial.
• For most functions, $f(x) = \mathrm{Taylor}(f(x), a)$ for some $a$.
• If $f(x) = \mathrm{Taylor}(f(x), a)$ for any $a$, the function is entire

Taylor series of exp(x)

• Which $N$ to take?
• What is the error?
  • $|\mathrm{Taylor}_4(exp(2), 0)-exp(2)|\approx 0.002$
  • $|\mathrm{Taylor}_4(exp(4), 0)-exp(4)|\approx 1.2$

Remez algorithm

1. Find optimal polynomial: $P(x)$, for $x_1, x_2, ..., x_{N+2}$: $$P(x_1) - f(x_1) = +\varepsilon$$ $$P(x_2) - f(x_2) = -\varepsilon$$ $$...$$
2. Move $x_i$ such that error $\le \varepsilon$ for any $x_1 \le x \le x_{N+2}$