/*	Copyright (c) 2004 Simon Brown

	Permission is hereby granted, free of charge, to any person obtaining a 
	copy of this software and associated documentation files (the "Software"),
	to deal in the Software without restriction, including without limitation
	the rights to use, copy, modify, merge, publish, distribute, sublicense,
	and/or sell copies of the Software, and to permit persons to whom the
	Software is furnished to do so, subject to the following conditions:

	The above copyright notice and this permission notice shall be included in
	all copies or substantial portions of the Software.

	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 
	DEALINGS IN THE SOFTWARE.
*/

#include <iostream>
#include <map>
#include <cassert>
#include <fstream>
#include <vector>

int GCDRecurse( int a, int b )
{
	assert( a >= b );
	if( b == 0 )
		return a;
	else 
		return GCDRecurse( b, a % b );
}

int GCD( int a, int b )
{
	if( a < 0 )
		a = -a;
	if( b < 0 )
		b = -b;

	if( b > a )
		std::swap( a, b );

	assert( a != 0 );
	return GCDRecurse( a, b );
}

int LCM( int a, int b )
{
	return std::abs( a*b )/GCD( a, b );
}

class Rational
{
public:
	Rational() : m_num( 0 ), m_denom( 1 ) {}

	explicit Rational( int num, int denom =  1 ) : m_num( 0 ), m_denom( 1 )
	{
		Assign( num, denom );
	}

	int GetNumerator() const { return m_num; }
	int GetDenominator() const { return m_denom; }

	Rational operator-() const
	{
		return Rational( -m_num, m_denom );
	}
	
	Rational& operator+=( Rational const& arg )
	{
		Assign( arg.m_denom*m_num + m_denom*arg.m_num, m_denom*arg.m_denom );
		return *this;
	}
	
	Rational& operator-=( Rational const& arg )
	{
		Assign( arg.m_denom*m_num - m_denom*arg.m_num, m_denom*arg.m_denom );
		return *this;
	}

	Rational& operator*=( Rational const& arg )
	{
		Assign( m_num*arg.m_num, m_denom*arg.m_denom );
		return *this;
	}
	
	Rational& operator/=( Rational const& arg )
	{
		Assign( m_num*arg.m_denom, m_denom*arg.m_num );
		return *this;
	}

	bool operator<( Rational const& arg ) const
	{
		return m_num*arg.m_denom < arg.m_num*m_denom;
	}

	bool operator==( Rational const& arg ) const
	{
		return m_num == arg.m_num && m_denom == arg.m_denom;
	}
	
private:
	void Assign( int num, int denom )
	{
		m_num = num;
		m_denom = denom;

		assert( m_denom != 0 );
		if( m_denom < 0 )
		{
			m_num = -m_num;
			m_denom = -m_denom;
		}

		int gcd = GCD( m_num, m_denom );
		m_num /= gcd;
		m_denom /= gcd;			
	}
	
	int m_num;
	int m_denom;
};

Rational operator+( Rational const& left, Rational const& right )
{
	Rational copy( left );
	copy += right;
	return copy;
}

Rational operator-( Rational const& left, Rational const& right )
{
	Rational copy( left );
	copy -= right;
	return copy;
}

Rational operator*( Rational const& left, Rational const& right )
{
	Rational copy( left );
	copy *= right;
	return copy;
}

Rational operator/( Rational const& left, Rational const& right )
{
	Rational copy( left );
	copy /= right;
	return copy;
}

std::ostream& operator<<( std::ostream& ostr, Rational const& arg )
{
	if( arg.GetDenominator() != 1 )
		ostr << "\\frac{" << arg.GetNumerator() << "}{" << arg.GetDenominator() << "}";
	else
		ostr << arg.GetNumerator();
	return ostr;	
}

class SpherePower
{
public:
	SpherePower() 
	{
		m_coeffs[0] = 0;
		m_coeffs[1] = 0;
		m_coeffs[2] = 0;
		m_coeffs[3] = 0;
	}

	SpherePower( int a, int b, int c, int d )
	{
		m_coeffs[0] = a;
		m_coeffs[1] = b;
		m_coeffs[2] = c;
		m_coeffs[3] = d;
	}

	int operator[]( int index ) const 
	{ 
		assert( 0 <= index && index < 4 );
		return m_coeffs[ index ]; 
	}

	int& operator[]( int index )
	{ 
		assert( 0 <= index && index < 4 );
		return m_coeffs[ index ]; 
	}

	SpherePower& operator*=( SpherePower const& arg )
	{
		for( int i = 0; i < 4; ++i )
			m_coeffs[i] += arg.m_coeffs[i];
		return *this;
	}

	bool operator<( SpherePower const& arg ) const
	{
		for( int i = 0; i < 4; ++i )
		{
			if( m_coeffs[i] < arg.m_coeffs[i] )
				return true;
			if( m_coeffs[i] > arg.m_coeffs[i] )
				return false;
		}
		return false;
	}

private:
	int m_coeffs[4];	
};

