how to approximate hyperbolic sine using its taylor series?

Question

I'm trying to approximate hyperbolic sine. I need to do it without using math.h library functions.

[don't provide me full solutions, just some hint, because I need to figure it out by myself]

here's what I did:

given the hyperbolic sine taylor series, I need to compute the factorial of (2*n + 1). to do it, I need to do only this step:

fact *= (2*i +1); // inside a for-loop. 

I need to compute the power x^(2*n +1), and I did this way:

double power(double x, unsigned int y) {
    double result = 1;
    for (unsigned int i = 0; i < y; i++) {
        result *= x;
    }
    return result;
}

now, I have every pieces, the taylor series is implemented as follows:

#include <stdio.h> 

double power(double x, unsigned int y) {
    double result = 1;
    for (unsigned int i = 0; i < y; i++) {
        result *= x;
    }
    return result;
}

double hyp_sin(double x) {
    double result = 0;
    double fact = 1;
    double pow = 0;
    for (unsigned int i = 0; i != 21; i++) {
        fact *= (2 * i + 1);
        pow = power(x, 2 * i + 1);
        result += ((1 / fact) * pow);
    }
    return result;
}

int main(void) {
    double result = hyp_sin(89.9878);
    printf("%lf", result);
    return 0;
}

The result is completely wrong, it should have been 6.028024141598018316924203992363e+38 (with 21 iterations)

Answer 1

You don't really need to calculate powers and factorials separately for each element.
As explained in the comments these values get big very fast (especially the factorial).

Instead you can update the next element in the series incrementally:
You need to multiply by x*x (because each element adds 2 to the power of x), and divide by the next values required for the factorial.

This way you can get good precision using doubles, practically for as many iterations as you need (assuming the series converges).

This method is also a lot more efficient since for each iteration the calculations are minimal.

#include <stdio.h> 

/* Calculate the sum of the series: 
   x^1/1! + x^3/3! + x^5/5! + ... 
 */
double hyp_sin(double x, unsigned int num_iters)
{
    double res = 0;
    double elem = x; /* first element is: x^1/1! == x */
    for (unsigned int i = 1; i <= num_iters; ++i)
    {
        res += elem;
        /* update for next iteration: */
        elem = elem * x * x / (2*i) / (2*i+1);
    }
    return res;
}

int main(void) 
{
    unsigned int iterations = 100;
    double x = 89.9878;
    double result = hyp_sin(x, iterations);
    printf("%lf", result);
    return 0;
}

---

UPDATE:
Following the comment's below by @DavidC.Rankin, below you can see a version that performs as many iterations as needed until the series converges.
This is done by comparing the difference of the current element with the last one against some epsilon based on x.
You can tune CONVERGENCE_EPSILON_FACTOR according to your needs.

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

/* Calculate the sum of the series:
   x^1/1! + x^3/3! + x^5/5! + ...
 */
double hyp_sin(double x)
{
    const double CONVERGENCE_EPSILON_FACTOR = 0.0001;

    double convergence_epsilon = fabs(x * CONVERGENCE_EPSILON_FACTOR);
    double res = 0;
    double last_elem = 0;
    double elem = x; /* first element is: x^1/1! == x */
    int i = 1;
    do
    {
        res += elem;
        /* update for next iteration: */
        last_elem = elem;
        elem = elem * x * x / (2 * i) / (2 * i + 1);
        i++;
    } while (fabs(elem - last_elem) > convergence_epsilon);

    printf("hyp_sin(%f): converged after %d iterations.\n", x, i);
    return res;
}

int main(void)
{
    double x = 89.9878;
    double result = hyp_sin(x);
    printf("result: %lf\n", result);
    return 0;
}

Output:

hyp_sin(89.987800): converged after 125 iterations.
result: 602802414159797056498630200284374106112.000000