order=300;
precision=25;

(*makeRoots[α_,β_]:=Table[Sort[Re[x/.{ToRules[NRoots[0==JacobiP[n,α,β,2x-1],x,PrecisionGoal->precision]]}]],{n,1,order}];*)
makeRoots[α_,β_]:=Table[Sort[Re[N[x/.{ToRules[Roots[0==JacobiP[n,α,β,2x-1],x]]},precision]]],{n,1,order}];

funName[name_]:="shifted"<>name<>"Roots";
toHaskell[roots_] :=StringReplace[ToString[roots],{"{{"->"  [ [","{"->"[","},"->"]\n  ,","}}"->"]\n  ]\n"}];

(*specialize[name_,type_]:="{-# SPECIALIZE "<>funName[name]<>" :: [["<>type<>"]] #-}\n";*)
specialize[name_,type_]:="";

doc[name_,α_,β_]:="-- | roots of shifted Jacobi polynomials with alpha="<>ToString[α]<>", beta="<>ToString[β]<>"\n";

(*typeSignature[name_]:=funName[name]<>" :: Fractional a => [[a]]\n"*)
typeSignature[name_]:=funName[name]<>" :: [[Double]]\n"

writeFun[name_,α_,β_]:=doc[name,α,β]<>specialize[name,"Double"]<>specialize[name,"Float"]<>typeSignature[name]<>funName[name]<>" =\n"<>toHaskell[makeRoots[α,β]];

header = "{-# OPTIONS_GHC -Wall #-}\n\nmodule JacobiRoots ( "<>funName["Legendre"]<>", "<>funName["Radau"]<>" ) where\n";

t0=AbsoluteTime[];
radauFun=writeFun["Radau",1,0];
t1=AbsoluteTime[];
"Radau roots computed in "<>ToString[t1-t0]<>" seconds"

t0=AbsoluteTime[];
legendreFun=writeFun["Legendre",0,0];
t1=AbsoluteTime[];
"Legendre roots computed in "<>ToString[t1-t0]<>" seconds"

module=header<>"\n"<>legendreFun<>"\n"<>radauFun;

Export[NotebookDirectory[]<>"JacobiRoots.hs",module,"String"];