/* calculate inner product matrix element $int_{-1}^1 Pm(x)Tn(x) dx$
 * of Legendre polynomials Pm and Chebyshev polynomials Tn
 */

#define D double
#define F float

#include <math.h>
#include "../iolib/message.h"
 
float PxT(m,n,ifnorm)
int m,n;
int ifnorm; /* ifnorm==0 -> P & T are not normalized (see ShuXueShouCe)
	       ifnorm==1 -> P & T normalized: $int_{-1}^1 Pm*Pm =1$
	       				      $int_{-1}^1 Tn*Tn =1$
	     */
{
float sum;
int N,J,K,j,k,jmax,kmax,i,nmm;
float fac,p,t;

	if(ifnorm!=0 && ifnorm!=1) STOP(PxT: illegal normalization);
	if(m<0||n<0) STOP(PxT: illegal argument);
	
	if((m+n)%2) return 0.0;
	
	jmax=n/2;
	kmax=m/2;
	
	sum=0.0;
	for(j=0;j<=jmax;j++)
	for(k=0;k<=kmax;k++)
	{ fac=1.0; if((k+j)%2) fac=-1.0; for(i=0;i<j;i++) fac/=4.0;
	
	  p=1.0; 
	  for(i=2*k;i<=m+k-1;i++) 
	  { if(i<=m-1) p*=(F)(2*m-i)/(F)(m-i);
	    else p*=(F)(2*m-i)/(F)(m+k-i);
	  }
	  
	  t=1.0;
	  for(i=1;i<=j-1;i++) t*=(F)(n-j-i)/(F)(j+1-i);
	  if(!j) { if(n) t=1./(F)n; else t=1;}

	  sum+=2.*fac*p*t/(F)(m+n-2*k-2*j+1);
	}
	
	N=n; if(!n) N=2;
	nmm=n-m-1;
	fac=1.0;
	if(nmm<0) for(i=0;i>nmm;i--) fac/=2.0;
	if(nmm>0) for(i=0;i<nmm;i++) fac*=2.0;
	sum*=fac*N;
	
	if(ifnorm) sum*=sqrt( (D)((2*m+1)*(4*n*n-1))/(D)(4*(2*n*n-1)) );
	
	return sum;
}