Fast Decompositions 
by Gegenbauer Family Orthogonal Polynomials

Package of programs  
                    NGegenbauerCListCoefficients, NGegenbauerCSeries, NChebyshevTListCoefficients,
                NChebyshevTSeries, NChebyshevUListCoefficients, NChebyshevUSeries,    
                NLegendrePListCoefficients,NLegendrePSeries

(*:Summary: This package calculates numerical  values of coefficients in the series expansions by Gegenbauer, LegendreP, ChebyshevT, ChebyshevU orthogonal polynomials and constructs the corresponding expansions *)
  
  
  
  
  
   
(*:Sources:
 1.  A . Nersessian, Periodic Decomposition of Orthogonal Systems, (In Russian), Preprint N21 - Ar00, Deposited in ArmNIINTI 09.03  .00, pp  .1-28. 
 2. A . Nersessian .   Quasiexponential Function and Fourier Quasitransforms (in Russian). Reports of NAS RA (Dokladi Nacionalnoj Akademii Nauk Armenii), vol. 101, N1, 2001, pp 5-11.
*)
(* Limitations: Programs  NChebyshevTListCoefficients[],  NChebyshevUListCoefficients[],NChebyshevUSeries[] and  NChebyshevTSeries[] actually have no any limitations;  
Programs NGegenbauerCListCoefficients[f, x, λ, n] and  NGegenbauerCSeries[
                  f, x, λ, n] are best when  - 
              0.3 ≤ λ ≤ 9, n ≤ 150 and PrecisionGoal ≤ 
            11; 
 
          Programs NLegendrePSeries[f, x, 
                n] and NLegendrePListCoefficients[f, x, 
                n] are best when n < 150 and PrecisionGoal ≤ 
            11  *)

BeginPackage["NumericalMath`Coefficients`"];
Off[General::"spell1"];
Off[
        General::"spell"];
Off[
        FindRoot::"frsec"];
Options[
          NGegenbauerCListCoefficients] = 
Options[
            NGegenbauerCSeries] = 
Options[
              NChebyshevTListCoefficients] = Options[
                NChebyshevUListCoefficients] = Options[
                  NLegendrePSeries] = 
Options[
                    NLegendrePListCoefficients] = 
Options[
                      NLegendrePSeries] = 
Options[
                        NChebyshevTSeries] = 
Options[
                          NChebyshevUSeries] = 
Options[
                          NIntegrate];
\(NGegenbauerCListCoefficients::usage = "NGegenbauerCListCoefficients[expr,t,λ,n] gives a list that contains the numerical  values of  Gegenbauer coefficients from 0 to n-1 in the Gegenbauer  series expansion on the interval (-1,1).";\)

NGegenbauerCSeries::usage = "NGegenbauerCSeries[expr,t,\
λ, n] gives the nth order Gegenbauer series
approximation to the  function of t  that is equal to expr for -1 <= t <= \
1.";

\(NLegendrePListCoefficients::usage = "NLegendrePListCoefficients[expr,t,n] gives a list that contains the numerical values of the LegendreP coefficients from 0 to n-1 in the LegendreP series expansion on the interval (-1,1).";\)

NLegendrePSeries::usage = "NLegendrePSeries[expr,t,n] gives \
the nth order LegendreP series
approximation to the  function of t  that is equal to expr for -1 <= t <= \
1.";

\(NChebyshevTListCoefficients::usage = "NChebyshevTListCoefficients[expr,t,n] gives a list that contains the numerical values of the ChebyshevT coefficients from 0 to n-1 in the ChebyshevT series expansion on the interval (-1,1).";\)

NChebyshevTSeries::usage = "NChebyshevTSeries[expr,t,n] gives \
the nth order ChebyshevT series
approximation to the  function of t  that is equal to expr for -1 <= t <= \
1.";

\(NChebyshevUListCoefficients::usage = "NChebyshevUListCoefficients[expr,t,n] gives a list that contains the numerical  values of the ChebyshevU coefficients from 0 to n-1 in the ChebyshevU  series expansion on the interval (-1,1).";\)

NChebyshevUSeries::usage = "NChebyshevUSeries[expr,t,n] gives \
the nth order ChebyshevT series
approximation to the  function of t  that is equal to expr for -1 <= t <= \
1.";

NGegenbauerCSeries::par = "Parameter λ in NGegenbauerCSeries[f,x,\
λ,n] must be greater than -1/2.";
NGegenbauerCSeries::par1 = "n≥70 in \
NGegenbauerCSeries[f,x,λ,n] is not recommended due to roundoff \
errors.";

NGegenbauerCListCoefficients::par = "Parameter λ in \
NGegenbauerCListCoefficients[f,x,λ,n] must be greater than \
-1/2.";
\(NGegenbauerCListCoefficients::par1 = "λ>15 in NGegenbauerCListCoefficients[f,x,λ,n] is not recommended. Try to increase WorkingPrecision or change other options";\)
\(NGegenbauerCListCoefficients::par2 = "λ<-0.3 in NGegenbauerCListCoefficients[f,x,λ,n] is not recommended. Try to increase WorkingPrecision or change other options";\)
\(NGegenbauerCListCoefficients::par3 = "PrecisionGoal greater than 11 is not recommended. Try to increase WorkingPrecision or change other options.";\)
NGegenbauerCListCoefficients::par4 = "n≥100 in \
NGegenbauerCListCoefficients[f,x,λ,n] for 9<λ<= 15 is not \
recommended. Try to increase WorkingPrecision or change other options.";
NGegenbauerCListCoefficients::par5 = "n≥100 in \
NGegenbauerCListCoefficients[f,x,λ,n] for -0.3≤λ<0 \
is not recommended. Try to increase WorkingPrecision or change other options.\
";

