/*********************************************************************
   Computes Legendre Polynomials P_l^m(x), where m and l are integers
   satisfying 0<= m <= l, while x lies in the range -1 <= x <= l.
   C.A. Bertulani        May/23/2000
*********************************************************************/
typedef double Number;

#include <math.h>
#include<iostream.h>

int main(){
	Number plgndr(int, int, Number);
	Number x, pl, p00, p10, p11, p20, p21, p22;
	int l,m,j;
	j=1;
	cout << "\n\n Enter x.\n";
	cin >> x;
    for(;;){
		if(j>5){ cout << "\n My patience is over. Stop, please!\n";
			break;
		}
		if(j!=1){
			cout << "\n\n Enter x.\n";
			cin >> x;
		}
		for(l=0;l<=2;l++){
			for(m=0;m<=l;m++){
				pl = plgndr(l, m, x);
				p00=1;
				p11=-sqrt(1-x*x);
				p10=x;
				p22=3*(1-x*x);
				p21=-3*sqrt(1-x*x)*x;
				p20=0.5*(3.*x*x-1);
				cout << "P_"<<l<<"^"<<m<<"(x)=" << pl<<endl ;

			}
		}
		cout << "\nCompare with exact values:\n" << p00 
			<<"\t" << p10 << '\t' << p11 << '\t' << p20 << '\t' 
			 << p21 << '\t' << p22 << endl; 
		j=j+1;
	}
	return 0;
}
/*********************************************************************
   Returns Legendre Polynomials P_l^m(x), where m and l are integers
   satisfying 0<= m <= l, while x lies in the range -1 <= x <= l.
   C.A. Bertulani        May/23/2000
*********************************************************************/
Number plgndr(int l, int m, Number x)
{
	Number fact,pll,pmm,pmmp1,somx2;
	int i,ll;

	try{
		if (m < 0 || m > l || fabs(x) > 1.0){
			throw "Bad arguments in routine plgndr\n";
		}
	}
	catch (char* message) {
          cerr << message;return 0.;
	}
	pmm=1.0;   /*  Compute P_m^m */
	if (m > 0) {
		somx2=sqrt((1.0-x)*(1.0+x));
		fact=1.0;
		for (i=1;i<=m;i++) {
			pmm *= -fact*somx2;
			fact += 2.0;
		}
	}
	if (l == m)
		return pmm;
	else {			   /*  Compute P_(m+1)^m */
		pmmp1=x*(2*m+1)*pmm;
		if (l == (m+1))
			return pmmp1;
		else {			   /*  Compute P_l^m, l > m+1 */
			for (ll=m+2;ll<=l;ll++) {
				pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
				pmm=pmmp1;
				pmmp1=pll;
			}
			return pll;
		}
	}
}
/*********************************************************************/