Taylor's Series by Horner's method using recursion

Question

Our professor has taught us to find Taylor's series by recursion using a static summation variable that accumulates the summation of each term and then multiples with itself with each iteration.

#include <stdio.h>

double t(int n,int m){
    static double a = 1.0;
    if(m==0)
        return a;
    a = 1+(double)n/m*a;
    return t(n,m-1);
}

int main(){
    printf("%lf\n",t(4,1500));
}

This, yields the desired accuracy and result.

---

I tried to implement the code using a single return statement:

#include <stdio.h>

double e(int x, int m){
    if(m==0)
        return 1.0;
    return 1+((double)x/m)*e(x,m-1);
}

int main(){
    printf("%lf\n",e(4,1500));
}

This, yields an output of 1.002674.

What is the difference between the two? What causes this code to lead this inaccuracy?

Answer 1

• List item

This is not a matter of an inaccuracy. The expression that is being calculated is different.

Let's replace the second argument (1500) by just 3. So t(4, 3) and e(4, 3). The logic of t(4, 3) calculates this expression:

a = 1
a = 1 + (4 / 3) * a
a = 1 + (4 / 2) * a
a = 1 + (4 / 1) * a

This expands to:

1 + (4 / 1) * (1 + (4 / 2) * (1 + (4 / 3) * 1))
= 23.66666...

while e(4, 3) is calculated like this:

e(4, 0) = 1
e(4, 1) = 1 + (4 / 1) * e(4, 0)
e(4, 2) = 1 + (4 / 2) * e(4, 1)
e(4, 3) = 1 + (4 / 3) * e(4, 2)

which expands to:

1 + (4 / 3) * (1 + (4 / 2) * (1 + (4 / 1) * 1))
=15.66666...

With this small example, you can see how the fractions are calculated in the opposite order. Because of the 1+ operation at each intermediate result, a different order yields a different result.

It might be important to realise that in the first version the variable a has reached its final value by the time the base case hits.

Correction

I can understand you'd want to refactor the first function, as the use of a static variable is going to turn sour when the function is called a second time at the top level. In that case you'd want that a to reset to 1, but it doesn't.

You could replace it with an extra parameter which you pass the starting value of 1.0. Not an issue, but you can omit the explicit cast to double by putting a as first operand of the expression that needs a floating point division. I named the function f:

double f_recur(int n, int m, double a) {
    if (m == 0)
        return a;
    return f_recur(n, m - 1, 1 + a * n / m);
}

double f(int n, int m) {
    return f_recur(n, m, 1.0);
}

And you could also do this in an iterative way, which is preferable as it doesn't use that additional stack space. This time called g:

double g(int n, int m) {
    double a = 1.0;
    for (int i = m; i > 0; i--) {
        a = 1 + a * n / i;
    }
    return a;
}