NGegenbauerCListCoefficients::par6 = "n>200 in \
NGegenbauerCListCoefficients[f,x,λ,n] for 0<λ<1  with \
PrecisionGoal≤8 is not recommended. Try to increase \
WorkingPrecision or change other options.";
NGegenbauerCListCoefficients::par7 = "n>120 in \
NGegenbauerCListCoefficients[f,x,λ,n] for 0<λ<1 with \
8<PrecisionGoal≤11 is not recommended. Try to increase \
WorkingPrecision or change other options.";
NGegenbauerCListCoefficients::par8 = "n>200 in \
NGegenbauerCListCoefficients[f,x,λ,n] for 1<λ<9 with \
PrecisionGoal≤8 is not recommended. Try to increase \
WorkingPrecision or change other options.";
NGegenbauerCListCoefficients::par9 = "n>120 in \
NGegenbauerCListCoefficients[f,x,λ,n] for 1<λ<9 with \
8<PrecisionGoal≤11 is not recommended. Try to increase \
WorkingPrecision or change other options.";


Begin["`Private`"];
mm[v_, prec_] := 
        Block[{ls1, ls2, ls, 
            int}, 
ls1 = {-0.5, 0, 3, 5, 10, 15, 20, 
              25, 30}; 
ls2 = {1, 1, 7, 12, 22, 31, 39, 45, 
              50}; 
ls = 
            Transpose[{ls1, ls2}]; 
int = 
            Interpolation[ls, 
              InterpolationOrder -> 2]; 
Min[30, 
            Ceiling[int[
                  v] (1 + 
                    (20/((v + 1.5) (prec + 2))))]]
];
NGegenbauerCListCoefficientsDirect[f_, x_, v_, nn_, opt___] := 
        Module[{d, r, l1, l2, l3}, 
If[
            v ≤  -1/2, 
            Message[NGegenbauerCListCoefficients::par]; 
            Break[]]; 
Table[(Sqrt[(((2^(2  v - 
                              1)) (n + v) n! 
                        Gamma[v]^2)/( π Gamma[
                          2  v + n]))]) 
              NIntegrate[
                Evaluate[
                  f  GegenbauerC[n, v, 
                      x] (((1 - (x^2)))^(v - 1/2))], {x, -1, 1}, 
                Evaluate[opt]], {n, 0, 
              nn - 1}]
];
NGegenbauerCListCoefficients[f_, x_, v_, nn_, opt___] := 
        Module[{d, r, l1, l2, 
            l3}, 
{prec} = {PrecisionGoal} /. {opt} /. 
              Options[NGegenbauerCListCoefficients]; 
If[
            prec === Automatic, 
            prec = $MachinePrecision - 
                11]; 
If[(9 < v ≤  
                  15) && (nn ≥ 100), 
            Message[NGegenbauerCListCoefficients::par4]]; 
If[prec > 11, 
            Message[NGegenbauerCListCoefficients::par3]]; 
If[v ≤  -1/2, 
            Message[NGegenbauerCListCoefficients::par]; 
            Break[]]; 
If[v > 15, 
            Message[NGegenbauerCListCoefficients::par1]]; 
If[v < -0.3, 
            Message[NGegenbauerCListCoefficients::par2]]; 
If[(-0.3 ≤  v < 
                  0) && (nn ≥ 100), 
            Message[NGegenbauerCListCoefficients::par5]]; 
If[
            0 < v < 1, 
     Which[(prec ≤  8) && (nn > 200), 
              Message[NGegenbauerCListCoefficients::par6], 
                 (prec ≤ 
                    11) && (nn > 120), 
              Message[NGegenbauerCListCoefficients::par7]]
  ]; 
If[
            1 < v ≤  
              9, 
     Which[(prec ≤  
                    8) && (nn > 200), 
              Message[NGegenbauerCListCoefficients::par8], 
                 (prec ≤ 
                    11) && (nn > 120), 
              Message[NGegenbauerCListCoefficients::par9]]
  ]; 

Which[
            Abs[v] ≤  (10^(-16)), 
            NChebyshevTListCoefficients[f, x, nn, opt], 
            Abs[v - 1] ≤ (10^(-16)), 
            NChebyshevUListCoefficients[f, x, nn, 
              opt], ((Abs[v] > (10^(-16))) || (Abs[v - 1] > 
                    (10^(-16)))), 
If[(nn ≤  30) || (v > 25), 
              NGegenbauerCListCoefficientsDirect[f, x, v, nn, 
                opt], 
If[prec ≤  5, 
                prec = 5]; 
p1 = 
                Max[ .5,  (prec - 3)]; 
r = 
                Which[v < 0.4, 
                  (10^(-(1/(4  (4^(-(nn/200))) nn)))), 
                  0.4 ≤  v < 1, 
                  (10^(-((Sqrt[p1])/(4  (4^(-(nn/200))) nn)))), 
                  1 ≤  v < 7, 
                  (nn^(-((( .4 - ((Sqrt[p1 + 1])/4.35)  .29) (p1) (Sqrt[v + 1]))/((nn - 1.7  Log[2, nn]) Log[4  (v^2) + 2 + ((2  nn)/(nn + 60))])))), True, 
                  (10^(-((Sqrt[p1])/(4  (4^(-(nn/200))) nn))))]; 
 
              l1 = Which[v < 1/2, 
                  Table[NIntegrate[
                      Evaluate[(f /. 
                              x -> 
                                Cos[x]) (((1 - 
                                  2  Cos[
                                      x] r (E^(((2  I π  k)/nn))) + (r^2) (E^(((4  I π  k)/nn)))))^(-v)) 
                          Sin[x]^(2 v)], {x, 0, π}, 
                      Evaluate[{opt}]], {k, 0, nn - 1}], 
                  1/2 ≤  v < 7, 
                  Table[NIntegrate[
                      Evaluate[
                        f (((1 - 
                                  2  r x (E^(((2  I π  k)/nn))) + (r^2) (E^(((4  I π  k)/nn)))))^(-v)) (((1 - (x^2)))^(v - 1/2))], {x, -(Sqrt[1 - (10^(-(45/(v + 2))))]), -(Sqrt[1 - (10^(-(25/(5 (v + 2)))))]), (Sqrt[1 - (10^(-(25/(5 (v + 2)))))]), (Sqrt[1 - (10^(-(45/(v + 2))))])}, Evaluate[{opt}]], {k, 0, nn - 1}], True, 
                  Table[NIntegrate[
                      Evaluate[(f /. 
                              x -> 
                                Cos[x]) (((1 - 
                                  2  r Cos[
                                      x] (E^(((2  I π  k)/nn))) + (r^2) (E^(((4  I π  k)/nn)))))^(-v)) 
                          Sin[x]^(2 v)], {x, -ArcCos[(Sqrt[1 - (10^(-(45/(v + 2))))])], -ArcCos[(Sqrt[1 - (10^(-(25/(5 (v + 2)))))])], ArcCos[(Sqrt[1 - (10^(-(25/(5 (v + 2)))))])], 
                        ArcCos[(Sqrt[1 - (10^(-(45/(v + 2))))])]}, 
                      Evaluate[{opt}]], {k, 0, 
                      nn - 1}]]; 

l2 =  
                InverseFourier[l1]; 
l3 = 
                Prepend[
                  Table[((r^(-n))/(Sqrt[nn])) (Sqrt[(((2^(2  v - 1)) (n + v) n! Gamma[v]^2)/( π Gamma[2  v + n]))]), {n, 1, 
                      nn - 1}], ((Sqrt[((v Gamma[v])/(nn Gamma[(1/2) + v]))])/(π^(1/4)))]; 
aux = 
                l3   l2; 
If[
                v ≤  
                  2,  
           Return[
                  aux], 
        Join[ 
                  Take[NGegenbauerCListCoefficientsDirect[f, x, v, nn, opt], 
                    mm[v, prec]], 
                  Drop[aux, 
                    mm[v, prec]]]
  ]]]
];
NGegenbauerCSeries[f_, 
          x_, v_, nn_, opt___] := 
        Module[{cof, 
            geglist}, 
{prec} = {PrecisionGoal} /. {opt} /. Options[
                NGegenbauerCListCoefficients]; 

If[v ≤  -1/2, 
            Message[NGegenbauerCSeries::par]; Break[]]; 
If[
            nn ≥  70, 
            Message[NGegenbauerCSeries::par1]]; 
If[
            prec === Automatic, 
            prec = $MachinePrecision - 11]; 
Which[
            Abs[v] ≤  (10^(-16)), 
            NChebyshevTSeries[f, x, nn, opt], 
            Abs[v - 1] ≤  (10^(-16)), 
            NChebyshevUSeries[f, x, nn, 
              opt], 
True, 
cof = 
              NGegenbauerCListCoefficients[f, x, v, nn, 
                opt]; 
geglist = 
              Table[(Sqrt[(((2^(2  v - 1)) (n + v) n! 
                            Gamma[v]^2)/( π Gamma[
                              2  v + n]))]) GegenbauerC[n, v, x], {n, 
                  0,  nn - 1}]; 
cof . 
              geglist]
];
NLegendrePListCoefficients[f_, x_, nn_, opt___] := 
        Module[{}, 
NGegenbauerCListCoefficients[f, x, 1/2, 
            nn, opt]
];
NLegendrePSeries[f_, x_, 
          nn_, opt___] := 
        Module[{}, 
NGegenbauerCSeries[f, x, 1/2, nn, 
            opt]
];
cofT[n_, f_, x_, opt___] := 
        2  NIntegrate[
            Evaluate[(f /. x -> Cos[π x]) Cos[
                  n π x] ], {x, 0, 1}, 
            Evaluate[{opt}]];
