Error writing recursive function for Bonnet's recursion formula in C

[PrakashPhy]

Bonnet’s recursion formula for Legendre polynomials is
P(n,x)=1 for n=0
P(n,x)=x for n=1 and
P(n,x)= (2n-1)/n*x*P(n-1,x)-(n-1)/n*P(n-2,x)

I tried to write recursive function in C to calculate the value of polynomial for a given n and x

float legpol(int n, float x)
    {
        if(n==0)
            return 1;
        else if(n==1)
            return x;
        else
            return(2*n-1)/n*x*legpol(n-1,x)-(n-1)/n*legpol(n-2,x);
    }

But this doesn't give the correct result. What's wrong with the function??

Thanks in advance

 

[bcbwilla]

Can you give a sample of the output vs. what the correct answer should be?

 

[PrakashPhy]

The actual legendre's polynomials are
for n = 0 1
for n= 1 x
for n = 2

for n= 3

for n =4

and so onif n= 2 and x= .1 the output must be -0.485 but my recursive function gives 0.01 which is wrong

 

[Mark44]

The problem is, I believe, that your formula is using int division in a couple of places, and they give you incorrect function values. Inexperienced C and C++ programmers often don't realize that there are two kinds of division - integer and floating point.

Your code has this line:

return(2*n-1)/n*x*legpol(n-1,x)-(n-1)/n*legpol(n-2,x);

If n == 3, for example, the (2*n - 1)/n expression evaluates to 5/3 == 1. The other int term is multiplied by 2/3 == 0.

An easy fix is to change the line above to

return(2*n-1.0)/n*x*legpol(n-1,x)-(n-1.0)/n*legpol(n-2,x);

This change causes the numerators of both fractions to be promoted to double before the division actually occurs, forcing floating point division to be performed, rather than integer division.

 

[PrakashPhy]

Thank you! This change worked fine.

I tried to figure out the error and checked the value of 'n' after each function call and thought there might have been changes in the value of 'n' within the recursive function call. So i tried to preserve the value of 'n' by creating another variable 'b' just to use it in nested function call. I modified the above function like this::

float legpol(int n, float x)
    {
        int b;
        if(n==0)
            return 1;
        else if(n==1)
            return x;
        else
        {
             b=n;
             return(2*b-1)/b*x*legpol(b-1,x)-(b-1)/b*legpol(b-2,x);
         }
     }

Amazingly it worked. In this case why shouldn't i write (2*b-1.0) ??

 

[Mark44]

I don't know why your revision should work. For some int divisions you get the answer you would expect, such as 6/3 or 10/2. It might be that what got evaluated was something like this.

In any case, I would get rid of the b variable and use code that forces floating point division, such as 2nd example in my previous post.

 

[PrakashPhy]

I'm sorry for the second revision i posted. I had actually typed

float b;

in the third line of my function not

 int b

.

I now know what the problem was. Thank you very much for the solution. It really helped me.

 

[Mark44]

I'm sorry for the second revision i posted. I had actually typed

float b;

in the third line of my function not

 int b

.

I now know what the problem was. Thank you very much for the solution. It really helped me.

Yes, that would have made a difference, by causing floating point division to be performed instead of integer division.