Own asin() function (with Taylor series) not accurate

Demystifying Inaccurate asin() Implementations with Taylor Series

Explore the challenges of implementing an accurate asin() function using Taylor series, understand common pitfalls, and discover strategies for improving precision in C.

Implementing mathematical functions from scratch, especially inverse trigonometric functions like asin(), can be a fascinating and educational exercise. Often, the Taylor series expansion is the first method that comes to mind due to its direct mathematical definition. However, as many developers discover, a direct Taylor series implementation for asin(x) can quickly lead to significant inaccuracies, particularly as x approaches its domain boundaries of -1 or 1. This article delves into why this happens and what can be done to achieve better precision.

The Taylor Series for asin(x)

The Taylor series expansion for asin(x) around x = 0 is given by:

asin(x) = x + (1/2) * (x^3/3) + (1*3)/(2*4) * (x^5/5) + (1*3*5)/(2*4*6) * (x^7/7) + ...

This can be written more compactly using factorials and powers:

asin(x) = sum_{n=0}^{inf} ( (2n)! / ( (4^n) * (n!)^2 * (2n+1) ) ) * x^(2n+1)

While mathematically sound, this series converges slowly, especially for values of x far from 0. The number of terms required to achieve a certain level of precision grows dramatically as |x| approaches 1. This slow convergence is the primary culprit behind the inaccuracies observed in simple implementations.

flowchart TD
    A[Input x] --> B{Is |x| > 0.5?}
    B -- Yes --> C[Apply identity: asin(x) = pi/2 - asin(sqrt(1-x^2))]
    B -- No --> D[Compute asin(x) using Taylor series]
    C --> E[Return Result]
    D --> E

Decision flow for improving asin(x) accuracy

Common Pitfalls and Sources of Inaccuracy

Several factors contribute to the inaccuracy of a naive Taylor series asin() implementation:

1. Slow Convergence: As mentioned, the series converges very slowly for |x| close to 1. This means you need an extremely large number of terms to get reasonable precision, which is computationally expensive and can lead to floating-point precision issues.
2. Floating-Point Precision Limits: Even with double precision, there's a limit to how accurately numbers can be represented. Summing many small terms can accumulate rounding errors, especially if terms become very small relative to the running sum.
3. Cancellation Error: When x is very small, x^3/3 and subsequent terms become much smaller than x. If x itself is very close to a representable floating-point number, subtracting or adding these tiny terms might not change the value due to the limited precision of the mantissa.
4. Domain Handling: The Taylor series is valid for |x| <= 1. Incorrect handling of inputs outside this range or at the boundaries can lead to NaN or inf results.

#include <stdio.h>
#include <math.h>

// Naive Taylor series implementation for asin(x)
double my_asin_taylor(double x, int terms) {
    if (x < -1.0 || x > 1.0) {
        return NAN; // Domain error
    }
    if (x == 0.0) return 0.0;
    if (x == 1.0) return M_PI_2; // pi/2
    if (x == -1.0) return -M_PI_2; // -pi/2

    double result = x;
    double term = x;
    double num_coeff = 1.0;
    double den_coeff = 2.0;

    for (int n = 1; n < terms; ++n) {
        num_coeff *= (2.0 * n - 1.0);
        den_coeff *= (2.0 * n);
        term *= (x * x);
        result += (num_coeff / den_coeff) * (term / (2.0 * n + 1.0));
    }
    return result;
}

int main() {
    double x_values[] = {0.1, 0.5, 0.8, 0.9, 0.95, 0.99, 0.999, 1.0};
    int num_terms = 20; // Number of terms to use

    printf("x\t\tmy_asin_taylor(%d terms)\tstd_asin()\tDifference\n", num_terms);
    printf("--------------------------------------------------------------------------------\n");

    for (int i = 0; i < sizeof(x_values) / sizeof(x_values[0]); ++i) {
        double x = x_values[i];
        double my_val = my_asin_taylor(x, num_terms);
        double std_val = asin(x);
        printf("%.3f\t\t%.10f\t\t%.10f\t%.10f\n", x, my_val, std_val, fabs(my_val - std_val));
    }

    return 0;
}

A naive C implementation of asin(x) using a fixed number of Taylor series terms. Note the M_PI_2 constant from <math.h>.

Strategies for Improving Accuracy

To overcome the limitations of a direct Taylor series expansion, professional math libraries employ several techniques:

1. Range Reduction: The most crucial technique is to reduce the input x to a smaller, more manageable range where the series converges quickly. For asin(x), this often involves using trigonometric identities. For example, if x is close to 1, asin(x) is close to pi/2. We can use the identity asin(x) = pi/2 - asin(sqrt(1-x^2)) when x > 0.5. The argument sqrt(1-x^2) will be much smaller, allowing the series to converge faster.
2. Chebyshev Polynomials or Pade Approximants: Instead of Taylor series, these methods provide polynomial or rational function approximations that are optimized for uniform accuracy across a given interval, rather than just at a single point (like Taylor series). They are generally more efficient for achieving high precision over a wider range.
3. Lookup Tables: For very specific ranges or for initial approximations, precomputed values in lookup tables can be used, followed by interpolation or a small number of series terms for refinement.
4. Iterative Methods: Some functions can be computed using iterative algorithms that refine an approximation until a desired precision is met.
5. Extended Precision Arithmetic: In cases where double precision is insufficient, some libraries might use custom extended-precision arithmetic types, though this comes with a significant performance cost.

#include <stdio.h>
#include <math.h>

// Improved asin(x) using range reduction for better accuracy
double my_asin_improved(double x, int terms) {
    if (x < -1.0 || x > 1.0) {
        return NAN; // Domain error
    }
    if (x == 0.0) return 0.0;
    if (x == 1.0) return M_PI_2;
    if (x == -1.0) return -M_PI_2;

    double abs_x = fabs(x);
    double result;
    int sign = (x < 0) ? -1 : 1;

    if (abs_x > 0.70710678118) { // Approximately sin(pi/4) or 1/sqrt(2)
        // Use identity: asin(x) = pi/2 - asin(sqrt(1-x^2))
        // This reduces the argument to a smaller value, closer to 0
        double y = sqrt(1.0 - abs_x * abs_x);
        result = M_PI_2 - my_asin_taylor(y, terms); // Reuse the Taylor series for the smaller argument
    } else {
        result = my_asin_taylor(abs_x, terms);
    }

    return sign * result;
}

int main() {
    double x_values[] = {0.1, 0.5, 0.8, 0.9, 0.95, 0.99, 0.999, 1.0};
    int num_terms = 20; // Number of terms to use

    printf("x\t\tmy_asin_improved(%d terms)\tstd_asin()\tDifference\n", num_terms);
    printf("------------------------------------------------------------------------------------\n");

    for (int i = 0; i < sizeof(x_values) / sizeof(x_values[0]); ++i) {
        double x = x_values[i];
        double my_val = my_asin_improved(x, num_terms);
        double std_val = asin(x);
        printf("%.3f\t\t%.10f\t\t%.10f\t%.10f\n", x, my_val, std_val, fabs(my_val - std_val));
    }

    return 0;
}

An improved asin(x) implementation incorporating range reduction for better accuracy, especially for x values closer to 1.