cofU[n_, f_, x_, opt___] := 
        2  NIntegrate[
            Evaluate[(f /. 
                    x -> 
                      Cos[π x]) Sin[(n  + 
                        1) π x]  Sin[π x]], {x, 0, 1}, 
            Evaluate[{opt}]];
NChebyshevTListCoefficients[f_, x_, 
          nn_, opt___] := 
Block[{l1, l2, 
            l3}, 
norm = 
            Flatten[{π, 
                Table[π/2, {i, 1, nn - 1}]}]; 
If[
            nn ≤ 30, 
            l1 = (Sqrt[norm ]) 
                Table[cofT[n, f, x, opt], {n, 0, nn - 1}]; 
            Return[ReplacePart[l1, l1[[1]]/2, 
                1]], 
l1 = 
              NFourierListCoefficients[(f /. x -> Cos[π x]), x, 
                nn - 1, FourierParameters -> {1, 1/2}, 
                Evaluate[opt]]; 
l2 = 
              Reverse[l1]; 
l3 = 
              Take[(l1 + l2)/2, -nn]]; 
ReplacePart[
              l3, l3[[1]]/2, 
              1] (Sqrt[norm])
];
NChebyshevUListCoefficients[f_, x_, nn_, 
          opt___] := 
Block[{l1, l2}, 
If[
            nn ≤ 30, 
            l1 = Table[(Sqrt[(π/2)]) cofU[n, f, x, opt], {n, 0, 
                  nn - 1}]; 
Return[
              l1], 
l1 = 
              NFourierListCoefficients[
                Evaluate[(f /. x ->  Cos[π x]) Sin[π x]], 
                x, nn, FourierParameters -> {1, 1/2}, 
                Evaluate[opt]]; 
l2 = 
              Reverse[l1]; 
Take[-I (l1 - l2)/
                    2, -nn] 
              Table[(Sqrt[(π/2)]), {n, 0, 
                  nn - 1}]]
];
NChebyshevUSeries[f_, x_, nn_, opt___] := 
        Module[{cof, l1}, 
cof = 
            NChebyshevUListCoefficients[f, x, nn, 
              opt]; 
l1 = 
            Table[(Sqrt[(2/π)]) ChebyshevU[n, x], {n, 0, 
                nn - 1}]; 
l1 . 
            cof
];
NChebyshevTSeries[f_, x_, nn_, opt___] := 
        Module[{cof, l1}, 
norm = 
            Flatten[{π, 
                Table[π/2, {i, 1, nn - 1}]}]; 
cof = 
            NChebyshevTListCoefficients[f, x, nn, 
              opt]; 
l1 = (1/(Sqrt[norm])) 
              Table[Cos[n ArcCos[x]], {n, 0, 
                  nn - 1}]; 
l1 . 
            cof
];
\(per[l_List, 
          p_Integer] := 
	
        Block[{nn, lb, le}, 
		nn = Length[l]; 
		lb = 
            Take[l, {1, p}]; 
		le = 
            Take[l, {nn - p + 1, nn}]; 
		Flatten[{le, l, 
              lb}]
];\)
nll[ll_List, m_] := 
        If[Length[ll] >= m, 
          ll, 
Flatten[{bB[Length[ll], m], ll, 
              eE[Length[ll], m]}]];
