module Main  where{
import List;
import Ratio;

newtype Power_of_pi = Power_of_pi(Int) deriving (Show,Eq);
data Pure_rational = Pure_rational(Rational)(Imaginess) deriving (Show);
data Imaginess = Real | Imaginary deriving (Show);
get_number_of_pure::Pure_rational->Rational;
get_number_of_pure x = (case x of {
  (Pure_rational(a) (_))->a
});
data Poly_term = Zero_poly_term | Poly_term(Pure_rational)(Power_of_pi) deriving (Show);
type B = [](Poly_term);
zero_B::B;
zero_B = [Zero_poly_term];
one_B::B;
one_B = [(Poly_term (Pure_rational 1 Real) (Power_of_pi 0))];
type B_list = [](B);
long_all_bs::B_list;
long_all_bs = ((:) zero_B ((:) zero_B ((:) one_B (zipWith4 adder long_all_bs (tail long_all_bs) (drop 2 long_all_bs) (enumFrom 0)))));
make_Real_pure::Rational->Pure_rational;
make_Real_pure r = (Pure_rational r Real);
simpleadder::B->B->B->Int->B;
simpleadder bn_2 bn_1 bn n = [(Poly_term (make_Real_pure ((%) (toInteger n) 1)) (Power_of_pi 0))];
adder::B->B->B->Int->B;
adder bn_2 bn_1 bn n = (let {
numerat::B;
numerat = (subtract_b special bn_2);
special::B;
special = (scalar_mult (Pi_mult_B 1) (scalar_mult (I_mult_B ) (scalar_mult (Rational_mult_B ((%) (toInteger ((+) 1 ((*) 2 n))) 1)) bn)));
scals_denon::B;
scals_denon = (scalar_mult (Rational_mult_B ((%) 1 (toInteger ((*) 4 ((+) 1 n))))) numerat);
throwpi::B;
throwpi = (scalar_mult (Pi_mult_B (negate 2)) scals_denon)
}
 in (scalar_mult W_mult_B throwpi));
data Scalar_to_mult_B = I_mult_B | Pi_mult_B(Int) | W_mult_B | Rational_mult_B(Rational) deriving (Show);
scalar_mult::Scalar_to_mult_B->B->B;
scalar_mult s b = (case s of {
  (I_mult_B)->(map i_mult_poly_term b);
  (Pi_mult_B(pow))->(map (add_pi_pow pow) b);
  (Rational_mult_B(x))->(map (rational_mult_poly_term x) b);
  (W_mult_B)->((:) Zero_poly_term b)
});
rational_mult_poly_term::Rational->Poly_term->Poly_term;
rational_mult_poly_term x p = (manip_poly_term (rational_mult_pure_rational x) p);
rational_mult_pure_rational::Rational->Pure_rational->Pure_rational;
rational_mult_pure_rational s x = (case x of {
  (Pure_rational(a) (i))->(Pure_rational ((*) s a) i)
});
negate_pure_rational::Pure_rational->Pure_rational;
negate_pure_rational x = (rational_mult_pure_rational ((%) (negate 1) 1) x);
manip_poly_term::(Pure_rational->Pure_rational)->Poly_term->Poly_term;
manip_poly_term f p = (case p of {
  (Poly_term(frac) (pow))->(Poly_term (f frac) pow);
  (_)->p
});
add_pi_pow::Int->Poly_term->Poly_term;
add_pi_pow n p = (case p of {
  (Poly_term(frac) (Power_of_pi(pow)))->(Poly_term frac (Power_of_pi ((+) pow n)));
  (Zero_poly_term)->p
});
i_mult_poly_term::Poly_term->Poly_term;
i_mult_poly_term p = (manip_poly_term i_mult p);
i_mult::Pure_rational->Pure_rational;
i_mult x = (case x of {
  (Pure_rational(a) (Real))->(Pure_rational a Imaginary);
  (Pure_rational(a) (Imaginary))->(Pure_rational (negate a) Real)
});
add_b::B->B->B;
add_b x y = (subtract_b x (negate_Bs y));
subtract_b::B->B->B;
subtract_b x y = (zipWith_lengthen subtract_poly_term x y);
zipWith_lengthen::(Poly_term->Poly_term->Poly_term)->B->B->B;
zipWith_lengthen f x y = (case (x,y) of {
  (((:)(x1) (xr)),((:)(y1) (yr)))->((:) (f x1 y1) (zipWith_lengthen f xr yr));
  (([]),((:)(y1) (yr)))->((:) (f Zero_poly_term y1) (zipWith_lengthen f [] yr));
  (((:)(x1) (xr)),([]))->((:) (f x1 Zero_poly_term) (zipWith_lengthen f xr []));
  (([]),([]))->[]
});
subtract_poly_term::Poly_term->Poly_term->Poly_term;
subtract_poly_term x y = (case (x,y) of {
  ((Poly_term(frac1) (pow1)),(Poly_term(frac2) (pow2)))->(case ((==) pow1 pow2) of {
  (True)->(Poly_term (subtract_pure_rational frac1 frac2) pow1);
  (_)->(error (concat ["tried to subtract-poly-term ",(show x)," and ",(show y)]))
});
  ((_),(Zero_poly_term))->x;
  ((Zero_poly_term),(_))->(negate_poly_term y)
});
negate_poly_term::Poly_term->Poly_term;
negate_poly_term y = (rational_mult_poly_term ((%) (negate 1) 1) y);
negate_Bs::B->B;
negate_Bs y = (map negate_poly_term y);
subtract_pure_rational::Pure_rational->Pure_rational->Pure_rational;
subtract_pure_rational x y = (case (x,y) of {
  ((Pure_rational(a) (Real)),(Pure_rational(b) (Real)))->(Pure_rational ((-) a b) (Real ));
  ((Pure_rational(a) (Imaginary)),(Pure_rational(b) (Imaginary)))->(Pure_rational ((-) a b) (Imaginary ))
});
add_pure_rational::Pure_rational->Pure_rational->Pure_rational;
add_pure_rational x y = (subtract_pure_rational x (negate_pure_rational y));
all_bs::B_list;
all_bs = (drop 2 long_all_bs);
all_cs::B_list;
all_cs = (map c_calc (enumFrom 0));
pretty_b_matlab::(Poly_term,W_power)->String;
pretty_b_matlab b = (case b of {
  ((Poly_term(rat) (Power_of_pi(x))),(W_power(w)))->(concat [(pretty_rat_matlab rat x),"*omega^",(show w)])
});
pretty_b::(Poly_term,W_power)->String;
pretty_b b = (case b of {
  ((Zero_poly_term),(_))->"0";
  ((Poly_term(rat) (Power_of_pi(x))),(W_power(w)))->(concat [(pretty_rat rat x),"\\omega^",(show w)])
});
pretty_c_phi::(Poly_term,Phi_power)->String;
pretty_c_phi b = (case b of {
  ((Zero_poly_term),(_))->"0";
  ((Poly_term(rat) (Power_of_pi(x))),(Phi_power(w)))->(concat [(pretty_rat rat x),"\\Psi^{(",(show w),")}(p)"])
});
pretty_c_phi_matlab::(Poly_term,Phi_power)->String;
pretty_c_phi_matlab b = (case b of {
  ((Zero_poly_term),(_))->"0";
  ((Poly_term(rat) (Power_of_pi(x))),(Phi_power(w)))->(concat [(pretty_rat_matlab rat x),"*sym('P",(show w),"')"])
});
pretty_big_C::BC_term->String;
pretty_big_C b = (case b of {
  (BC_term(rat) (Power_of_pi(x)) (_) (Phi_power(w)))->(concat [(pretty_rat rat x),"\\Psi^{(",(show w),")}(p)"])
});
pretty_rat::Pure_rational->Int->String;
pretty_rat x pi = (case x of {
  (Pure_rational(x) (i))->(concat ["\\frac{",(show (numerator x)),"}{",(show (denominator x)),(maybe_pi_denom pi),"}",(maybe_pi_numerator pi),(maybe_imag i)])
});
pretty_rat_matlab::Pure_rational->Int->String;
pretty_rat_matlab x pi = (case x of {
  (Pure_rational(x) (i))->(concat ["sym('",(show (numerator x)),"')/sym('",(show (denominator x)),"')*sym('pi')^(",(show pi),")*(1",(maybe_imag i),")"])
});
fold_plus::String->String->String;
fold_plus a b = (concat [a," + ",b]);
maybe_imag::Imaginess->String;
maybe_imag i = (case i of {
  (Real)->"";
  (Imaginary)->"i"
});
maybe_pi_numerator::Int->String;
maybe_pi_numerator pi = (case (compare pi 0) of {
  (EQ)->"";
  (GT)->(concat [" \\pi^{",(show pi),"}"]);
  (_)->""
});
maybe_pi_denom::Int->String;
maybe_pi_denom pi = (case (compare pi 0) of {
  (EQ)->"";
  (LT)->(concat [" \\pi^{",(show (negate pi)),"}"]);
  (_)->""
});
not_Zero_poly::(Poly_term,a)->Bool;
not_Zero_poly x = (case (fst x) of {
  (Zero_poly_term)->False;
  (_)->True
});
show_a_B_matlab::Int->String;
show_a_B_matlab n = (concat [(foldr1 fold_plus (map pretty_b_matlab (reverse (filter not_Zero_poly (annotate_power W_power ((!!) all_bs n)))))),"\n"]);
show_a_B::Int->String;
show_a_B n = (concat ["b_{",(show n),"} &=& ",(foldr1 fold_plus (map pretty_b (reverse (filter not_Zero_poly (annotate_power W_power ((!!) all_bs n)))))),"\\\\\n"]);
show_a_C_matlab::Int->String;
show_a_C_matlab n = (concat [(foldr1 fold_plus (map pretty_c_phi_matlab (reverse (filter not_Zero_poly (annotate_power Phi_power ((!!) all_cs n)))))),"\n"]);
show_a_C::Int->String;
show_a_C n = (concat ["c_{",(show n),"} &=& ",(foldr1 fold_plus (map pretty_c_phi (reverse (filter not_Zero_poly (annotate_power Phi_power ((!!) all_cs n)))))),"\\\\\n"]);
show_a_big_C::Int->String;
show_a_big_C n = (concat ["C_{",(show n),"} &=& ",(foldr1 fold_plus (map pretty_big_C (reverse (terms_of_C 60 60 n)))),"\\\\\n"]);
c_calc::Int->B;
c_calc n = (extract_id (do{
  x<-(Id (foldl1 add_b (map (inner_sum n) (enumFromTo 0 (floor ((%) (toInteger n) 2))))));
  x<-(Id (scalar_mult (Rational_mult_B ((%) ((!!) fact n) ((^) 2 n))) x));
  (return x);
}));
raw_polynomial_bc::Int->[](BC_term);
raw_polynomial_bc m = (filter not_zero_BC (concat (take m (zipWith mult_B_and_C all_bs all_cs))));
is_W_pow::Int->BC_term->Bool;
is_W_pow pow bc = (case bc of {
  (BC_term(_) (_) (W_power(w)) (_))->((==) pow w)
});
is_pi_pow::Int->BC_term->Bool;
is_pi_pow pow bc = (case bc of {
  (BC_term(_) (Power_of_pi(w)) (_) (_))->((==) pow w)
});
raw_C::Int->[](BC_term)->[](BC_term);
raw_C n terms = (filter (is_W_pow n) terms);
sum_pow_pi::[](BC_term)->Int->BC_term;
sum_pow_pi terms which_pi = (foldl sum_identical_BC_terms Zero_BC_term (filter (is_pi_pow which_pi) terms));
sum_identical_BC_terms::BC_term->BC_term->BC_term;
sum_identical_BC_terms a b = (case (a,b) of {
  ((Zero_BC_term),(_))->b;
  ((BC_term(rat1) (xx) (yy) (zz)),(BC_term(rat2) (xx2) (yy2) (zz2)))->(case ((==) (xx,yy,zz) (xx2,yy2,zz2)) of {
  (True)->(BC_term (add_pure_rational rat1 rat2) xx yy zz)
})
});
collected_C::Int->[](BC_term)->[](BC_term);
collected_C n terms = (map (sum_pow_pi (raw_C n terms)) (enumFromThen 0 (negate 1)));
fact::[](Integer);
fact = (let {
fact1::Integer->Integer->[](Integer);
fact1 v n = ((:) v (fact1 ((*) v n) ((+) 1 n)))
}
 in (fact1 1 1));
compose_pow::Int->(B->B)->B->B;
compose_pow n f b = ((!!) (iterate f b) n);
newtype Id(a) = Id(a) deriving (Show);
extract_id::Id(a)->a;
extract_id x = (case x of {
  (Id(y))->y
});
instance Monad (Id) where {
return v = (Id v);
(>>=) m f = (f (extract_id m))
}
;
i::((a->b)->(a->Id(b)));
i = ((.) Id);
inner_sum::Int->Int->B;
inner_sum n j = (let {
the_pow::Int;
the_pow = ((*) 1 j);
n_minus_2j::Int;
n_minus_2j = ((-) n ((*) 2 j))
}
 in (extract_id (do{
  x<-(Id (phi_pow n_minus_2j));
  x<-(Id (compose_pow the_pow (scalar_mult I_mult_B) x));
  x<-(Id (scalar_mult (Pi_mult_B the_pow) x));
  x<-(Id (scalar_mult (Rational_mult_B ((%) ((^) 2 the_pow) 1)) x));
  x<-(Id (scalar_mult (Rational_mult_B ((%) 1 ((!!) fact ((*) 1 j)))) x));
  x<-(Id (scalar_mult (Rational_mult_B ((%) 1 ((!!) fact n_minus_2j))) x));
  (return x);
})));
data BC_term = Zero_BC_term | BC_term(Pure_rational)(Power_of_pi)(W_power)(Phi_power) deriving (Show);
type BC = [](BC_term);
newtype W_power = W_power(Int) deriving (Show,Eq);
newtype Phi_power = Phi_power(Int) deriving (Show,Eq);
annotate_power::(Int->b)->[](a)->[]((a,b));
annotate_power f l = (zip l (map f (enumFrom 0)));
mult_B_and_C::B->B->BC;
mult_B_and_C b c = (do{
  y<-(annotate_power Phi_power c);
  x<-(annotate_power W_power b);
  (return (mult_W_Phi x y));
});
not_zero_BC::BC_term->Bool;
not_zero_BC x = (case x of {
  (Zero_BC_term)->False;
  (BC_term(rat) (_) (_) (_))->(not (is_rational_zero rat))
});
is_rational_zero::Pure_rational->Bool;
is_rational_zero x = ((==) 0 (get_number_of_pure x));
mult_W_Phi::(Poly_term,W_power)->(Poly_term,Phi_power)->BC_term;
mult_W_Phi w phi = (case ((fst w),(fst phi)) of {
  ((Zero_poly_term),(_))->Zero_BC_term;
  ((_),(Zero_poly_term))->Zero_BC_term;
  ((Poly_term(rat1) (Power_of_pi(pow1))),(Poly_term(rat2) (Power_of_pi(pow2))))->(BC_term (mult_pure_rational rat1 rat2) (Power_of_pi ((+) pow1 pow2)) (snd w) (snd phi))
});
mult_pure_rational::Pure_rational->Pure_rational->Pure_rational;
mult_pure_rational x y = (let {
real_y::Rational;
real_y = (get_number_of_pure y);
almost::Pure_rational;
almost = (rational_mult_pure_rational real_y x)
}
 in (case y of {
  (Pure_rational(_) (Real))->almost;
  (Pure_rational(_) (Imaginary))->(i_mult almost)
}));
phi_pow::Int->B;
phi_pow j = (case (compare j 0) of {
  (EQ)->one_B;
  (GT)->((:) Zero_poly_term (phi_pow ((-) j 1)))
});
terms_of_C::Int->Int->Int->[](BC_term);
terms_of_C cut0 cut n = (filter not_zero_BC (take cut (collected_C n (raw_polynomial_bc cut0))));
main::IO(());
main = (do{
  (putStrLn "syms omega");
  (putStr "b=");
  (beg_eqnarray );
  (putStr (concatMap show_a_B_matlab (enumFromTo 0 200)));
  (end_eqnarray );
  (putStr "c=");
  (beg_eqnarray );
  (putStr (concatMap show_a_C_matlab (enumFromTo 0 200)));
  (end_eqnarray );
});
beg_eqnarray::IO(());
beg_eqnarray = (putStr "[");
end_eqnarray::IO(());
end_eqnarray = (putStrLn "];");
show_list::(Show a)=>[](a)->String;
show_list l = (let {
a_show elem = ((++) (show elem) "\n")
}
 in (foldr (++) "" (map a_show l)))
}