template< typename C, typename P >
class Polynomial
{
private:
	typedef std::map< P, C > Coeffs;
	typedef typename Coeffs::const_iterator ConstIterator;
	typedef typename Coeffs::iterator Iterator;
	
public:
	C operator[]( P const& power ) const 
	{
		ConstIterator it = m_coeffs.find( power );
		if( it != m_coeffs.end() )
			return it->second;
		return C();
	}
	
	C& operator[]( P const& power )
	{
		Iterator it = m_coeffs.find( power );
		if( it == m_coeffs.end() )
			it = m_coeffs.insert( std::make_pair( power, C() ) ).first;
		return it->second;
	}

	Polynomial< C, P >& operator+=( Polynomial< C, P > const& arg )
	{
		for( ConstIterator src = arg.m_coeffs.begin();
			 src != arg.m_coeffs.end(); ++src )
		{
			if( !m_coeffs.insert( *src ).second )
				m_coeffs[ src->first ] += src->second;
		}
		return *this;
	}
	
	Polynomial< C, P >& operator-=( Polynomial< C, P > const& arg )
	{
		for( ConstIterator src = arg.m_coeffs.begin();
			 src != arg.m_coeffs.end(); ++src )
		{
			if( !m_coeffs.insert( std::make_pair( src->first, -src->second ) ).second )
				m_coeffs[ src->first ] -= src->second;
		}
		return *this;
	}
	
	Polynomial< C, P >& operator*=( C const& arg )
	{
		for( Iterator it = m_coeffs.begin(); it != m_coeffs.end(); ++it )
			it->second *= arg;
		return *this;
	}
	
	Polynomial< C, P >& operator*=( P const& arg )
	{
		Coeffs coeffs;
		for( ConstIterator it = m_coeffs.begin(); it != m_coeffs.end(); ++it )
		{
			coeffs.insert( std::make_pair( it->first * arg, it->second ) );
		}
		m_coeffs.swap( coeffs );
		return *this;
	}

	Polynomial< C, P >& operator*=( Polynomial< C, P > const& arg )
	{
		Polynomial< C, P > result;
		for( ConstIterator it = m_coeffs.begin(); it != m_coeffs.end(); ++it )
		{
			Polynomial< C, P > copy( arg );
			copy *= it->first;
			copy *= it->second;

			result += copy;
		}
		m_coeffs.swap( result.m_coeffs );
		return *this;
	}
	
	friend std::ostream& operator<<( std::ostream& ostr, Polynomial< C, P > const& arg )
	{
		for( ConstIterator it = arg.m_coeffs.begin();; )
		{
			bool zero = ( it->second == C() );
			if( !zero )
				ostr << it->second << it->first;
			if( ++it == arg.m_coeffs.end() )
				break;
			if( !zero && C() < it->second )
				ostr << " + ";
		}
		return ostr;
	}

	template< typename InputIterator >
	void GetPowers( InputIterator input ) const
	{
		for( ConstIterator it = m_coeffs.begin(); it != m_coeffs.end(); ++it )
			*input++ = it->first;
	}

private:
	Coeffs m_coeffs;
};

template< typename C, typename P >
Polynomial< C, P > operator+( Polynomial< C, P > const& left, 
							  Polynomial< C, P > const& right )
{
	Polynomial< C, P > copy( left );
	copy += right;
	return copy;
}

template< typename C, typename P >
Polynomial< C, P > operator-( Polynomial< C, P > const& left, 
							  Polynomial< C, P > const& right )
{
	Polynomial< C, P > copy( left );
	copy -= right;
	return copy;
}

template< typename C, typename P >
Polynomial< C, P > operator*( Polynomial< C, P > const& left, 
							  Polynomial< C, P > const& right )
{
	Polynomial< C, P > copy( left );
	copy *= right;
	return copy;
}

SpherePower operator*( SpherePower const& left, SpherePower const& right )
{
	SpherePower copy( left );
	copy *= right;
	return copy;
}

std::ostream& operator<<( std::ostream& ostr, SpherePower const& index )
{
	char const* prefixes[] = { "\\cos", "\\sin", "\\cos", "\\sin" };
	char const* postfixes[] = { "\\phi", "\\phi", "\\theta", "\\theta" };

	for( int i = 0; i < 4; ++i )
	{	
		if( index[i] == 1 )
			ostr << prefixes[i] << postfixes[i];
		else if( index[i] > 1 )
			ostr << prefixes[i] << "^" << index[i] << postfixes[i];
	}
	return ostr;
}

#define CONSTANT	SpherePower()
#define COS_PHI		SpherePower( 1, 0, 0, 0 )
#define SIN_PHI		SpherePower( 0, 1, 0, 0 )
#define COS_THETA	SpherePower( 0, 0, 1, 0 )
#define SIN_THETA	SpherePower( 0, 0, 0, 1 )

typedef Polynomial< Rational, SpherePower > TrigPolynomial;