\(flst[l_List] := 
        Block[{fy, nn = Length[l]}, 
	fy = RotateLeft[l, Floor[nn/2]]; 

          RotateRight[Fourier[fy]/\@nn, Floor[nn/2]]
	];\)
\(invflst[
          l_List] := 
        Block[{fy, nn = Length[l]}, 
	fy = RotateLeft[l, Floor[nn/2]]; 

          RotateRight[InverseFourier[fy]/\@nn, 
            Floor[nn/2]]
	];\)
roG[x_, 
          al_] = (((Sqrt[π]) (E^(-(((π^2) (x^2))/(4 al)))))/(Sqrt[al]));
<< Calculus`FourierTransform`;
<< 
        Utilities`FilterOptions`;
NFourierListCoefficients[f_, x_, nm_, opt___] := 
        Module[{pg, a, b, m, it, 
            int}, 
{prec, {a, 
                b}} = {PrecisionGoal, FourierParameters} /. {opt} /. 
              Options[NGegenbauerCListCoefficients]; 
If[
            prec === Automatic, 
            prec = $MachinePrecision - 11]; 
If[
            nm ≤  30, 
            Return[Table[
                NFourierCoefficient[f, x, n, opt], {n, -nm, 
                  nm}]],     
it = {{30, 
                  3.7 (1 + (prec/100))}, {100, 
                  3.5 (1 + (prec/100))}, {250, 
                  3 (1 + (prec/100))}, {500, 
                  2.6 (1 + (prec/100))}, {750, 
                  2.4` (1 + (prec/100))}, {1000, 
                  2.3 (1 + (prec/100))}, {1500, 
                  2.1 (1 + (prec/100))}, {2000, 
                  1.93 (1 + (prec/100))}, {2500, 
                  1.82 (1 + (prec/100))}, {3000, 
                  1.75 (1 + (prec/100))}, {4000, 
                  1.67 (1 + (prec/100))}, {5000, 
                  1.65 (1 + (prec/100))}, {10000, 
                  1.6 (1 + (prec/100))}}; 
int = 
              Interpolation[it, 
                InterpolationOrder -> 1]; 
m = 
              If[nm < 10000, Rationalize[int[nm]], 
                16/10 (1 + (prec/100))]; 
If[
              b > 0, 
      Return[Abs[
                        b]^(-((a + 1)/2)) 
                  Reverse[
                    NFourierListCoefficientsH[(f /. 
                          x -> x/(2  Abs[b])), x, pg, {m, nm}, 
                      opt]]]
 ,  
              Return[Abs[b]^(-((a + 1)/2)) 
                  NFourierListCoefficientsH[(f /. 
                        x -> x/(2  Abs[b])), x, pg, {m, nm}, 
                    opt]]
   ]
]
];
\(NFourierListCoefficientsH[f_, x_, 
          pc_, {m_, nm_}, opt___] := 
	
        Block[{cof, cof1, ifo, lo, fco0, prec = pc, 
            nn = 2  Ceiling[m\ nm] + 1, α0, α, d0, p, h, tb1, 
            tb2, h1, t1, t2, 
            optI}, 
{prec} = \({PrecisionGoal} /. {opt}\) /. 
              Options[NGegenbauerCListCoefficients]; 
If[
            prec === Automatic, 
            prec = $MachinePrecision - 11]; 
optI = 
            FilterOptions[NIntegrate, opt]; 
d0 = 
            Abs[\(1\/\(nn\ \ π\)\) \((\((f /. 
                        x -> 1.  -  .3\/nm)\) - \((f /. 
                        x -> \(-1. \) +  .3\/nm)\))\)]; 
α0 = \(1\/\(m - 
                    1\)\) \((m \((Log[Max[2, 5\/2 - Round[prec]\/10]] + \ 
                      Log[10]\ prec + 
                      Log[d0 + 1\/100])\))\); 
α = 
            al /. FindRoot[
                al == 
                  SetPrecision[\(m\/\(m - 1\)\) \((Log[
                            Max[2, 5\/2 - Round[prec]\/10]] + \ 
                          Log[10]\ prec + Log[d0 + 1\/100] + Log[m/al]\/2)\), 
                    100], {al,  .85  α0,  .86  α0}, 
                MaxIterations -> 3, 
                WorkingPrecision -> 100]; 
h = 
            hh /. FindRoot[
                SetPrecision[\(E\^\(\(α\/4\) \((1\/\(2  m\))\)\^2\)\) \((\(-\(\((\@\(π\/2\)\ \((\(-1\) + Erf[\(\@2\ hh\ π\)\/\@α])\))\)/\((nn\ \@α)\)\)\))\)\^\(1/8\), 
                    100] == 
                  SetPrecision[10\^\(\(-prec\)/4\), 100], {hh,  .7  
                    prec,  .73  prec}, MaxIterations -> 3, 
                WorkingPrecision -> 100]; 
h = 
            If[prec < 3, prec\/\(m - 1. \), 
              Chop[h, 10\^\(-14\)]]; 
h1 = 
            Ceiling[h]; 
h = \(11\/10\) 
              h; 
pp = 
            Floor[\ nn/2]; 
p = 
            pp + h1; 
cof = 
            Table[NIntegrate[
                Evaluate[roG[nn\ x - 2  s, α] f], {x, 
                  Max[\(2  s - 2  h\)\/nn, \(-1\)], 
                  Min[\(2  s + 2  h\)\/nn, 1]}, 
                AccuracyGoal -> prec + 1, 
Evaluate[
                  optI]], {\ s, \(-p\), 
                nn - pp - 1 + h1}]; 
cof1 = 
            Take[cof, {h1 + 1, nn + h1}]; 
ifo = 
            nn\ invflst[cof1]; 
lo = 
            Table[\(1\/2\) E\^\(α \((\ n/nn)\)\^2\), {n, \(-pp\), 
                nn - pp - 1}]; 
fco0 = 
            lo\ ifo; 
tb2 = \(1\/2\) 
              Table[\ Evaluate[\(E\^\(α \((n/nn)\)\^2\)\) \((∑\+\(s = nn - pp\)\%\(nn - pp - 1 + h1\)\(E\^\(\(\(-2\) 
                                      I\ π\ s\ n\)\/nn\)\) 
                            cof[\([s + p + 
                                  1]\)] + ∑\+\(s = \(-p\)\)\%\(\(-pp\) - 1\)\(E\^\(\(\(-2\) I\ π\ s\ n\)\/nn\)\) 
                            cof[\([s + p + 1]\)])\)], {n, \(-pp\), 
                  nn - pp - 1}]; 
Take[
            fco0 + tb2, {Ceiling[m\ \ nm] + 1 - nm, 
              Ceiling[m\ nm] + 1 + 
                nm}]
];\)
End[];
EndPackage[];

Introduction

The efficient expansions by classical orthogonal polynomials is a well known and hard solving problem for applied mathematics.  At first sight there is no problem as  MATHEMATICA has a package for the fast generation of all known classical orthogonal polynomials and has a powerful package for numerical integration. However  when a user makes an expansion of a smooth function the answers come very slowly and with significant errors. For example, it happens when it is necessary to have more  than 20 terms in Gegenbauer expansion in a one-dimensional case  (see [1]).   The main reason of this phenomenon is high oscillation of all known functional orthogonal systems if the corresponding index is rather big. 
The second reason is connected with main coefficients of far polynomials. For example, the view of the polynomial  (GegenbauerC[75,2.1,x] by MATHEMATICA notation) is 

-4091.293483305945` x + 4.046834760787354`*^6 (x^3) - 
          1.1975393424121938`*^9 (x^5) + «48» + 
          6.753027227430814`*^26 (x^71) - 7.719213771385703`*^25 (x^73) + 
          4.233745354972603`*^24 (x^75)                                              

It is easy to see that if we use a 32-bit computer for calculation of a value of this polynomial for x close to one, then a big roundoff error may appears.
   This package is based on the papers [1,2] and contains an absolutely new approach for calculation of expansions by orthogonal  polynomials.
      
   [1] 	 A.Nersessian, Periodic Decomposition of Orthogonal Systems, (In Russian), Preprint N21-Ar00, Deposited in ArmNIINTI 09.03.00, pp.1-28. 
  [2]   A.Nersessian.  Quasiexponential Function and Fourier Quasitransforms (in Russian). Reports of NAS RA (Dokladi Nacionalnoj Akademii Nauk Armenii), vol. 101, N1, 2001, pp 5-11.

Implementation Notes

◼The package contains eight programs:
NGegenbauerCListCoefficients, NGegenbauerCSeries, NChebyshevTListCoefficients, NChebyshevTSeries, NChebyshevUListCoefficients, NChebyshevUSeries, NLegendrePSeries, NLegendrePListCoefficients.

◼ The  program NChebyshevTListCoefficients    gives the list of coefficients for ChebyshevT
orthonormal system
◼ The  program NChebyshevUListCoefficients    gives the list of coefficients for ChebyshevU
orthonormal system

◼ The  program NLegendrePListCoefficients    gives the list of coefficients for LegendreP
orthonormal system
◼ The  program NGegenbauerCListCoefficients   gives the list of coefficients for GegenbauerC orthonormal system
◼The programs which names are ended by the word Series gives corresponding trancated
orthogonal expansions.
◼ For values  =0,1/2 and 1 (see (8) above)  programs NGegenbauerCListCoefficients
and NGegenbauerCSeries automatically transfer to corresponding programs for ChebyshevT,
LegendreP and ChebyshevU polynomials.
◼ Programs  support all  NIntegrate \(\(options\)\(.\)\(
\)\)

◼Programs  support all  NIntegrate options.
◼When N≤30 exact evaluation of the corresponding integrals works.
◼The programs NChebyshevTListCoefficients[], NChebyshevTSeries[], NChebyshevUListCoefficients[], NChebyshevUSeries[] actually have no limitations.
◼Programs NGegenbauerCListCoefficients[] and  NGegenbauerCSeries[] are best when -0.3 ≤ λ ≤ 9, N ≤ 150 with PrecisionGoal ≤ 11. When  λ > 25 direct evaluation of the corresponding integrals works. For other values of λ,N and PrecisionGoal we recommend to increase WorkingPrecision or change other options.
◼ Programs NLegendrePSeries[] and NLegendrePListCoefficients[] are best when N < 150 and PrecisionGoal ≤ 11. For other values of  N and PrecisionGoal we recommend to increase WorkingPrecision or change other options.

Demonstration of the package

All examples were solved on the computer
 PENTIUM  4,  1.5 GHz,512 RAM

Gegenbauer

Here we  calculate  the Gegenbauer coefficients for the function  (E^(I a x )) (for this function we know exact coefficients)  and construct the corresponding  orthogonal decomposition. We compare two methods: the first (classic) is the direct evaluation of the corresponding integrals, the second (new) are  the programs NGegenbauerCListCoefficients[] and  NGegenbauerCSeries[]. We compare both the times and the precisions ((Subscript[L, 2])-error).

Now  λ = -0.1,  PrecisionGoal = 9 and  35≤N≤275.  The message show that the program  NGegenbauerCSeries[] doesn't give satisfactory precision for N > 70 due to roundoff errors.  Besides the program  NGegenbauerCListCoefficients[] is not recommended for N > 100 and  λ = -0.1 despite that for this example the program works best up to N=275.  For N=275  the precisions of the coefficients are uncomparable but they have spent actually the same time for calculations.

\(Do[\ 
Print["***** nn->", 
        nn]; 
Block[{λ = \(-0.1\), 
          p = 9}, 
Print["***** λ->", λ]; 
Print["***** PrecisionGoal->", 
          p]; 
a = 5.231; 
f[x_] = 
          E\^\(\(I\)\(\ \)\(a\)\(\ \)\(x\)\(\ \)\); 
cofexact = 
          Sign[λ]\ Table[\(\@\(\(2\^\(2  λ - 
                            1\)\) \((n + λ)\)\)\) \(2\^λ\/\@\(n!\)\) 
                Gamma[λ + 
                    1\/2] \(\(Gamma[λ] \@Gamma[2\ λ + n]\)\/Gamma[2\ λ]\) \(I\^n\) \(a\^\(-λ\)\) 
                BesselJ[λ + n, a], {n, 0, 
                nn - 1}]; 
{t1, s1} = 
          NGegenbauerCListCoefficients[f[x], x, λ, nn, 
              PrecisionGoal -> p] // Timing; 
{t2, 
            s2} = Table[\(\@\(\(\((n + λ)\)\ Gamma[λ]\)\/\(\@π\ Gamma[1\/2 + λ]\ GegenbauerC[n, λ, 1]\)\)\) 
                NIntegrate[
                  Evaluate[
                    f[Cos[x]] GegenbauerC[n, λ, Cos[x]] 
                      Sin[x]\^\(2\ λ\)], {x, 0, π}, 
                  PrecisionGoal -> p], {n, 0, nn - 1}] // 
            Timing; 
errornew = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s1]\^2]]; 
errorclassic = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s2]\^2]]; 
Print[\*""(Subscript[L, 2])-error->{classic,new,classic/new}->"", {errorclassic, errornew, 
            errorclassic/
              errornew}]; 
Print["Time->{classic,new,classic/new}->", {t2, t1, 
            t2/t1}]; 
\ 
        ort = Table[\(\@\(\(\((n + λ)\)\ Gamma[λ]\)\/\(\@π\ Gamma[1\/2 + λ]\ GegenbauerC[
                          n, λ, 1]\)\)\) 
              GegenbauerC[n, λ, x], {n, 0, 
              nn - 1}]; 
s3 = 
          NGegenbauerCSeries[f[x], x, λ, nn, 
            PrecisionGoal -> p]; 
s4 = 
          s2 . \ ort; 
g1 = 
          Plot[Evaluate[
              Abs[s3 - f[x]] \((1 - x\^2)\)\^\(λ/2 - 1/4\)], {x, 0, 
              1}, PlotRange -> All, DisplayFunction -> Identity, 
            PlotStyle -> {Hue[1]}, 
            PlotLabel -> "New"]; 
g2 = 
          Plot[Evaluate[
              Abs[s4 - f[x]] \((1 - x\^2)\)\^\(λ/2 - 1/4\)], {x, 0, 
              1}, PlotRange -> All, DisplayFunction -> Identity, 
            PlotLabel -> "Classic"]; 
Show[
          GraphicsArray[{g2, 
              g1}]]

], {nn, 35, 275, 
        20}]\)

Now  λ = 2.2,  PrecisionGoal = 11 and  35≤N≤275.  The situation is similar

\(Do[\ 
Print["***** nn->", 
        nn]; 
Block[{λ = 2.2, 
          p = 11}, 
Print["***** λ->", λ]; 
Print["***** PrecisionGoal->", 
          p]; 
a = 5.231; 
f[x_] = 
          E\^\(\(I\)\(\ \)\(a\)\(\ \)\(x\)\(\ \)\); 
cofexact = 
          Sign[λ]\ Table[\(\@\(\(2\^\(2  λ - 
                            1\)\) \((n + λ)\)\)\) \(2\^λ\/\@\(n!\)\) 
                Gamma[λ + 
                    1\/2] \(\(Gamma[λ] \@Gamma[2\ λ + n]\)\/Gamma[2\ λ]\) \(I\^n\) \(a\^\(-λ\)\) 
                BesselJ[λ + n, a], {n, 0, 
                nn - 1}]; 
{t1, s1} = 
          NGegenbauerCListCoefficients[f[x], x, λ, nn, 
              PrecisionGoal -> p] // Timing; 
{t2, 
            s2} = Table[\(\@\(\(\((n + λ)\)\ Gamma[λ]\)\/\(\@π\ Gamma[1\/2 + λ]\ GegenbauerC[n, λ, 1]\)\)\) 
                NIntegrate[
                  Evaluate[
                    f[Cos[x]] GegenbauerC[n, λ, Cos[x]] 
                      Sin[x]\^\(2\ λ\)], {x, 0, π}, 
                  PrecisionGoal -> p], {n, 0, nn - 1}] // 
            Timing; 
errornew = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s1]\^2]]; 
errorclassic = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s2]\^2]]; 
Print[\*""(Subscript[L, 2])-error->{classic,new,classic/new}->"", {errorclassic, errornew, 
            errorclassic/
              errornew}]; 
Print["Time->{classic,new,classic/new}->", {t2, t1, 
            t2/t1}]; 
\ 
        ort = Table[\(\@\(\(\((n + λ)\)\ Gamma[λ]\)\/\(\@π\ Gamma[1\/2 + λ]\ GegenbauerC[
                          n, λ, 1]\)\)\) 
              GegenbauerC[n, λ, x], {n, 0, 
              nn - 1}]; 
s3 = 
          NGegenbauerCSeries[f[x], x, λ, nn, 
            PrecisionGoal -> p]; 
s4 = 
          s2 . \ ort; 
g1 = 
          Plot[Evaluate[
              Abs[s3 - f[x]] \((1 - x\^2)\)\^\(λ/2 - 1/4\)], {x, 0, 
              1}, PlotRange -> All, DisplayFunction -> Identity, 
            PlotStyle -> {Hue[1]}, 
            PlotLabel -> "New"]; 
g2 = 
          Plot[Evaluate[
              Abs[s4 - f[x]] \((1 - x\^2)\)\^\(λ/2 - 1/4\)], {x, 0, 
              1}, PlotRange -> All, DisplayFunction -> Identity, 
            PlotLabel -> "Classic"]; 
Show[
          GraphicsArray[{g2, 
              g1}]]

], {nn, 35, 275, 
        20}]\)

Now  λ = 16,  PrecisionGoal = 8 and  35≤N≤95.  For λ > 15 we recommend to increase WorkingPrecision to acheive the required precision.  Here we set WorkingPrecision=17.   The time of calculations increases but the precision of the coefficients is satisfactory up to N=95.  For N=95 new method works 3 times slow  but by precision classic and new methods are uncomparable.

\(Do[\ 
Print["***** nn->", 
        nn]; 
Block[{λ = 16, 
          p = 8}, 
Print["***** λ->", λ]; 
Print["***** PrecisionGoal->", 
          p]; 
a = 5.231; 
f[x_] = 
          E\^\(\(I\)\(\ \)\(a\)\(\ \)\(x\)\(\ \)\); 
cofexact = 
          Sign[λ]\ Table[\(\@\(\(2\^\(2  λ - 
                            1\)\) \((n + λ)\)\)\) \(2\^λ\/\@\(n!\)\) 
                Gamma[λ + 
                    1\/2] \(\(Gamma[λ] \@Gamma[2\ λ + n]\)\/Gamma[2\ λ]\) \(I\^n\) \(a\^\(-λ\)\) 
                BesselJ[λ + n, a], {n, 0, 
                nn - 1}]; 
{t1, s1} = 
          NGegenbauerCListCoefficients[f[x], x, λ, nn, 
              PrecisionGoal -> p, WorkingPrecision -> 17] // 
            Timing; 
{t2, s2} = 
          Table[\(\@\(\(\((n + λ)\)\ Gamma[λ]\)\/\(\@π\ Gamma[1\/2 + λ]\ GegenbauerC[n, λ, 1]\)\)\) 
                NIntegrate[
                  Evaluate[
                    f[Cos[x]] GegenbauerC[n, λ, Cos[x]] 
                      Sin[x]\^\(2\ λ\)], {x, 0, π}, 
                  PrecisionGoal -> p, WorkingPrecision -> 17], {n, 
                0, nn - 1}] // Timing; 
errornew = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s1]\^2]]; 
errorclassic = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s2]\^2]]; 
Print[\*""(Subscript[L, 2])-error->{classic,new,classic/new}->"", {errorclassic, errornew, 
            errorclassic/
              errornew}]; 
Print["Time->{classic,new,classic/new}->", {t2, t1, 
            t2/t1}];
], {nn, 35, 95, 20}]\)

LegendreP

Here we  calculate  the LegendreP  coefficients for the function  (E^(I a x )) (for this function we know exact coefficients)  and construct the corresponding  orthogonal decomposition. We compare two methods: the first (classic) is the direct evaluation of the corresponding integrals, the second (new) are  the programs NLegendrePListCoefficients[] and  NLegendrePSeries[]. We compare both the times and the precisions ((Subscript[L, 2])-error). Here we set PrecisionGoal=11 and 35≤N≤275. The new method is much more precise and fast than classic.

\(Do[\ 
Print["***** nn->", 
        nn]; 
Block[{p = 11}, 
a = 
          5.231; 
f[x_] = 
          E\^\(\(I\)\(\ \)\(a\)\(\ \)\(x\)\(\ \)\); 
Print["***** PrecisionGoal->", 
          p]; 
cofexact = 
          Table[\(\@\((n + 1/2)\)\) \(\@2\/\@\(n!\)\) 
              Gamma[1/2] \(\@Gamma[1 + n]\) \(I\^n\) \(a\^\(\(-1\)/2\)\) 
              BesselJ[1/2 + n, a], {n, 0, nn - 1}]; 
{t1, 
            s1} = NLegendrePListCoefficients[f[x], x, nn, 
              PrecisionGoal -> p] // Timing; 
{t2, 
            s2} = Table[\(\@\(\(\((n + 1/2)\) \(n!\) 
                          Gamma[1/2]\^2\)\/\(\(\ \)\(π\ Gamma[
                            1 + n]\)\)\)\) 
                NIntegrate[Evaluate[f[x] LegendreP[n, x]], {x, \(-1\), 1}, 
                  PrecisionGoal -> p], {n, 0, nn - 1}] // 
            Timing; 
errornew = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s1]\^2]]; 
errorclassic = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s2]\^2]]; 
Print[\*""(Subscript[L, 2])-error->{classic,new,classic/new}->"", {errorclassic, errornew, 
            errorclassic/
              errornew}]; 
Print["Time->{classic,new,classic/new}->", {t2, t1, 
            t2/t1}]; 
ort = 
          Prepend[Table[\(\@\(\(\((n + 1/2)\) \(n!\) 
                          Gamma[1/2]\^2\)\/\(\(\ \)\(π\ Gamma[
                            1 + n]\)\)\)\) LegendreP[n, x], {n, 1, 
                nn - 1}], \@\(\(1/2\ Gamma[1/2]\)\/Gamma[1]\)\/π\^\(1/4\)]; 
s3 = 
          NLegendrePSeries[f[x], x, nn, 
            PrecisionGoal -> p]; 
s4 = 
          s2 . \ ort; 
g1 = 
          Plot[Evaluate[Abs[s3 - f[x]]], {x, \(-1\), 1}, 
            PlotRange -> All, DisplayFunction -> Identity, 
            PlotStyle -> {Hue[1]}, 
            PlotLabel -> "New"]; 
g2 = 
          Plot[Evaluate[Abs[s4 - f[x]]], {x, \(-1\), 1}, 
            PlotRange -> All, DisplayFunction -> Identity, 
            PlotLabel -> "Classic"]; 
Show[
          GraphicsArray[{g2, 
              g1}]]

], {nn, 35, 275, 
        20}]\)

ChebyshevU

Here we  calculate  the  ChebyshevU  coefficients for the function  (E^(I a x )) (for this function we know exact coefficients)  and construct the corresponding  orthogonal decomposition. We compare two methods: the first (classic) is the direct evaluation of the corresponding integrals, the second (new) are  the programs NChebyshevUListCoefficients[] and  NChebyshevUSeries[]. We compare both the times and the precisions ((Subscript[L, 2])-error). Here we set PrecisionGoal=11 and 35≤N≤175. The new method is much more precise and faster than classic.

\(Do[\ 
Print["***** nn->", 
        nn]; 
Block[{p = 
            11}, 
Print["***** PrecisionGoal->", 
          p]; 
a = 5.231; 
f[x_] = 
          E\^\(\(I\)\(\ \)\(a\)\(\ \)\(x\)\(\ \)\); 
cofexact = 
          Table[\(\@\(2  π\)\) \((n + 1)\)\ \(I\^n\) \(a\^\(-1\)\) 
              BesselJ[1 + n, a], {n, 0, nn - 1}]; 
{t1, 
            s1} = NChebyshevUListCoefficients[f[x], x, nn, 
              PrecisionGoal -> p] // 
            Timing; 

{t2, s2} = 
          Table[\(\@\(2\/π\)\) 
                NIntegrate[
                  Evaluate[
                    f[x] ChebyshevU[n, 
                        x] \((1 - x\^2)\)\^\(1/4\)], {x, \(-1\), 1}, 
                  PrecisionGoal -> p], {n, 0, nn - 1}] // 
            Timing; 
errornew = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s1]\^2]]; 
errorclassic = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s2]\^2]]; 
Print[\*""(Subscript[L, 2])-error->{classic,new,classic/new}->"", {errorclassic, errornew, 
            errorclassic/
              errornew}]; 
Print["Time->{classic,new,classic/new}->", {t2, t1, 
            t2/t1}]; 
ort = 
          Table[\(\@\(2\/π\)\) ChebyshevU[n, x], {n, 0, 
              nn - 1}]; 
s3 = 
          NChebyshevUSeries[f[x], x, nn, 
            PrecisionGoal -> p]; 
s4 = 
          s2 . \ ort; 
g1 = 
          Plot[Evaluate[Abs[s3 - f[x]] \((1 - x\^2)\)\^\(1/4\)], {x, \(-1\), 
              1}, PlotRange -> All, DisplayFunction -> Identity, 
            PlotStyle -> {Hue[1]}, 
            PlotLabel -> "New"]; 
g2 = 
          Plot[Evaluate[Abs[s4 - f[x]] \((1 - x\^2)\)\^\(1/4\)], {x, \(-1\), 
              1}, PlotRange -> All, DisplayFunction -> Identity, 
            PlotLabel -> "Classic"]; 
Show[
          GraphicsArray[{g2, 
              g1}]]\ 

], {nn, 35, 175, 
        20}]\)

ChebyshevT

Here we  calculate  the ChebyshevU  coefficients for the function  (E^(I a x )) (for this function we know exact coefficients)  and construct the corresponding  orthogonal decomposition. We compare two methods: the first (classic) is the direct evaluation of the corresponding integrals, the second (new) are  the programs NChebyshevTListCoefficients[] and  NChebyshevTSeries[]. We compare both the times and the precisions ((Subscript[L, 2])-error). Here we set PrecisionGoal=11 and 35≤N≤175. The new method is much more precise and fast than classic.

\(Do[\ 
Print["***** nn->", 
        nn]; 
Block[{p = 
            11}, 
Print["***** PrecisionGoal->", 
          p]; 
a = 5.231; 
f[x_] = 
          E\^\(\(I\)\(\ \)\(a\)\(\ \)\(x\)\(\ \)\); 
cofexact = \ 
          Prepend[Table[\(\@\(2  π\)\) \(I\^n\) BesselJ[n, a], {n, 1, 
                nn - 1}], \@π\ BesselJ[0, a]]; 
{t1, 
            s1} = NChebyshevTListCoefficients[f[x], x, nn, 
              PrecisionGoal -> p] // Timing; 
norm = 
          Flatten[{π, 
              Table[π/2, {i, 1, nn - 1}]}]; 
{t2, s2} = 
          1\/\@norm\ Table[
                NIntegrate[
                  Evaluate[
                    f[x] Cos[
                        n\ ArcCos[
                            x]] \((1 - x\^2)\)\^\(\(-1\)/2\)], {x, \(-1\), 
                    1}, PrecisionGoal -> p], {n, 0, nn - 1}] // 
            Timing; 
errornew = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s1]\^2]]; 
errorclassic = 
          Sqrt[Apply[Plus, 
              Abs[cofexact - s2]\^2]]; 
Print[\*""(Subscript[L, 2])-error->{classic,new,classic/new}->"", {errorclassic, errornew, 
            errorclassic/
              errornew}]; 
Print["Time->{classic,new,classic/new}->", {t2, t1, 
            t2/t1}]; 
ort = \(1\/\@norm\) 
            Table[Cos[n\ ArcCos[x]], {n, 0, nn - 1}]; 
s3 = 
          NChebyshevTSeries[f[x], x, nn, 
            PrecisionGoal -> p]; 
s4 = 
          s2 . \ ort; 
g1 = 
          Plot[Evaluate[
              Abs[s3 - f[x]] \((1 - x\^2)\)\^\(\(-1\)/4\)], {x, \(-1\), 1}, 
            PlotRange -> All, DisplayFunction -> Identity, 
            PlotStyle -> {Hue[1]}, 
            PlotLabel -> "New"]; 
g2 = 
          Plot[Evaluate[
              Abs[s4 - f[x]] \((1 - x\^2)\)\^\(\(-1\)/4\)], {x, \(-1\), 1}, 
            PlotRange -> All, DisplayFunction -> Identity, 
            PlotLabel -> "Classic"]; 
Show[
          GraphicsArray[{g2, 
              g1}]]

], {nn, 35, 175, 
        20}]\)