TrigPolynomial P( int l, int m )
{
	assert( 0 <= m && m <= l );
	TrigPolynomial result;
	if( l == 0 && m == 0 )
	{
		result[ CONSTANT ] = Rational( 1 );
	}
	else if( l == m )
	{
		result = P( m - 1, m - 1 );
		result *= SIN_THETA;
		result *= Rational( 2*m - 1 );
	}
	else if( l == m + 1 )
	{
		result = P( m, m );
		result *= COS_THETA;
		result *= Rational( 2*m + 1 );
	}
	else
	{
		TrigPolynomial left = P( l - 1, m );
		left *= COS_THETA;
		left *= Rational( 2*l - 1 );
			
		TrigPolynomial right = P( l - 2, m );
		right *= Rational( l + m - 1 );

		result = left - right;
		result *= Rational( 1, l - m );
	}
	return result;
}

TrigPolynomial SinMPhi( int m );
TrigPolynomial CosMPhi( int m );
	
TrigPolynomial Q( int m )
{
	TrigPolynomial result;
	if( m == 0 )
	{
		result[ CONSTANT ] = Rational( 1 );
	}
	else if( m < 0 )
	{
		result = SinMPhi( -m );		
	}
	else // m > 0
	{
		result = CosMPhi( m );
	}
	return result;
}

TrigPolynomial SinMPhi( int m )
{
	assert( m > 0 );
	TrigPolynomial result;
	if( m == 1 )
	{
		result[ SIN_PHI ] = Rational( 1 );
	}
	else
	{
		TrigPolynomial left = CosMPhi( m - 1 );
		left *= SIN_PHI;

		TrigPolynomial right = SinMPhi( m - 1 );
		right *= COS_PHI;

		result = left + right;
	}
	return result;
}

TrigPolynomial CosMPhi( int m )
{
	assert( m > 0 );
	TrigPolynomial result;
	if( m == 1 )
	{
		result[ COS_PHI ] = Rational( 1 );
	}
	else
	{
		TrigPolynomial left = CosMPhi( m - 1 );
		left *= COS_PHI;

		TrigPolynomial right = SinMPhi( m - 1 );
		right *= SIN_PHI;

		result = left - right;
	}
	return result;
}

int Factorial( int i )
{
	assert( i >= 0 );
	return ( i < 2 ) ? 1 : ( i*Factorial( i - 1 ) );
}

int LargestSquareFactor( int value )
{
	int factor = 1;
	for( int i = 2; i*i <= value; ++i )
	{
		int const ii = i*i;
		while( ( value % ii ) == 0 )
		{
			factor *= i;
			value /= ii;
		}			
	}
	return factor;
}

int main()
{
	static int const order_count = 4;
	
	std::ofstream file( "output.tex" );
	file << "\\documentclass{article}\n" 
		<< "\\pagestyle{empty}\n"
		<< "\\begin{document}\n\n";
	
	for( int l = 0; l <= order_count; ++l )
	{
		for( int m = -l; m <= l; ++m )
		{
			int abs_m = std::abs( m );
			
			// compute the components
			TrigPolynomial p = P( l, abs_m );
			TrigPolynomial q = Q( m );
			TrigPolynomial y = p * q;

			int qscale = ( m == 0 ) ? 1 : 2;			
			Rational n( 
				( 2*l + 1 )*Factorial( l - abs_m )*qscale, 
				4*Factorial( l + abs_m ) 
			);

			// simplify the coefficients
			typedef std::vector< SpherePower > Powers;
			Powers powers;
			y.GetPowers( std::back_inserter( powers ) );

			// get the LCM of the denominators (so as to remove them)
			int lcm = 1;
			for( Powers::const_iterator it = powers.begin(); it != powers.end(); ++it )
				lcm = LCM( lcm, y[*it].GetDenominator() );

			n *= Rational( 1, lcm*lcm );
			for( Powers::const_iterator it = powers.begin(); it != powers.end(); ++it )
				y[*it] *= Rational( lcm );

			// get the GCD of the numerators (so as to simplify them)
			int gcd = 0;
			for( Powers::const_iterator it = powers.begin(); it != powers.end(); ++it )
				gcd = GCD( gcd, y[*it].GetNumerator() );
						
			n *= Rational( gcd*gcd, 1 );
			for( Powers::const_iterator it = powers.begin(); it != powers.end(); ++it )
				y[*it] *= Rational( 1, gcd );
			
			// simplify the squareroot term
			Rational n2( LargestSquareFactor( n.GetNumerator() ), 
						 LargestSquareFactor( n.GetDenominator() ) );
			n /= n2*n2;

			// print it out
			file << "$Y_" << l << "^{" << m << "} = " << n2 
				<< "\\frac{\\sqrt{" << n.GetNumerator() <<"}}{\\sqrt{";
			if( n.GetDenominator() != 1 )
				file << n.GetDenominator();
			file << "\\pi}}(" << y << ")$\n\n";
		}
	}

	file << "\\end{document}\n";
	file.close();

	return 0;
}