
(* :Name: System`Series` *)

(* :Title: InverseSeries and Series of Special Functions *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
	Defines InverseSeries[ ] using Newton's method symbolically.
	Provides series of special functions using special properties
	of those functions.  Also allows series expansions of special
	functions about singularities of certain simple types.
*)

(* :Context: System` *)

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1990  Wolfram Research, Inc.

	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :History:
	Updated by Jerry B. Keiper, December 1990.
	Version 1.2 by Jerry B. Keiper, 1988.
	InverseSeries[ ] was originally by Daniel R. Grayson.
*)

(* :Keywords:
	Series, InverseSeries, special functions.
*)

(* :Source:
	M. Abramowitz & I. Stegun, Handbook of Mathematical Functions With
	Formulas, Graphs, and Mathematical Tables, National Bureau of
	Standards, Applied Mathematics Series, Vol. 55.  Reprinted by
	Dover.
*)

(* :Mathematica Version: 2.0 *)

(* :Limitation:
	InverseSeries[ ] will not find the inverse series of a series
	expansion about a transcendental singularity.
*)

 (* --------------------- reversion of series --------------------------- *)

Begin["System`"]

Unprotect[ InverseSeries, Series, SeriesData]

Map[ Clear, {InverseSeries, Series, SeriesData}]

InverseSeries::usage = "InverseSeries[s] takes the series s generated by
Series, and gives a series for the inverse of the function represented by s.
InverseSeries[s, y] gives the inverse series in terms of the variable y."

Begin["InverseSeries`Private`"]

InverseSeries[f_SeriesData] := InverseSeries[f, f[[1]]];

InverseSeries[f_SeriesData, y_Symbol] :=
    Module[{s},
	    s /; (s = f; s[[2]] = 0; s = 1/InverseSeries[s, y];
		    Head[s] === SeriesData)
	] /; Head[f[[2]]] === DirectedInfinity

InverseSeries[f_SeriesData, y_Symbol] :=
	Which[
	    f[[4]] === 0,
		invser[f[[3,1]], 0, 1, f - f[[3,1]], y],
	    f[[4]] > 0,
		invser[0, 0, 1, f, y],
	    f[[4]] < 0,
		invser[DirectedInfinity[ ], 0, 1, 1/f, y]
	] /; (Length[f[[3]]] > If[ f[[4]] === 0, 1, 0] && IntegerQ[f[[4]]] &&
		Head[f[[1]]] === Symbol && (FreeQ[f, y] || f[[1]] === y) &&
		Head[f[[2]]] =!= DirectedInfinity)

InverseSeries[y0_. + (x_ b_. + mx0_.)^(a_) f_SeriesData, y_Symbol] :=
	invser[y0,a,b,f,y] /;
	(f[[1]] === x && Head[x] === Symbol && TrueQ[f[[2]] b == -mx0] &&
	FreeQ[y0,x] && FreeQ[mx0,x] && FreeQ[a,x] && FreeQ[b,x] &&
	(Length[f[[3]]] > 0) && (x === y ||
		(FreeQ[y0, y] && FreeQ[mx0,y] && FreeQ[a,y] && FreeQ[b,y] &&
			FreeQ[f,y])))

InverseSeries[y0_. + b_ (xmx0_)^(a_) f_SeriesData, y_Symbol] :=
    Module[{s = b f}, InverseSeries[y0 + (xmx0)^a s, y]] /; FreeQ[b, x]

invser[yc_, a_, b_, f_, y_] :=
    Module[{ba = f[[6]]/(a f[[6]]+f[[4]]), y0, z, s, w, newts},
	    y0 = If[Head[yc] === DirectedInfinity, 0, yc];
	    s = b z SeriesData[z,0,f[[3]]/f[[3,1]],0,f[[5]]-f[[4]],f[[6]]]^ba;
	    s[[3]] = Together[s[[3]]];
	    s = Series[Normal[s], {z, 0, Floor[(s[[5]]-1)/s[[6]]]}];
	    w = Series[(y-y0)/f[[3,1]], {y, y0, Length[s[[3]]]+1}]^ba;
	    newts = newton[s];
	    If[Head[w] === SeriesData, 
		w = ComposeSeries[newts, w] + f[[2]];
		If[Head[yc] === DirectedInfinity, w[[2]] = yc],
	      (* else *)
		(* Another horrible kludge in Series[ ].  This one is because
			ComposeSeries[ ] doesn't. *)
		w = (newts + f[[2]]) /. z -> ((y-y0)/f[[3,1]])^ba;
		];
	    w] /; Length[f[[3]]] > 0

newton[f_SeriesData] :=
    Module[{n=f[[5]], d=f[[6]], g, var=f[[1]], i, ff, gg, hh, kk},
	    g = SeriesData[var,0,{1/f[[3,1]]},d,n,d];
	    i = 2d;
	    While[ i < n,
		i *= 2;
		ff = f + SeriesData[var,0,{},i,i,d];
		gg = g + SeriesData[var,0,{},i,i,d];
		kk = (SeriesData[var,0,{1},d,i,d] - ComposeSeries[ff,gg])/
				ComposeSeries[D[ff,var],gg];
		kk[[5]] = n;
		g = g + kk
		];
	    g]

Protect[InverseSeries];

End[]  (* "InverseSeries`Private`" *)

 (* --------------------- Special Functions ---------------------------- *)

StieltjesGamma::usage =
"StieltjesGamma[n] is the nth Stieltjes constant and is related to the
Laurent series expansion of Zeta[s] about s == 1."

Unprotect[AiryAi,AiryAiPrime,AiryBi,AiryBiPrime,BesselI,BesselJ,BesselK,
BesselY,Beta,BetaRegularized,ChebyshevT,ChebyshevU,ExpIntegralE,ExpIntegralEi,
Factorial,FresnelC,FresnelS,Gamma,GammaRegularized,GegenbauerC,HermiteH,
Hypergeometric0F1,Hypergeometric0F1Regularized,Hypergeometric1F1,
Hypergeometric1F1Regularized,Hypergeometric2F1,Hypergeometric2F1Regularized,
HypergeometricU,LaguerreL,LegendreP,LegendreQ,LerchPhi,LogIntegral,LogGamma,
Pochhammer,PolyGamma,PolyLog,SphericalHarmonicY,Zeta,SinIntegral,CosIntegral,
SinhIntegral,CoshIntegral,SeriesData,Series];

Begin["SpecialFunctions`Series`Private`"]


OnQ[mess_] := MemberQ[{MessageName, String}, Head[mess]];

onSDsdatc = OnQ[SeriesData::sdatc];
If[onSDsdatc, Off[SeriesData::sdatc]];

 (* -----------------  utilities for SeriesData  -------------------------- *)

	(* Evaluate the limit of a series if the limit point is
	 the same as the expansion point. *)

SeriesData/: Literal[Limit[s_SeriesData, x_ -> a_]] :=
	faas[s] /; (simser[s] && s[[1]] === x && s[[2]] == a)
 
	(* Coefficient of a SeriesData object *)

SeriesData/: Literal[Coefficient[s_SeriesData, x___]] :=
	Coefficient[Normal[s], x] /; simser[s]

	(* CoefficientList of a SeriesData object *)

SeriesData/: Literal[CoefficientList[s_SeriesData, x___]] :=
	CoefficientList[Normal[s], x] /; simser[s]

 (* --------------------- powers of series ----------------------------- *)

(*
Commented out by DJW
SeriesData /: Power[s_SeriesData, v_?NumberQ] :=
	Module[{vv},(s^vv) /. vv -> v]
*)

SeriesData /:
Power[SeriesData[z_, z0_, cl_List, nmin_, nmax_, den_], v_] :=
    Module[{cla, a, na, p, exp = v nmin/den},
	If[Length[cl] == 0, a = 1; cla = cl, a = cl[[1]]; cla = cl/a];
	na = N[a];
	If[NumberQ[na] && Re[na] < 0,
	    (-a)^v If[Head[z0]===DirectedInfinity, (-z)^(-exp), (z0-z)^exp],
	    a^v If[Head[z0]===DirectedInfinity, z^(-exp), (z-z0)^exp]] *
		SeriesData[z, z0, cla, 0, nmax-nmin, den]^v
	] /; (nmin =!= 0 && !NumberQ[v])

 (* -------------------- Log expansion about Infinity ------------------ *)

  (* this should also deal with Log[a_. s_SeriesData], but then it couldn't
     be attached to SeriesData (SeriesData would be inside of a Times) *)

SeriesData /:
Log[SeriesData[z_, z0_DirectedInfinity, cl_List, nmin_, nmax_, den_]] :=
    Log[cl[[1]] z^(-nmin/den)] +
	ComposeSeries[Series[Log[z], {z, 1, mterms[1, nmin, nmax, den]}],
	    SeriesData[z, z0, cl/cl[[1]], 0, nmax-nmin, den]] /;
    Length[cl] > 0 && simser[z, cl]

 (* ------- Various simplification rules and support functions --------- *)

logrule = ((Log[n_?((NumberQ[#] && #<0)&)] + Log[-1 + y_]) ->
		(Log[-n] + Log[1-y]));

logsplit[Plus[logpart_ + rest___]] :=
	{logpart,Plus[rest]} /; (!FreeQ[logpart,Log] && FreeQ[{rest},Log]);

logsplit[rest_] := {0,rest} /; FreeQ[{rest},Log];

logsplit[logpart_] := {logpart,0} /; !FreeQ[logpart,Log];

intQ[m_] := (NumberQ[m] && (m == Round[m]));

    (* faa[] takes a coefficient list and the exponent associated with the
    first term in the coefficient list and returns the value of the series
    evaluated at the center of expansion *)

faas[s_SeriesData] := faa[s[[1]], s[[2]], s[[3]], s[[4]]] /; Length[s] == 6

faa[x_, a_, cl_, nmin_] :=
	If[nmin == 0,
		If[Length[cl] > 0, Limit[cl[[1]], x -> a], Indeterminate],
	    (* else *)
		If[nmin < 0, DirectedInfinity[ ], 0]];

    (* mterms[ ] takes a number (either 0 or nonzero: presumably the value of
    the derivative of our function f at the expansion point), nmin, nmax,
    and den and returns a guess as to the number of terms needed in the
    expansion of f w.r.t. a temporary variable in order to get the required
    number of terms in the composed series *)

mterms[fp_, s_SeriesData] := mterms[fp, s[[4]], s[[5]], s[[6]]]

mterms[fp_,nmin_,nmax_,den_] :=
	If[nmin == 0,
		nmax,
	    (* else *)
		If[nmin > 0,
			Floor[nmax/nmin],
		    (* else *)
			2 - Ceiling[nmax/nmin]
			]
		] + If[TrueQ[fp==0], 1, 0];

    (* simser[ ] tests whether the series is a simple series.  Series in x
    containing things like Log[x] in their coefficients cannot be handled
    properly by some of the functions below. *)

simser[s_SeriesData] :=
	Length[s] == 6 && (Symbol === Head[s[[1]]]) && FreeQ[s[[3]], s[[1]]]

simser[x_, cl_] := (Symbol === Head[x]) && FreeQ[cl, x]

Attributes[esssQ] = HoldRest

esssQ[f0_, s_] := If[Head[f0] === DirectedInfinity,
			Message[Series::esss, HoldForm[s]]; False,
			True]

 (* ----------------------- Airy functions ----------------------------- *)

	(* special cases: expansions about +/- Infinity *)

airyc[k_] := Module[{j}, Product[j, {j, 2k+1, 6k-1, 2}]/(216^k k!)];

airyd[k_] := -(6k+1) airyc[k]/(6k-1);

airycdsum[h_, zeta_, n_, i_] :=
    Module[{k}, {Sum[i^(k/2) h[k] zeta^(-k), {k, 0, 2n/3 + 2, 2}],
		Sum[i^((k-1)/2) h[k] zeta^(-k), {k, 1, 2n/3 + 2, 2}]}
	];

AiryAi /: Series[AiryAi[z_], {z_, Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{s1, s2} = airycdsum[airyc, zeta, n, 1];
	s1 = 1/2 Pi^(-1/2) z^(-1/4) E^(-zeta) Series[s1-s2, {z, Infinity, n}];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryAiPrime /: Series[AiryAiPrime[z_], {z_, Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{s1, s2} = airycdsum[airyd, zeta, n, 1];
	s1 = -1/2 Pi^(-1/2) z^(1/4) E^(-zeta) Series[s1-s2, {z, Infinity, n}];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryBi /: Series[AiryBi[z_], {z_, Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{s1, s2} = airycdsum[airyc, zeta, n, 1];
	s1 = Pi^(-1/2) z^(-1/4) E^zeta Series[s1+s2, {z, Infinity, n}];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryBiPrime /: Series[AiryBiPrime[z_], {z_, Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{s1, s2} = airycdsum[airyd, zeta, n, 1];
	s1 = Pi^(-1/2) z^(1/4) E^zeta Series[s1+s2, {z, Infinity, n}];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryAi /: Series[AiryAi[z_], {z_, -Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	onGivar = OnQ[General::ivar];
	If[onGivar, Off[General::ivar]];
	{s1, s2} = airycdsum[airyc, zeta, n, -1];
	s1 = Sin[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
	s2 = Cos[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
	s1 = Pi^(-1/2) (-z)^(-1/4) (s1 - s2);
	If[onGivar, On[General::ivar]];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryAiPrime /: Series[AiryAiPrime[z_], {z_, -Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	onGivar = OnQ[General::ivar];
	If[onGivar, Off[General::ivar]];
	{s1, s2} = airycdsum[airyd, zeta, n, -1];
	s1 = Cos[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
	s2 = Sin[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
	s1 = -Pi^(-1/2) (-z)^(1/4) (s1 + s2);
	If[onGivar, On[General::ivar]];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryBi /: Series[AiryBi[z_], {z_, -Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	onGivar = OnQ[General::ivar];
	If[onGivar, Off[General::ivar]];
	{s1, s2} = airycdsum[airyc, zeta, n, -1];
	s1 = Cos[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
	s2 = Sin[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
	s1 = Pi^(-1/2) (-z)^(-1/4) (s1 + s2);
	If[onGivar, On[General::ivar]];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

AiryBiPrime /: Series[AiryBiPrime[z_], {z_, -Infinity, n_Integer}] :=
    Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	onGivar = OnQ[General::ivar];
	If[onGivar, Off[General::ivar]];
	{s1, s2} = airycdsum[airyd, zeta, n, -1];
	s1 = Sin[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
	s2 = Cos[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
	s1 = Pi^(-1/2) (-z)^(1/4) (s1 - s2);
	If[onGivar, On[General::ivar]];
	If[onSesss, On[Series::esss]];
	s1
	] /; n >= 0

	(* generic cases *)

AiryAi /: Series[AiryAi[g_],{z_,aa_,nn_Integer}] :=
	AiryAi[Series[g,{z,aa,nn}]] /; nn >= 0

AiryAiPrime /: Series[AiryAiPrime[g_],{z_,aa_,nn_Integer}] :=
	AiryAiPrime[Series[g,{z,aa,nn}]] /; nn >= 0

AiryBi /: Series[AiryBi[g_],{z_,aa_,nn_Integer}] :=
	AiryBi[Series[g,{z,aa,nn}]] /; nn >= 0

AiryBiPrime /: Series[AiryBiPrime[g_],{z_,aa_,nn_Integer}] :=
	AiryBiPrime[Series[g,{z,aa,nn}]] /; nn >= 0

AiryAi /: AiryAi[s_SeriesData] :=
    Module[{aa = faas[s]},
	diffeqAiry[AiryAi,False,s,aa] /; esssQ[aa, AiryAi[s]]] /; simser[s]

AiryAiPrime /: AiryAiPrime[s_SeriesData] :=
    Module[{aa = faas[s]},
	diffeqAiry[AiryAi,True,s,aa] /; esssQ[aa, AiryAiPrime[s]]] /; simser[s]

AiryBi /: AiryBi[s_SeriesData] :=
    Module[{aa = faas[s]},
	diffeqAiry[AiryBi,False,s,aa] /; esssQ[aa, AiryBi[s]]] /; simser[s]

AiryBiPrime /: AiryBiPrime[s_SeriesData] :=
    Module[{aa = faas[s]},
	diffeqAiry[AiryBi,True,s,aa] /; esssQ[aa, AiryBiPrime[s]]] /; simser[s]

diffeqAiry[head_, der_, s_SeriesData, aa_] :=
    Module[{i, m, f0, f1, f, z = s[[1]], w},
	If[SameQ[aa,Indeterminate], Return[SeriesData[z,aa,{},1,1,1]]];
	m = mterms[1, s] + 1;
	f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i,{i,2,m}],{z,aa,m}];
	f = MapAt[Together,
		(f /. Solve[D[f,{z,2}]==z f,Table[w[i],{i,2,m}]])[[1]],
		{3}];
	f = f /. {f0 -> head[aa], f1 -> head'[aa]};
	If[der, f = D[f, z]];
	ComposeSeries[f, s]
	];

 (* ------------------------ Bessel functions ----------------------------- *)

	(* special cases: expansion about Infinity *)

bessPQ[v_, z_, n_, i_] :=
    Module[{j, k, mu = 4v^2, p, q},
	(* AS 9.2.9, 9.2.10 *)
	p = Sum[i^(k/2) Product[mu-j^2,{j,1,2k-1,2}]/
		(k! 8^k z^k),{k,0,n+1,2}];
	q = Sum[i^((k-1)/2) Product[mu-j^2,{j,1,2k-1,2}]/
		(k! 8^k z^k),{k,1,n+1,2}];
	{p, q}
	];

BesselJ /: Series[BesselJ[v_, z_], {z_, Infinity, n_Integer}] :=
    Module[{p, q, x = z - (v/2 + 1/4)Pi, onSesss},
	(* AS 9.2.5 *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{p, q} = bessPQ[v, z, n, -1];
	p = Series[p, {z, Infinity, n}];
	q = Series[q, {z, Infinity, n}];
	p = Sqrt[2/(Pi z)] (p Cos[x] - q Sin[x]);
	If[onSesss, On[Series::esss]];
	p
	] /; n >= 0

BesselY /: Series[BesselY[v_, z_], {z_, Infinity, n_Integer}] :=
    Module[{p, q, x = z - (v/2 + 1/4)Pi, onSesss},
	(* AS 9.2.6 *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{p, q} = bessPQ[v, z, n, -1];
	p = Series[p, {z, Infinity, n}];
	q = Series[q, {z, Infinity, n}];
	p = Sqrt[2/(Pi z)] (p Sin[x] + q Cos[x]);
	If[onSesss, On[Series::esss]];
	p
	] /; n >= 0

BesselI /: Series[BesselI[v_, z_], {z_, Infinity, n_Integer}] :=
    Module[{p, q, onSesss},
	(* AS 9.7.1, corrected *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{p, q} = bessPQ[v, z, n, 1];
	p = Series[p, {z, Infinity, n}];
	q = Series[q, {z, Infinity, n}];
	p = Sqrt[1/(2 Pi z)] (E^z (p-q) + E^(I (v + 1/2) Pi - z) (p+q));
	If[onSesss, On[Series::esss]];
	p
	] /; n >= 0

BesselK /: Series[BesselK[v_, z_], {z_, Infinity, n_Integer}] :=
    Module[{p, q, onSesss},
	(* AS 9.7.2 *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	{p, q} = bessPQ[v, z, n, 1];
	p = Series[p+q, {z, Infinity, n}];
	p = Sqrt[Pi/(2 z)] E^(-z) p;
	If[onSesss, On[Series::esss]];
	p
	] /; n >= 0

	(* generic case *)

besser[v_,z_,a_,n_,fa_,fpa_,m_] := 
    Module[{i, f0, f1, df, f},
	f = Series[f0 + f1(z-a) + Sum[w[i] (z-a)^i, {i,2,n}],{z,a,n}];
	df = D[f,z];
	f = MapAt[Expand,
		(f /. Solve[z^2 D[df,z] + z df + (m z^2 - v^2) f == 0,
			Table[w[i],{i,2,n}]])[[1]],
		{3}];
	f /. {f0 -> fa, f1 -> fpa}
	];

BesselI /: Series[BesselI[v_,g_],{z_,aa_,nn_Integer}] :=
	BesselI[v,Series[g,{z,aa,nn}]] /; nn >= 0

BesselI /: Literal[BesselI[v_, zser_SeriesData]] :=
    Module[{m, fp, f, aa = faas[zser], z},
	z = zser[[1]];
	If[TrueQ[aa==0],
		m = mterms[1, zser];
		If[NumberQ[v], m -= Floor[Re[v]]];
		f = Series[Hypergeometric0F1Regularized[v+1,z^2/4], {z,aa,m}];
		(zser/2)^v ComposeSeries[f, zser],
	    (* else *)
		fp = D[BesselI[v,z],z] /. z->aa;
		m = mterms[fp, zser];
		f = besser[v,z,aa,m,BesselI[v,aa],fp,-1];
		ComposeSeries[f, zser]
		] /;
	esssQ[aa, BesselI[v, zser]]
	] /; simser[zser]

BesselJ /: Series[BesselJ[v_,g_],{z_,aa_,nn_Integer}] :=
	BesselJ[v,Series[g,{z,aa,nn}]] /; nn >= 0

BesselJ /: Literal[BesselJ[v_, zser_SeriesData]] :=
    Module[{m, fp, f, aa = faas[zser], z},
	z = zser[[1]];
	If[TrueQ[aa==0],
		m = mterms[1, zser];
		If[NumberQ[v], m -= Floor[Re[v]]];
		f = Series[Hypergeometric0F1Regularized[v+1,-z^2/4], {z,aa,m}];
		(zser/2)^v ComposeSeries[f, zser],
	    (* else *)
		fp = D[BesselJ[v,z],z] /. z->aa;
		m = mterms[fp, zser];
		f = besser[v,z,aa,m,BesselJ[v,aa],fp,1];
		ComposeSeries[f,zser]
		] /;
	esssQ[aa, BesselJ[v, zser]]
	] /; simser[zser]

BesselK /: Series[BesselK[v_,g_],{z_,aa_,nn_Integer}] :=
	BesselK[v,Series[g,{z,aa,nn}]] /; nn >= 0

BesselK /: Literal[BesselK[v_, zser_SeriesData]] :=
    Module[{m, k, nu = v, fp, f, aa = faas[zser], z},
	z = zser[[1]];
	(If[Re[nu] < 0, nu = -nu];
	If[TrueQ[aa==0],
	    If[intQ[nu],
		m = mterms[1, zser];
		f = Series[1/2 (2/z)^nu *
			Sum[(nu-k-1)! (-z^2/4)^k/k!, {k,0,nu-1}] -
			(-1)^nu (Log[z/2] BesselI[nu,z] - 1/2(z/2)^nu *
			Sum[(PolyGamma[k+1]+PolyGamma[k+1+nu])(z^2/4)^k/
			(k!(nu+k)!), {k,0,m-Abs[Round[nu]]}]),{z,0,m}];
		ComposeSeries[f,zser],
	      (* else (!intQ[nu]) *)
		Pi(BesselI[-nu,zser] - BesselI[nu,zser])/(2 Sin[nu Pi])
		],
	  (* else (aa != 0)*)
	    fp = D[BesselK[nu,z],z] /. z->aa;
	    m = mterms[fp, zser];
	    f = besser[nu,z,aa,m,BesselK[nu,aa],fp,-1];
	    ComposeSeries[f,zser]
	    ]) /;
	esssQ[aa, BesselK[v, zser]]
	] /; simser[zser]

BesselY /: Series[BesselY[v_,g_],{z_,aa_,nn_Integer}] :=
	BesselY[v,Series[g,{z,aa,nn}]] /; nn >= 0

BesselY /: Literal[BesselY[v_, zser_SeriesData]] :=
    Module[{m, k, nu = v, fp, f, aa = faas[zser], sign = 1, z},
	z = zser[[1]];
	(If[intQ[nu/2] && nu < 0, nu = -nu; sign = (-1)^nu];
	If[TrueQ[aa==0],
	    If[intQ[nu],
		m = mterms[1, zser];
		f = Series[-1/Pi (2/z)^nu *
			Sum[(nu-k-1)! (z^2/4)^k/k!, {k,0,nu-1}] +
			2/Pi Log[z/2] BesselJ[nu,z] - 1/Pi (z/2)^nu *
			Sum[(PolyGamma[k+1]+PolyGamma[k+1+nu]) *
			(-z^2/4)^k/(k!(nu+k)!), {k,0,m-Abs[Round[nu]]}],
			{z,0,m}];
		ComposeSeries[f,zser],
	      (* else (!intQ[nu]) *)
		(BesselJ[v,zser] Cos[v Pi] - BesselJ[-v,zser])/Sin[v Pi]
		],
	  (* else (aa != 0)*)
	    fp = D[BesselY[nu,z],z] /. z->aa;
	    m = mterms[fp, zser];
	    f = besser[nu,z,aa,m,BesselY[nu,aa],fp,1];
	    ComposeSeries[f,zser]
	    ]) /;
	esssQ[aa, BesselY[v, zser]]
	] /; simser[zser]

 (* ------------------------- Beta functions ----------------------------- *)

betaseries[sf1_,sf2_,a_,b_,regflag_,z_,z0_] :=
    Module[{aa=a,bb=b,m1,m2,f1=0,f2=0,f01,f02,t,anint,bnint,logs1,logs2},
	If[intQ[a],aa = Round[aa]];
	anint = IntegerQ[aa] && aa<=0;
	If[intQ[b],bb = Round[bb]];
	bnint = IntegerQ[bb] && bb<=0;
	If[sflag1 = (Head[sf1] === SeriesData),
		f01 = faas[sf1];
		m1 = mterms[1, sf1],
	    (* else *)
		f01 = sf1;
		m1 = Infinity;
		];
	If[sflag2 = (Head[sf2] === SeriesData),
		f02 = faas[sf2];
		m2 = mterms[1,sf2],
	    (* else *)
		f02 = sf2;
		m2 = Infinity;
		];
	logs1 = sflag1 && TrueQ[f01==0];
	logs2 = sflag2 && TrueQ[f02==1];
	If[sflag1, m1 += Max[If[TrueQ[f01==0], -Round[aa],0],
		If[TrueQ[f01==1], -Round[bb],0]]];
	If[sflag2, m2 += Max[If[TrueQ[f02==0], -Round[aa],0],
		If[TrueQ[f02==1], -Round[bb],0]]];
	If[anint && bnint && (logs1 || logs2),
		If[regflag, Return[0]];
		(* this is slow, but I don't know how else to do it *)
		f1 = Integrate[t^(aa-1)(1-t)^(bb-1),t] /. Log[-1+t]->Log[1-t];
		f2 = If[sflag2, ComposeSeries[Series[f1,{t,f02,m2}],sf2],
		    (* else *) f1 /. t->f02];
		f1 = If[sflag1, ComposeSeries[Series[f1,{t,f01,m1}],sf1],
		    (* else *) f1 /. t->f01];
		Return[f2-f1]
		];
	If[anint && logs1,
		Return[betahyper[1-f02,1-sf2,m2,1-f01,1-sf1,m1,bb,aa,regflag]]];
	If[bnint && logs2,
		Return[betahyper[f01,sf1,m1,f02,sf2,m2,aa,bb,regflag]]];
	If[sflag1, {f01,f1} = betser0[aa,bb,sf1,f01,m1]];
	If[sflag2, {f02,f2} = betser0[aa,bb,sf2,f02,m2]];
	If[regflag,
		BetaRegularized[f01,f02,a,b] + (f2 - f1)/Beta[a,b],
		Beta[f01,f02,a,b] + f2 - f1
		]
	]

betahyper[f01_,sf1_,m1_,f02_,sf2_,m2_,a_,b_,regflag_] :=
    Module[{f1,f2,t,beta = Beta[a,b]},
	If[regflag && (Head[beta] === DirectedInfinity), Return[0]];
	(* AS 6.6.8; a is not a non-positive integer! *)
	f2 = ComposeSeries[t^a/a Series[Hypergeometric2F1[
		a,1-b,a+1,t],{t,f02,m2}],sf2];
	f1 = If[!SameQ[Head[sf1],SeriesData],
		f01^a/a Hypergeometric2F1[a,1-b,a+1,f01],
	    (* else *)
		ComposeSeries[t^a/a Series[Hypergeometric2F1[
			a,1-b,a+1,t],{t,f01,m1}],sf1]];
	If[regflag, (f2 - f1)/beta, f2 - f1] /. logrule
	];

betser0[a_,b_,sf_,f0_,m_] :=
    Module[{t, f},
	f = Integrate[Series[(t+f0)^(a-1)(1-f0-t)^(b-1),{t,0,m}],t];
	{f0,ComposeSeries[f, sf-f0]}
	];

Beta /: (* complete beta function *)
Series[Beta[f1z_,f2z_],{z_,aa_,nn_Integer}] :=
    Module[{sf1 = Series[f1z,{z,aa,nn}],sf2 = Series[f2z,{z,aa,nn}]},
	Gamma[sf1] Gamma[sf2] / Gamma[sf1+sf2]] /; nn >= 0

Beta /: (* complete beta function *)
Literal[Beta[sf1_SeriesData, sf2_SeriesData]]:=
	Gamma[sf1] Gamma[sf2] / Gamma[sf1+sf2] /;
	(simser[sf1] && simser[sf2] &&
	sf1[[1]] === sf2[[1]] && sf1[[2]] == sf2[[2]])

Beta /: (* complete beta function *)
Literal[Beta[sf_SeriesData, a_]]:=
    Module[{sa},
	sa = Series[a, {sf[[1]], sf[[2]],
			Round[(sf[[5]] - Min[0, sf[[4]]])/sf[[6]]]}];
	Gamma[sf] Gamma[sa] / Gamma[sf+sa]
	] /; simser[sf]

Beta /: (* complete beta function *)
Literal[Beta[a_, sf_SeriesData]] :=
    Module[{sa},
	sa = Series[a, {sf[[1]], sf[[2]],
			Round[(sf[[5]] - Min[0, sf[[4]]])/sf[[6]]]}];
	Gamma[sf] Gamma[sa] / Gamma[sf+sa]
	] /; simser[sf]

Beta /: (* incomplete beta function *)
Series[Beta[f1z_:0,f2z_,a_,b_],{z_,aa_,nn_Integer}] :=
	Beta[Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}],a,b] /;
		(nn >= 0 && FreeQ[a,z] && FreeQ[b,z])

Beta /: (* incomplete beta function *)
Literal[Beta[s_SeriesData, a_, b_]] :=
	betaseries[0, s, a, b, False, s[[1]], s[[2]]] /;
	(simser[s] && FreeQ[a,s[[1]]] && FreeQ[b,s[[1]]])

Beta /: (* incomplete beta function *)
Literal[Beta[s1_SeriesData, s2_SeriesData, a_, b_]] :=
	betaseries[s1, s2, a, b, False, s1[[1]], s1[[2]]] /;
	(simser[s1] && simser[s2] && s1[[1]] === s2[[1]] &&
	s1[[2]] == s2[[2]] && FreeQ[a,s1[[1]]] && FreeQ[b,s1[[1]]])

Beta /: (* incomplete beta function *)
Literal[Beta[f1_SeriesData,f2_,a_,b_]] :=
    betaseries[f1,
	Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
	a,b,False,f1[[1]],f1[[2]]] /;
    (simser[f1] && FreeQ[a,f1[[1]]] && FreeQ[b,f1[[1]]])

Beta /: (* incomplete beta function *)
Literal[Beta[f1_,f2_SeriesData,a_,b_]] :=
    betaseries[
	Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
	f2,a,b,False,f2[[1]],f2[[2]]] /;
	(simser[f2] && FreeQ[a,f2[[1]]] && FreeQ[b,f2[[1]]])

BetaRegularized /: (* incomplete beta function *)
Series[BetaRegularized[f1z_:0,f2z_,a_,b_],{z_,aa_,nn_Integer}] :=
	BetaRegularized[Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}],a,b] /;
	(nn >= 0 && FreeQ[a,z] && FreeQ[b,z])

BetaRegularized /: (* incomplete beta function *)
Literal[BetaRegularized[s2_SeriesData, a_, b_]] :=
	betaseries[0, s2, a, b, True, s2[[1]], s2[[2]]] /;
	(simser[s2] && FreeQ[a,s2[[1]]] && FreeQ[b,s2[[1]]])

BetaRegularized /: (* incomplete beta function *)
Literal[BetaRegularized[s1_SeriesData, s2_SeriesData, a_, b_]] :=
	betaseries[s1, s2, a, b, True, s1[[1]], s1[[2]]] /;
	(simser[s1] && simser[s2] && s1[[1]] === s2[[1]] &&
	s1[[2]] == s2[[2]] && FreeQ[a,s1[[1]]] && FreeQ[b,s1[[1]]])

BetaRegularized /: (* incomplete beta function *)
Literal[BetaRegularized[f1_SeriesData,f2_,a_,b_]] :=
    betaseries[f1,
	Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
	a,b,True,f1[[1]],f1[[2]]] /;
    (simser[f1] && FreeQ[a,f1[[1]]] && FreeQ[b,f1[[1]]])

BetaRegularized /: (* incomplete beta function *)
Literal[BetaRegularized[f1_,f2_SeriesData,a_,b_]] :=
    betaseries[
	Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
	f2,a,b,True,f2[[1]],f2[[2]]] /;
    (simser[f2] && FreeQ[a,f2[[1]]] && FreeQ[b,f2[[1]]])

 (* --------------------- Binomial function ---------------------------- *)

Series[Binomial[x_,y_],{z_,aa_,nn_Integer}] :=
    Module[{ys = Series[y, {z, aa, nn}]},
	1/(ys Beta[ys, 1+Series[x, {z, aa, nn}] - ys])] /; (nn >= 0)

SeriesData /:
Literal[Binomial[xs_SeriesData, y_]] :=
    Module[{ys = Series[y, {xs[[1]], xs[[2]],
		Round[(xs[[5]] - Min[0, xs[[4]]])/xs[[6]]]}]},
	1/(ys Beta[ys, xs + (1-ys)])
	] /; simser[xs]

SeriesData /:
Literal[Binomial[x_, ys_SeriesData]] :=
    Module[{xs = Series[x, {ys[[1]], ys[[2]],
		Round[(ys[[5]] - Min[0, ys[[4]]])/ys[[6]]]}]},
	1/(ys Beta[ys, 1+xs-ys])
	] /; simser[ys]

SeriesData /:
Literal[Binomial[xs_SeriesData, ys_SeriesData]] :=
	1/(ys Beta[ys, 1+xs-ys]) /; simser[xs] && simser[ys]

 (* -------------------- Chebyshev functions --------------------------- *)

ChebyshevT /: Series[ChebyshevT[n_,fz_],{z_,aa_,nn_Integer}] :=
	ChebyshevT[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])

ChebyshevT /: ChebyshevT[n_,s_SeriesData] := Cos[n ArcCos[s]] /; simser[s]

ChebyshevU /: Series[ChebyshevU[n_,fz_],{z_,aa_,nn_Integer}] :=
	ChebyshevU[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])

ChebyshevU /: ChebyshevU[n_,s_SeriesData] :=
    Module[{theta = ArcCos[s]}, Sin[(n+1) theta]/Sin[theta]] /; simser[s]

 (* ------------------- Exponential integrals -------------------------- *)

	(* special cases: expansions about Infinity *)

ExpIntegralE /: Series[ExpIntegralE[n_, z_], {z_, Infinity, nn_Integer}] :=
    Module[{j, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	j = E^(-z) Series[Sum[(-1)^j Pochhammer[n, j] z^(-j-1), {j, 0, nn}],
			{z, Infinity, nn}];
	If[onSesss, On[Series::esss]];
	j
	] /; (nn >= 0 && FreeQ[n, z])

ExpIntegralEi /: Series[ExpIntegralEi[z_], {z_, Infinity, nn_Integer}] :=
    Module[{j, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	j = E^z Series[Sum[j! z^(-j-1), {j, 0, nn}], {z, Infinity, nn}];
	If[onSesss, On[Series::esss]];
	j
	] /; nn >= 0

	(* generic cases *)

ExpIntegralE /: Series[ExpIntegralE[n_,fz_],{z_,aa_,nn_Integer}] :=
	ExpIntegralE[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])

ExpIntegralE /:
Literal[ExpIntegralE[n_, fz_SeriesData]] :=
	fz^(n-1) Gamma[1-n,fz] /; (simser[fz] && FreeQ[n,fz[[1]]])

ExpIntegralEi /: Series[ExpIntegralEi[fz_],{z_,aa_,nn_Integer}] :=
	ExpIntegralEi[Series[fz,{z,aa,nn}]] /; nn >= 0

ExpIntegralEi /:
Literal[ExpIntegralEi[s_SeriesData]] :=
    Module[{f0 = faas[s], m = mterms[1, s], temp, fs},
	(fs = If[TrueQ[f0==0],
		EulerGamma+Log[temp]+Sum[temp^k/(k k!),{k,m}] + O[temp]^(m+1),
	    (* else *)
		Integrate[Series[Exp[temp+f0]/(temp+f0),{temp,0,m}],temp] +
			+ ExpIntegralEi[f0]
		];
	ComposeSeries[fs, s-f0]) /; esssQ[f0, ExpIntegralEi[s]]
	] /; simser[s]

 (* ----------------- Factorial and Gamma  functions ---------------------- *)

	(* special cases: expansions about Infinity *)

Gamma /: Series[Gamma[z_],{z_, Infinity, nn_Integer}] :=
    Module[{s, i, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s = Sum[BernoulliB[i]/(i(i-1) z^(i-1)), {i, 2, nn+2, 2}];
	s = z^z E^(-z) Sqrt[2 Pi/z] E^Series[s, {z, Infinity, nn}];
	If[onSesss, On[Series::esss]];
	s
	] /; nn >= 0

Factorial/: Series[Factorial[z_],{z_, Infinity, nn_Integer}] :=
    Module[{s, i, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s = Sum[BernoulliB[i]/(i(i-1) (z+1)^(i-1)), {i, 2, nn+2, 2}];
	s = (z+1)^z E^(-z-1) Sqrt[2 Pi (z+1)] E^Series[s, {z, Infinity, nn}];
	If[onSesss, On[Series::esss]];
	s
	] /; nn >= 0

	(* generic cases *)

gammaseries[a_,sf1_,sf2_,regflag_,z_,z0_] :=
    Module[{f1,f2,f01,f02},
	If[Head[sf1] === SeriesData,
		{f01,f1} = gamser0[a,sf1], f01 = sf1; f1 = 0];
	If[Head[sf2] === SeriesData,
		{f02,f2} = gamser0[a,sf2], f02 = sf2; f2 = 0];
	If[f01 === $Failed || f02 === $Failed,
		HoldForm[Head[a, sf1, sf2]] /.
			Head -> If[regflag, GammaRegularized, Gamma],
	    (* else *)
		If[regflag, GammaRegularized[a,f01,f02] + (f2 - f1)/Gamma[a],
			Gamma[a,f01,f02] + f2 - f1]
		]
	];

gamser0[a_, sf_] :=
    Module[{t, f, ma, i, m = mterms[1,sf], f0 = faas[sf]},
	If[Head[f0] === DirectedInfinity,
		If[TrueQ[f0 == Infinity],
			Return[{f0, 0}],
		    (* else *)
			Message[Series::esss, HoldForm[Gamma[a, sf]]];
			Return[{$Failed, $Failed}]
			]
		];
	If[TrueQ[f0==0] && intQ[a] && a <= 0,
		ma = -Round[a];
		m += ma;
		f = Sum[If[i==ma,0,(-1)^i t^(i-ma)/((i-ma)i!)],{i,0,m}] +
				O[t]^(m+1);
		f = ComposeSeries[f,sf] + (-1)^ma((Log[sf] /. logrule) -
			PolyGamma[ma+1])/ma!;
		{Infinity,f},
	    (* else *)
		f = Integrate[Series[Exp[-(t+f0)] (t+f0)^(a-1), {t,0,m}],t];
		{f0,ComposeSeries[f, sf-f0]}
		]
	];

Series[Factorial[fz_],{z_,aa_,nn_Integer}] :=
	Gamma[Series[fz+1,{z,aa,nn}]] /; nn >= 0

SeriesData /: Literal[Factorial[s_SeriesData]] := Gamma[1+s] /; simser[s]

Series[Gamma[fz_],{z_,aa_,nn_Integer}] :=
	Gamma[Series[fz,{z,aa,nn}]] /; nn >= 0

SeriesData /:  (* complete gamma function *)
Literal[Gamma[fz_SeriesData]] :=
    Module[{f0 = faas[fz], i},
	If[NumberQ[f0] && f0 < 1/2,
		Pi/(Sin[Pi fz] Gamma[1-fz]),
	    (* else *)
		Sum[Derivative[i][Gamma][f0](fz-f0)^i/(i!),{i,0,mterms[1,fz]}]
		] /;
	esssQ[f0, Gamma[fz]]
	] /; simser[fz]

Series[Gamma[a_,f1z_,f2z_:Infinity],{z_,aa_,nn_Integer}] :=
	Gamma[a,Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[a,z])

SeriesData /:  (* incomplete gamma function *)
Literal[Gamma[a_, s_SeriesData]]:=
	gammaseries[a,s,Infinity,False,s[[1]],s[[2]]] /;
	(simser[s] && FreeQ[a,s[[1]]])

SeriesData /:  (* incomplete gamma function *)
Literal[Gamma[a_, s1_SeriesData, s2_SeriesData]]:=
    gammaseries[a, s1, s2, False, s1[[1]], s1[[2]]] /;
	(FreeQ[a,z] && s1[[1]] === s1[[1]] && s1[[2]] == s2[[2]] &&
	simser[z, cl1] && simser[z, cl2])

SeriesData /:  (* incomplete gamma function *)
Literal[Gamma[a_,f1_,f2_SeriesData]] :=
    gammaseries[a,
	Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
	f2, False, f2[[1]], f2[[2]]] /;
    (simser[f2] && FreeQ[a,f2[[1]]])

SeriesData /:  (* incomplete gamma function *)
Literal[Gamma[a_, f1_SeriesData, f2_:Infinity]] :=
    gammaseries[a, f1,
	Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
	False, f1[[1]], f1[[2]]] /;
    (simser[f1] && FreeQ[a,f1[[1]]])

GammaRegularized /: (* incomplete gamma function *)
Series[GammaRegularized[a_,f1z_,f2z_:Infinity],{z_,aa_,nn_Integer}] :=
	GammaRegularized[a,Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[a,z])

GammaRegularized /: (* incomplete gamma function *)
Literal[GammaRegularized[a_, s_SeriesData]]:=
	gammaseries[a,s,Infinity,True,s[[1]],s[[2]]] /;
	(simser[s] && FreeQ[a,s[[1]]])

GammaRegularized /: (* incomplete gamma function *)
Literal[GammaRegularized[a_, s1_SeriesData, s2_SeriesData]]:=
    gammaseries[a,s1,s2,True,s1[[1]],s1[[2]]] /;
	(FreeQ[a,z] && s1[[1]] === s2[[1]] && s1[[2]] == s2[[2]] &&
	simser[z,cl1] && simser[z,cl2])

GammaRegularized /: (* incomplete gamma function *)
Literal[GammaRegularized[a_, f1_, f2_SeriesData]]:=
    gammaseries[a,
	Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
	f2, True, f2[[1]], f2[[2]]] /;
    (simser[f2] && FreeQ[a,f2[[1]]])

GammaRegularized /: (* incomplete gamma function *)
Literal[GammaRegularized[a_, f1_SeriesData, f2_:Infinity]] :=
    gammaseries[a, f1,
	Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
	True, f1[[1]], f1[[2]]] /;
    (simser[f1] && FreeQ[a,f1[[1]]])

 (* ----------------- Gegenbauer polynomials --------------------------- *)

GegenbauerC /: Series[GegenbauerC[n_,m_,fz_],{z_,aa_,nn_Integer}] :=
	GegenbauerC[n,m,Series[fz,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[n,z] && FreeQ[m,z])

GegenbauerC /:
Literal[GegenbauerC[n_, 0, s_SeriesData]] :=
	2/n ChebyshevT[n,s] /; (simser[s] && FreeQ[n,s[[1]]])

GegenbauerC /:
Literal[GegenbauerC[n_, m_, s_SeriesData]] :=
	Gamma[n+2m] Sqrt[Pi]/(n! 2^(2m-1) Gamma[m]) *
		Hypergeometric2F1Regularized[-n, n+2m, m+1/2, (1-s)/2] /;
	(simser[s] && FreeQ[n,s[[1]]] && FreeQ[m,s[[1]]])

 (* --------------------- Hermite polynomials --------------------------- *)

HermiteH /: Series[HermiteH[n_,fz_],{z_,aa_,nn_Integer}] :=
	HermiteH[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])

HermiteH /: Literal[HermiteH[n_, fz_SeriesData]] :=
	2^n fz HypergeometricU[(1-n)/2,3/2,fz^2] /;
	(simser[fz] && FreeQ[n,fz[[1]]])

 (* ------------------- Hypergeometric functions ------------------------- *)

	(* utility for expansions about Infinity *)

hyper2F0r[a_, b_, z_, nn_] :=
    Module[{n}, Sum[Pochhammer[a,n] Pochhammer[b,n] z^(-n)/n!, {n, 0, nn}]];

	(* utility for distinction between regularized and unregularized *)

regden[a_, i_, regflag_] := If[regflag, Gamma[a+i], Pochhammer[a, i]]

 (* --- Hypergeometric 0F1 --- *)

	(* special case: expansions about Infinity *)

hyper0F1inf[a_, z_, nn_, regflag_] :=
    Module[{s1, s2, zz = 4 Sqrt[z], aa = a - 1/2, onSesss},
	(* based on AS 13.5.1, with 0F1[a,z] == E^(-zz/2) 1F1[aa, 2aa, zz] *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s1 = Series[hyper2F0r[aa, 1-aa, -zz, nn], {z, Infinity, nn}];
	s2 = Series[hyper2F0r[aa, 1-aa, zz, nn], {z, Infinity, nn}];
	s1 = (E^(I Pi aa-zz/2)s1+E^(zz/2)s2) zz^(-aa) 4^(a-1)
				If[regflag, 1, Gamma[a]]/Sqrt[Pi];
	If[onSesss, On[Series::esss]];
	s1
	] /; (nn >= 0 && FreeQ[a, z])

Hypergeometric0F1 /:
Series[Hypergeometric0F1[a_, z_], {z_, Infinity, nn_Integer}] :=
    hyper0F1inf[a, z, nn, False (* not regularized *)] /;
	(nn >= 0 && FreeQ[a, z])

Hypergeometric0F1Regularized /:
Series[Hypergeometric0F1Regularized[a_, z_], {z_, Infinity, nn_Integer}] :=
    hyper0F1inf[a, z, nn, True (* regularized *)] /; (nn >= 0 && FreeQ[a, z])

	(* generic case *)

diffeq0F1[a_,fz_,aa_,z_,m_,head_] :=
    Module[{i, f0, f1, df, f},
	f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i, {i,2,m}],{z,aa,m}];
	df = D[f,z];
	f = MapAt[Expand,
		(f /. Solve[z D[df,z]+a df-f == 0,Table[w[i],{i,2,m}]])[[1]],
		{3}];
	ComposeSeries[f /. {f0 -> head[a,aa],
			    f1 -> Derivative[0,1][head][a,aa]},
			fz]
	];

hyper0F1[a_, fz_SeriesData, regflag_] :=
    Module[{i, m, f0 = faas[fz], tmp, z=fz[[1]]},
	If[intQ[a] && a <= 0,
	    If[regflag,
		m = 1-Round[a];
		fz^m hyper0F1[m+1,fz,True],
	      (* else *)
		ComplexInfinity
		],
	  (* else *)
	    If[Head[f0] === DirectedInfinity,
		tmp = 2 Sqrt[fz];
		If[regflag,
		    Hypergeometric1F1Regularized[a-1/2,2a-1,2tmp]*
				4^(a-1) Gamma[a-1/2]/Sqrt[Pi],
		  (* else *)
		    Hypergeometric1F1[a-1/2,2a-1,2tmp]
		    ] Exp[-tmp],
	      (* else *)
		m = mterms[1, fz];
		If[TrueQ[f0==0],
		    ComposeSeries[
			Sum[z^i/(i! regden[a,i,regflag]), {i,0,m}] + O[z]^(m+1),
			fz],
		  (* else *)
		    diffeq0F1[a,fz,f0,z,m, If[regflag,
			Hypergeometric0F1Regularized, Hypergeometric0F1]]
		    ]
		]
	    ]
	];

Hypergeometric0F1 /: Series[Hypergeometric0F1[a_,fz_],{z_,aa_,nn_Integer}] :=
    hyper0F1[a, Series[fz,{z,aa,nn}], False (* regularized *)] /;
	(nn >= 0 && FreeQ[a,z])

Hypergeometric0F1 /: Hypergeometric0F1[a_,s_SeriesData] :=
	hyper0F1[a, s, False (* regularized *)] /;
	(simser[s] && FreeQ[a,s[[1]]])

Hypergeometric0F1Regularized /:
Series[Hypergeometric0F1Regularized[a_,fz_],{z_,aa_,nn_Integer}] :=
	hyper0F1[a, Series[fz,{z,aa,nn}], True (* regularized *)] /;
	(nn >= 0 && FreeQ[a,z])

Hypergeometric0F1Regularized /:
Literal[Hypergeometric0F1Regularized[a_, s_SeriesData]] :=
    hyper0F1[a, s, True (* regularized *)] /; (simser[s] && FreeQ[a,s[[1]]])

 (* --- Hypergeometric 1F1 --- *)

	(* special cases: expansion about Infinity *)

Hypergeometric1F1 /:
Literal[ Series[Hypergeometric1F1[a_, b_, z_], {z_, Infinity, nn_Integer}]] :=
    Series[Hypergeometric1F1Regularized[a, b, z], {z, Infinity, nn}] Gamma[b] /;
	(nn >= 0 && FreeQ[a, z] && FreeQ[b, z])

Hypergeometric1F1Regularized /:
Series[Hypergeometric1F1Regularized[a_, b_, z_], {z_, Infinity, nn_Integer}] :=
    Module[{s1, s2, onSesss},
	(* AS 13.5.1 *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s1 = Series[hyper2F0r[a, 1+a-b, -z, nn], {z, Infinity, nn}];
	s2 = Series[hyper2F0r[b-a, 1-a, z, nn], {z, Infinity, nn}];
	s1 = E^(I Pi a) z^(-a) s1/Gamma[b-a] + E^z z^(a-b) s2/Gamma[a];
	If[onSesss, On[Series::esss]];
	s1
	] /; (nn >= 0 && FreeQ[a, z] && FreeQ[b, z])

	(* generic case *)

hyper1F1[a_,b_,fz_,regflag_,Uflag_] :=
    Module[{i, n = 1+a-b, m = mterms[1,fz], f0 = faas[fz], z = fz[[1]]},
	If[Head[f0] === DirectedInfinity,
	    If[Uflag,
		(* Abramowitz & Stegun 13.5.2 *)
		fz^(-a)Sum[Pochhammer[a,i]Pochhammer[n,i]/(i! (-fz)^i),{i,0,m}],
	      (* else *)
		(* Abramowitz & Stegun 13.5.1 *)
		(Exp[I Pi a] HypergeometricU[a,b,fz]/Gamma[b-a] +
				Exp[fz]/Gamma[a] HypergeometricU[b-a,b,-fz]) *
				If[regflag, 1, Gamma[b]]
		],
	  (* else *)
	    If[Uflag && TrueQ[f0==0],
		(* Abramowitz & Stegun 13.1.3 *)
		Pi/Sin[Pi b] (
		    diffeq1F1[a,b,fz,f0,z,m,True,False]/Gamma[n] -
		    diffeq1F1[n,2-b,fz,f0,z,m,True,False]/(Gamma[a] fz^(b-1))
		    ),
	      (* else *)
		diffeq1F1[a,b,fz,f0,z,m,regflag,Uflag]
		]
	    ]
	];

diffeq1F1[a_,b_,fz_,aa_,z_,m_,regflag_,Uflag_] :=
    Module[{i, f0, f1, df, f},
	If[TrueQ[aa==0] && !Uflag,
		f = Sum[Pochhammer[a,i] z^i/(i! regden[b,i,regflag]), {i,0,m}];
		Return[ComposeSeries[f + O[z]^(m+1), fz]]
		];
	f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i, {i,2,m}],{z,aa,m}];
	df = D[f,z];
	f = MapAt[Expand, (f /. Solve[z D[df,z] + (b-z) df - a f == 0,
					Table[w[i],{i,2,m}]])[[1]], {3}];
	If[regflag,
		f = f /. {f0 -> Hypergeometric1F1Regularized[a,b,aa],
		    f1 -> a Hypergeometric1F1Regularized[a+1,b+1,aa]},
	    (* else *)
		f = f /. If[Uflag,
			{f0 -> HypergeometricU[a,b,aa],
			f1 -> -a HypergeometricU[a+1,b+1,aa]},
		    (* else *)
			{f0 -> Hypergeometric1F1[a,b,aa],
			f1 -> a Hypergeometric1F1[a+1,b+1,aa]/b}]
		];
	ComposeSeries[f, fz]
	];

Hypergeometric1F1 /:
 Literal[ Series[Hypergeometric1F1[a_,b_,fz_],{z_,aa_,nn_Integer}]] :=
	Hypergeometric1F1[a,b,Series[fz,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[a,z] && FreeQ[b,z])

Hypergeometric1F1 /:
Literal[Hypergeometric1F1[a_,b_,s_SeriesData]] :=
	hyper1F1[a, b, s, False (* reg *), False (* U *)] /;
	(simser[s] && FreeQ[a,s[[1]]] && FreeQ[b,s[[1]]])

Hypergeometric1F1Regularized /:
Series[Hypergeometric1F1Regularized[a_,b_,fz_],{z_,aa_,nn_Integer}] :=
	Hypergeometric1F1Regularized[a,b,Series[fz,{z,aa,nn}]] /;
	(nn >= 0 && FreeQ[a,z] && FreeQ[b,z])

Hypergeometric1F1Regularized /:
Literal[Hypergeometric1F1Regularized[a_,b_,s_SeriesData]] :=
    Module[{m=1-Round[b]},
	Pochhammer[a,m] s^m * hyper1F1[a+m,m+1,s,True (* reg *),False (* U *)]
	]/; (simser[s] && FreeQ[a,s[[1]]] && intQ[b] && b <= 0)

Hypergeometric1F1Regularized /:
Literal[Hypergeometric1F1Regularized[a_,b_,s_SeriesData]] :=
	hyper1F1[a, b, s, True (* reg *), False (* U *)] /;
	(simser[s] && FreeQ[a,s[[1]]] && FreeQ[b,s[[1]]])

 (* --- Hypergeometric 2F1 --- *)

hyper2F1[a_,b_,c_,fz_,regflag_] :=
    Module[{m=b-a, ab, l, f0 = faas[fz], nn = mterms[1,fz], z = fz[[1]]},
	Which[
	    Head[f0] === DirectedInfinity,
		If[intQ[m],
		    m = Round[m];
		    ab = If[m > 0, b, a];
		    m = Abs[m];
		    l = c - a - m;
		    If[intQ[l], l = Round[l]];
		    If[IntegerQ[l] && l>=0,
			(* Erdelyi 2.1.4(19) p63 *)
			erdelyi21419[ab,m,l,fz,z,nn],
		      (* else *)
			(* Abramowitz & Stegun 15.3.14 *)
			as15314[ab,m,c,fz,z,nn]
			],
		  (* else *)
		    (* Abramowitz & Stegun 15.3.7 *)
		    as1537[a,b,c,fz,z]
		    ] * If[regflag, 1, Gamma[c]],
	    TrueQ[f0==1],
		m=c-a-b;
		If[intQ[m],
		    m = Round[m];
		    If[m >= 0,
			(* Abramowitz & Stegun 15.3.11 *)
			as15311[a,b,m,fz,z,nn],
		      (* else *)
			(* Abramowitz & Stegun 15.3.12 *)
			as15312[a,b,-m,fz,z,nn]
			],
		  (* else *)
		    (* Abramowitz & Stegun 15.3.6 *)
		    as1536[a,b,c,fz,z]
		    ] * If[regflag, 1, Gamma[c]],
	    TrueQ[f0==0],
		ComposeSeries[Sum[Pochhammer[a,l] Pochhammer[b,l]
		    z^l/(l! regden[c,l,regflag]), {l,0,nn}] + O[z]^(nn+1),fz],
	    True,
		diffeq2F1[a,b,c,fz,f0,z,nn,regflag]
	    ]
	];

diffeq2F1[a_,b_,c_,fz_,aa_,z_,m_,regflag_] :=
    Module[{i, f0, f1, df, f},
	f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i, {i,2,m}],{z,aa,m}];
	df = D[f,z];
	f = MapAt[Expand, (f /. Solve[z(1-z)D[df,z] + (c-(a+b+1)z) df
			- a b f == 0, Table[w[i],{i,2,m}]])[[1]], {3}];
	If[regflag,
	    f = f /. {f0 -> Hypergeometric2F1Regularized[a,b,c,aa],
		    f1 -> a b Hypergeometric2F1Regularized[a+1,b+1,c+1,aa]},
	  (* else *)
	    f = f /. {f0 -> Hypergeometric2F1[a,b,c,aa],
			f1 -> a b Hypergeometric2F1[a+1,b+1,c+1,aa]/c}
	    ];
	ComposeSeries[f, fz]
	]

as1537[a_,b_,c_,fz_,z_] :=
	Gamma[b-a]/(Gamma[b] Gamma[c-a]) (-fz)^-a *
			Hypergeometric2F1[a,a+1-c,a+1-b,1/fz] +
		Gamma[a-b]/(Gamma[a] Gamma[c-b]) (-fz)^-b *
			Hypergeometric2F1[b,b+1-c,b+1-a,1/fz];

as15314[a_,m_,c_,fz_,z_,nn_] :=
    Module[{n},
	1/Gamma[a+m] *
	(1/Gamma[c-a] (-fz)^(-a-m) *
		Sum[Pochhammer[a,n+m] Pochhammer[1-c+a,n+m]/(n! (n+m)! fz^n)
			((Log[-fz] /. logrule)+PolyGamma[1+m+n]+PolyGamma[1+n]-
			PolyGamma[a+m+n]-PolyGamma[c-a-m-n]),{n,0,nn-Re[a]-m}]+
		(-fz)^(-a) Sum[Gamma[m-n] Pochhammer[a,n]/
			(n! Gamma[c-a-n] fz^n), {n,0,Min[m-1,nn-Re[a]]}])
		];

erdelyi21419[a_,m_,l_,fz_,z_,nn_] :=
    Module[{n},
	1/Gamma[a+m] * ((-fz)^(-a)((-fz)^(-m) *
		((-1)^(m+l) Sum[Pochhammer[a,n+m](n-l)!/
				(n! (n+m)! fz^n),{n,l,nn-Re[a]-m}] +
		1/(l+m-1)! Sum[Pochhammer[n+m] Pochhammer[1-m-l,n+m]/
			(n! (n+m)! fz^n)((Log[-fz]/.logrule)+PolyGamma[1+m+n]+
			PolyGamma[1+n] - PolyGamma[a+m+n] - PolyGamma[l-n]),
			{n,0,Min[l-1,nn-Re[a]-m]}]) +
		Sum[(m-n-1)! Pochhammer[a,n]/(n! (m+l-n-1)! fz^n),
			{n,0,Min[m-1,nn-Re[a]]}]))
		];

as1536[a_,b_,c_,fz_,z_] :=
	Gamma[c-a-b]/(Gamma[c-a] Gamma[c-b]) *
			Hypergeometric2F1[a,b,a+b+1-c,1-fz] + 
		Gamma[a+b-c]/(Gamma[a] Gamma[b]) (1-fz)^(c-a-b) *
			Hypergeometric2F1[c-a,c-b,c+1-a-b,1-fz];

as15311[a_,b_,m_,fz_,z_,nn_] :=
    Module[{n,hyper=0},
	If[m!=0, hyper = Gamma[m]/(Gamma[a+m] Gamma[b+m]) *
			Sum[Pochhammer[a,n] Pochhammer[b,n] (1-fz)^n/
			(n! Pochhammer[1-m,n]), {n,0,Min[m-1,nn]}]];
	hyper - (fz-1)^m/(Gamma[a] Gamma[b]) *
			Sum[Pochhammer[a+m,n] Pochhammer[b+m,n]/(n! (n+m)!)
				(1-fz)^n ((Log[1-fz] /. logrule)-PolyGamma[n+1]-
				PolyGamma[n+m+1] + PolyGamma[a+n+m] +
				PolyGamma[b+n+m]), {n,0,nn}]
		];

as15312[a_,b_,m_,fz_,z_,nn_] :=
    Module[{n,hyper=0},
	If[m!=0, hyper = Gamma[m] (1-fz)^(-m)/(Gamma[a] Gamma[b]) *
			Sum[Pochhammer[a-m,n] Pochhammer[b-m,n] (1-fz)^n/
				(n! Pochhammer[1-m,n]), {n,0,Min[m-1,nn]}]];
	hyper - (-1)^m/(Gamma[a-m] Gamma[b-m]) *
			Sum[Pochhammer[a,n] Pochhammer[b,n]/(n! (n+m)!)
				(1-fz)^n ((Log[1-fz]/.logrule)-PolyGamma[n+1]-
				PolyGamma[n+m+1] + PolyGamma[a+n] +
				PolyGamma[b+n]), {n,0,nn}]
		];

Hypergeometric2F1 /:
Literal[ Series[Hypergeometric2F1[a_,b_,c_,fz_],{z_,aa_,nn_Integer}]] :=
	Hypergeometric2F1[a,b,c,Series[fz,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[a,z] && FreeQ[b,z] && FreeQ[c,z])

Hypergeometric2F1 /:
Literal[Hypergeometric2F1[a_,b_,c_,fz_SeriesData]] :=
    hyper2F1[a,b,c,fz,False] /;
    (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]] && FreeQ[c,fz[[1]]])

Hypergeometric2F1Regularized /:
Series[Hypergeometric2F1Regularized[a_,b_,c_,fz_],{z_,aa_,nn_Integer}] :=
	Hypergeometric2F1Regularized[a,b,c,Series[fz,{z,aa,nn}]] /;
	(nn >= 0 && FreeQ[a,z] && FreeQ[b,z] && FreeQ[c,z])

Hypergeometric2F1Regularized /:
Literal[Hypergeometric2F1Regularized[a_,b_,c_,fz_SeriesData]] :=
    Module[{m=1-Round[c]},
	Pochhammer[a,m] Pochhammer[b,m] fz^m hyper2F1[a+m,b+m,m+1,fz,True]
	]/; (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]] &&
		intQ[c] && c <= 0)

Hypergeometric2F1Regularized /:
Literal[Hypergeometric2F1Regularized[a_,b_,c_,fz_SeriesData]] :=
    hyper2F1[a,b,c,fz,True] /;
    (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]] && FreeQ[c,fz[[1]]])

 (* --- Hypergeometric U --- *)

	(* special case: expansion about Infinity *)

HypergeometricU /:
Series[HypergeometricU[a_,b_,z_],{z_, Infinity, nn_Integer}] :=
    Module[{n, onSesss}, (* AS 13.5.2 *)
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	n = z^(-a) Series[hyper2F0r[a, 1+a-b, -z, nn], {z, Infinity, nn}];
	If[onSesss, On[Series::esss]];
	n
	] /; (nn >= 0 && FreeQ[a, z] && FreeQ[b, z])

	(* generic cases *)

as1316[a_,b_,fz_,z_,nn_] :=
    Module[{k, n = Round[b]-1},
	(n-1)!/(Gamma[a] fz^n) *
		Sum[Pochhammer[a-n,k] fz^k/(k! Pochhammer[1-n,k]),
			{k,0,Min[n-1,nn+n]}] -
		(-1)^n/Gamma[a-n] *
		(Hypergeometric1F1Regularized[a,n+1,fz] (Log[fz]/.logrule) +
			Sum[Pochhammer[a,k] fz^k/(k! (n+k)!) (PolyGamma[a+k] -
				PolyGamma[1+k]-PolyGamma[1+n+k]),{k,0,nn}])
		];

HypergeometricU /: Series[HypergeometricU[a_,b_,fz_],{z_,aa_,nn_Integer}] :=
	HypergeometricU[a,b,Series[fz,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[a,z] && FreeQ[b,z])

HypergeometricU /:
Literal[HypergeometricU[a_, b_, fz_SeriesData]] :=
    Module[{k, n = a+1-b, nn = mterms[1,fz], f0 = faas[fz], z = fz[[1]]},
	If[intQ[b] && (f0 == 0),
		If[b < 0,
			(* Abramowitz & Stegun 13.1.29 *)
			fz^(1-b) as1316[n,2-b,fz,z,nn],
		    (* else *)
			(* Abramowitz & Stegun 13.1.6 *)
			as1316[a,b,fz,z,nn]
			],
	    (* else *)
		hyper1F1[a, b, fz, False (* reg *), True (* U *)]
		] /; (* exclude simple closed forms *)
			!((intQ[a] && a < .5) || (intQ[n] && n < .5)) &&
			esssQ[f0, HypergeometricU[a, b, fz]]
	]/; (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]])

 (* -------------------- Laguerre polynomials ------------------------ *)

LaguerreL /: Series[LaguerreL[n_,a_:0,fz_],{z_,aa_,nn_Integer}] :=
	LaguerreL[n,a,Series[fz,{z,aa,nn}]] /;
		(nn >= 0 && FreeQ[n,z] && FreeQ[a,z])

LaguerreL /: LaguerreL[n_,a_:0,s_SeriesData] :=
	Pochhammer[n+1,a] Hypergeometric1F1Regularized[-n,a+1,s] /;
	(simser[s] && FreeQ[n,s[[1]]] && FreeQ[a,s[[1]]])

 (* -------------- Legendre polynomials and functions ------------------ *)

legendreP[n_,m_,fz_,opt___Rule] :=
	If[SameQ[LegendreType /. {opt} /. Options[LegendreP], Real],
		(* Gradshteyn & Ryzhik 8.704(b) *)
		((1+fz)/(1-fz))^(m/2) *
			Hypergeometric2F1Regularized[-n,n+1,1-m,(1-fz)/2],
	    (* else *)
		(* Gradshteyn & Ryzhik 8.702 *)
		((fz+1)/(fz-1))^(m/2) *
			Hypergeometric2F1Regularized[-n,n+1,1-m,(1-fz)/2]
		];

legendreQ[n_,m_,fz_,opt___Rule] :=
	If[SameQ[LegendreType /. {opt} /. Options[LegendreQ], Real],
		If[TrueQ[m==Round[m]],
			(* Gradshteyn & Ryzhik 8.737.2, Erdelyi 3.4(14) *)
			Pi/(2 Sin[(m+n)Pi]) * (
			   Cos[(m+n)Pi] legendreP[n,m,fz,LegendreType->Real] -
				legendreP[n,m,-fz,LegendreType->Real]),
		    (* else *)
			(* Gradshteyn & Ryzhik 8.705(b), Erdelyi 3.4(17) *)
			Pi/(2 Sin[m Pi]) * (
			    Cos[m Pi] legendreP[n,m,fz,LegendreType->Real] -
				Gamma[m+n+1]/Gamma[n-m+1]
					legendreP[n,-m,fz,LegendreType->Real])
			],
	    (* else *)
		(* AS 8.1.3, Gradshteyn & Ryzhik 8.703 *)
		Exp[m Pi I] Gamma[m+n+1] Gamma[1/2] 2^-(n+1) (fz^2-1)^(m/2) *
		Hypergeometric2F1Regularized[(m+n+2)/2,(m+n+1)/2,n+3/2,fz^-2]/
		fz^(m+n+1)
		];

LegendreP /:
Series[Literal[LegendreP[n_,m_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
	LegendreP[n,m,Series[fz,{z,aa,nn}],opt] /;
		(nn >= 0 && FreeQ[n,z] && FreeQ[m,z])

LegendreP /: LegendreP[n_,m_,s_SeriesData,opt___Rule] :=
    legendreP[n,m,s,opt] /; (simser[s] && FreeQ[n,s[[1]]] && FreeQ[m,s[[1]]])

LegendreQ /:
Series[Literal[LegendreQ[n_,m_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
	LegendreQ[n,m,Series[fz,{z,aa,nn}],opt] /;
		(nn >= 0 && FreeQ[n,z] && FreeQ[m,z])

LegendreQ /: LegendreQ[n_,m_,s_SeriesData,opt___Rule] :=
    legendreQ[n,m,s,opt] /; (simser[s] && FreeQ[n,s[[1]]] && FreeQ[m,s[[1]]])

(* the following Legendre definitions are only needed because the pattern
matcher doesn't work right. the "m"'s above need to be ":0"'ed *)

LegendreP /:
Series[Literal[LegendreP[n_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
	LegendreP[n,0,Series[fz,{z,aa,nn}],opt] /;
	(nn >= 0 && FreeQ[n,z])

LegendreP /: LegendreP[n_,s_SeriesData,opt___Rule] :=
	legendreP[n,0,s,opt] /; (simser[s] && FreeQ[n,s[[1]]])

LegendreQ /:
Series[Literal[LegendreQ[n_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
	LegendreQ[n,0,Series[fz,{z,aa,nn}],opt] /;
	(nn >= 0 && FreeQ[n,z])

LegendreQ /: LegendreQ[n_,s_SeriesData,opt___Rule] :=
	legendreQ[n,0,s,opt] /; (simser[s] && FreeQ[n,s[[1]]])

 (* ----------------------- Lerch's transcendant ---------------------- *)

LerchPhi /:
Series[Literal[LerchPhi[fz_,s_,a_,opt___Rule]],{z_,aa_,nn_Integer}] :=
	LerchPhi[Series[fz,{z,aa,nn}],s,a,opt] /;
	(nn >= 0 && FreeQ[s,z] && FreeQ[a,z])

LerchPhi /:
Literal[LerchPhi[fz_SeriesData,s_,a_,opt___Rule]] :=
    Module[{k},
	Sum[If[TrueQ[a+k==0], 0, fz^k/((a+k)^2)^(s/2)],
	    {k,0,mterms[1,fz]}] /;
	!(DoublyInfinite /. {opt} /. Options[LerchPhi]) && (faas[fz] == 0)
	] /; (simser[fz] && FreeQ[s,fz[[1]]] && FreeQ[a,fz[[1]]])

 (* ------------------ Logarithmic integral ----------------------- *)

	(* special case: expansion about Infinity *)

LogIntegral /: Series[LogIntegral[z_], {z_, Infinity, nn_Integer}] :=
    Module[{j, s, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s = Series[Sum[j! z^(-j-1), {j, 0, nn}], {z, Infinity, nn}];
	s = z (s /. z -> Log[z]);
	If[onSesss, On[Series::esss]];
	s
	] /; nn >= 0

	(* generic case *)

LogIntegral /: Series[LogIntegral[fz_],{z_,aa_,nn_Integer}] :=
	LogIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0

LogIntegral /:
Literal[LogIntegral[s_SeriesData]] :=
    Module[{f0 =faas[s], i},
	ExpIntegralEi[Log[s]] /;
	    !TrueQ[f0 == 0] && esssQ[f0, LogIntegral[s]]
	] /; simser[s]

 (* ------------------------- Pochhammer --------------------------------- *)

Pochhammer/:
Series[Pochhammer[f1z_,f2z_],{z_,aa_,nn_Integer}] :=
    Module[{sf1 = Series[f1z,{z,aa,nn}],sf2 = Series[f2z,{z,aa,nn}]},
	Gamma[sf1+sf2] / Gamma[sf1]] /; nn >= 0

Pochhammer/:
Literal[Pochhammer[sf1_SeriesData, sf2_SeriesData]]:=
	Gamma[sf1+sf2] / Gamma[sf1] /; simser[sf1] && simser[sf1]

Pochhammer/:
Literal[Pochhammer[sf_SeriesData,a_]]:=
    Gamma[sf +
	Series[a,{sf[[1]],sf[[2]],Round[(sf[[5]]-Min[0,sf[[4]]])/sf[[6]]]}]] /
    Gamma[sf] /; simser[sf]

Pochhammer/:
Literal[Pochhammer[a_,sf_SeriesData]] :=
    Module[{sa},
	sa = Series[a, {sf[[1]], sf[[2]],
		Round[(sf[[5]] - Min[0, sf[[4]]])/sf[[6]]]}];
	Gamma[sf+sa] / Gamma[sa]] /; simser[sf]

 (* -------------------- PolyGamma functions -------------------------- *)

	(* special cases: expansion about Infinity *)

LogGamma /: Series[LogGamma[z_],{z_, Infinity, nn_Integer}] :=
    Module[{s, i, onSesss},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s = Sum[BernoulliB[i]/(i(i-1) z^(i-1)), {i, 2, nn, 2}];
	s = (z - 1/2) Log[z] - z + Log[2 Pi]/2 + Series[s, {z, Infinity, nn}];
	If[onSesss, On[Series::esss]];
	s
	] /; nn >= 0

PolyGamma /: Series[PolyGamma[n_Integer:0, z_],{z_, Infinity, nn_Integer}] :=
    Module[{s, i},
	onSesss = OnQ[Series::esss];
	If[onSesss, Off[Series::esss]];
	s = Sum[BernoulliB[i]/(i z^i), {i, 2, nn, 2}];
	s = Log[z] - 1/(2z) - Series[s, {z, Infinity, nn}];
	s = D[s, {z, n}];
	If[onSesss, On[Series::esss]];
	s
	] /; n >= 0 && nn >= 0

	(* generic cases *)

LogGamma /: Series[LogGamma[fz_],{z_,aa_,nn_Integer}] :=
	LogGamma[Series[fz,{z,aa,nn}]] /; nn >= 0

LogGamma /: LogGamma[s_SeriesData] :=
    Module[{f0 = faas[s]},
	(f0 = If[intQ[f0] && f0 < .5,
	    (Round[f0]-1) Pi I,
	    LogGamma[f0] - Log[Gamma[f0]]
	    ];
	Log[Gamma[s]] + f0) /; esssQ[f0, LogGamma[s]]
	] /; simser[s]

PolyGamma /: Series[PolyGamma[n_Integer:0,fz_],{z_,aa_,nn_Integer}] :=
	PolyGamma[n,Series[fz,{z,aa,nn+n}]] /; nn >= n

PolyGamma /: PolyGamma[-1,s_SeriesData] :=
	LogGamma[s] /; simser[s] && esssQ[faas[s], LogGamma[s]]

PolyGamma /: PolyGamma[n_Integer:0, s_SeriesData] :=
    Module[{a, pgs, f0 = faas[s], z = s[[1]]},
	(pgs = Gamma[Series[z, {z, f0, mterms[1, s]}]];
	pgs = D[D[pgs, z]/pgs, {z, n}];
	pgs = ComposeSeries[pgs, s] /. (a_ :> (Expand[a] /; FreeQ[a,z]))) /;
	Head[f0] =!= DirectedInfinity
	] /; n >= 0 && simser[s]

 (* ------------------- PolyLogarithm functions ------------------------- *)

PolyLog /: Series[PolyLog[n_,fz_],{z_,aa_,nn_Integer}] :=
	PolyLog[n,Series[fz,{z,aa,nn}]] /; (nn >= 0  && FreeQ[n,z])

polylog1[n_,fz_,m_] :=
    Module[{k, fpow=1, plder, logpart=0, t, sum=0, i},
	If[TrueQ[n==1], Return[-Log[1-fz] //. logrule]];
	plder = PolyLog[n,t];
	Do[sum += (plder /. t-> 1)fpow/i!;
		If[logpart =!= 0,
			logpart = logpart /. t -> 1-t;
			logpart = Normal[Series[logpart,{t,0,m-i}]];
			Do[logpart = -Integrate[logpart,t],{i}];
			sum += (logpart /. t -> 1-fz)
			];
		fpow *= fz-1;
		{logpart,plder} = logsplit[D[plder,t]],
		{i,0,m}];
	sum //. logrule
	]

PolyLog /: Literal[PolyLog[n_,fz_SeriesData]] :=
    Module[{k, f0=faas[fz], m = mterms[1, fz]},
	Which[TrueQ[f0==0], Sum[fz^k/k^n,{k,m}],
	    TrueQ[f0==1] && IntegerQ[n] && n>0, polylog1[n,fz,m],
	    True, Sum[Derivative[0,k][PolyLog][n,f0](fz-f0)^k/(k!),{k,0,m}]
	    ]
	] /; simser[fz] && FreeQ[n, fz[[1]]]

 (* ------------------- Spherical Harmonic functions --------------------- *)

SphericalHarmonicY /:
Series[SphericalHarmonicY[l_,m_,ftheta_,phi_],{theta_,aa_,nn_Integer}] :=
	SphericalHarmonicY[l,m,Series[ftheta,{theta,aa,nn}],phi] /;
	(nn >= 0 && FreeQ[l,theta] && FreeQ[m,theta] && FreeQ[phi,theta])

SphericalHarmonicY /:
Literal[SphericalHarmonicY[l_,m_,s_SeriesData,phi_]] :=
	Exp[m I phi] Sqrt[(2l+1)/(4 Pi Pochhammer[l-m+1,2m])] *
		LegendreP[l,m,Cos[s],LegendreType->Real] /;
	(simser[s] && FreeQ[l,s[[1]]] && FreeQ[m,s[[1]]] && FreeQ[phi,s[[1]]])

 (* ----------------------- Zeta function ---------------------------- *)

$stltjesgammas;
$NStltjesGamma;
$PrecNStltjesGamma;

(* values for debugging:
$NStltjesGamma[1] = -0.0728158454836767248605863758749013191377363383343
$NStltjesGamma[2] = -0.00969036319287231848453038604
$NStltjesGamma[3] = 0.0020538344203033458661600465
$NStltjesGamma[4] = 0.0023253700654673000574681702
$NStltjesGamma[5] = 0.000793323817301062701753335
$NStltjesGamma[6] = -0.000238769345430199609872422
*)

StieltjesGamma[0] = EulerGamma;

N[StieltjesGamma[n_Integer]] ^:= N[StieltjesGamma[n], $MachinePrecision]

N[StieltjesGamma[n_Integer], m_Integer] ^:=
    Module[{eglist, i},
	If[NumberQ[$NStltjesGamma[n]] && m <= $PrecNStltjesGamma[n],
	    N[$NStltjesGamma[n], m],
	  (* else *)
	    eglist = EGList[n, m];
	    AbortProtect[
		Do[ If[!NumberQ[$NStltjesGamma[i]] || m > $PrecNStltjesGamma[i],
		    {$NStltjesGamma[i], $PrecNStltjesGamma[i]} =
			{eglist[[i]], m}], {i, n}]];
	    $NStltjesGamma[n]]] /; (n > 0 && m > 0)

EGList[n_, m_] :=
    Module[{zser, i},
	zser = ($stltjesgammas - 1) N[EulerGamma, m] -
		Sum[$NSigmaZeta[i, m] (1 - $stltjesgammas)^i, {i, n+1, 2, -1}];
	zser = Exp[Series[zser, {$stltjesgammas, 1, n+1}]];
	zser = Drop[zser[[3]], 2];
	zser Table[(-1)^i i!, {i, n}]];

$NSigmaZeta;
$PrecNSigmaZeta;

(* values for debugging:
$NSigmaZeta[2] = 0.0937731164201826122986016923027207941919722317905
$NSigmaZeta[3] = 0.017229544011064297934002741028
$NSigmaZeta[4] = 0.00368791470636343601614505920
$NSigmaZeta[5] = 0.000904895577699075748249223220
$NSigmaZeta[6] = 0.000241132534087530523369496737
$NSigmaZeta[7] = 0.000067363439740772150048502904
*)

$NSigmaZeta[n_, m_] :=
    Module[{szlist, i},
	(* coefficients in the series expansion of Log[(s-1) Zeta[s]]
		about s == 1 *)
	If[NumberQ[$NSigmaZeta[n]] && m <= $PrecNSigmaZeta[n],
	    Return[N[$NSigmaZeta[n], m]]];
	szlist = SZList[n, m];
	AbortProtect[
	    Do[ If[!NumberQ[$NSigmaZeta[i]] || m > $PrecNSigmaZeta[i],
		{$NSigmaZeta[i], $PrecNSigmaZeta[i]} = {szlist[[i-1]], m}],
		{i, 2, n}]];
	$NSigmaZeta[n]]

$Nzeta32;
$PrecNzeta32;

(* values for debugging:
$Nzeta32[2] = 0.1168502750680849136771556874922594459571062129525
$Nzeta32[3] = 0.017266596754881666574923630441
$Nzeta32[4] = 0.00366950790104801363656336638
$Nzeta32[5] = 0.000904752559027923226702063001
$Nzeta32[6] = 0.000241179440157020317746431190
$Nzeta32[7] = 0.000067364093196650679301661642
*)

$Nzeta32[i_, m_] :=
    Module[{zeta},
	If[NumberQ[$Nzeta32[j]] && m <= $PrecNzeta32[j],
	    Return[N[$Nzeta32[j], m]]];
	zeta = N[Zeta[i, 3/2], m]/(i 2^i);
	AbortProtect[{$Nzeta32[i], $PrecNzeta32[i]} = {zeta, m}];
	zeta];

SZList[n_, m_] :=
    Module[{xiser, i},
	(* series expansion of 2*xi(s) about s == 1 *)
	xiser = 1 + Sum[($NXiBeta[i-2, m] + $NXiBeta[i-1, m]) *
		($stltjesgammas - 1)^i, {i, n}];
	(* series expansion of log(2*xi(s)) about s == 1 *)
	xiser = Log[Series[xiser, {$stltjesgammas, 1, n}]];
	xiser = Rest[xiser[[3]]] Table[(-1)^i, {i, n-1}];
	(* coefficients in the series expansion of Log[(s-1) Zeta[s]]
		about s == 1 *)
	Table[$Nzeta32[i, m], {i, 2, n}] + xiser];

$NXiBeta;
$PrecNXiBeta;

$NXiBeta[-1, _] := 0;
$NXiBeta[0, m_] := N[1 + EulerGamma/2 - Log[2 Sqrt[Pi]], m];

(* values for debugging:
$NXiBeta[1] = 0.000248155568105149320572108979315189423336753288730593
$NXiBeta[2] = 0.0002498282818177993517790627950742989298
$NXiBeta[3] = 3.3534484987276532770582890420745483*10^-6
$NXiBeta[4] = 1.69680629349152089252694618293863650*10^-6
$NXiBeta[5] = 2.418074837001468525320879096545763*10^-8
$NXiBeta[6] = 8.197666248796084350271868231276427*10^-9
*)

betatheta[t_, m_] :=
    Module[{sum, term, cterm, dterm, eps = 0.1^m},
	(* evaluate the theta function Sum[E^(-n^2 Pi t), {n, 1, Infinity}] *)
	sum = term = Exp[-N[Pi, m] t];
	cterm = term;
	dterm = cterm^2;
	While[term > sum eps,
		cterm *= dterm;
		term *= cterm;
		sum += term
		];
	sum
	];

$NXiBeta[j_, m_] :=
    Module[{prec = m+10, t, beta, onNIslwcon},
	If[NumberQ[$NXiBeta[j]] && m <= $PrecNXiBeta[j],
	    Return[N[$NXiBeta[j], m]]];
	onNIslwcon = OnQ[NIntegrate::slwcon];
	If[onNIslwcon, Off[NIntegrate::slwcon]];
	beta = NIntegrate[betatheta[t, prec] (Log[t]/2)^j (Sqrt[t] + (-1)^j)/t,
		{t, 1, Infinity}, Method -> DoubleExponential,
		WorkingPrecision -> prec, PrecisionGoal -> m,
		MaxRecursion -> 20]/j!;
	If[onNIslwcon, On[NIntegrate::slwcon]];
	AbortProtect[{$NXiBeta[j], $PrecNXiBeta[j]} = {beta, m}];
	beta];


Zeta /: Series[Zeta[fz_],{z_,aa_,nn_Integer}] :=
	Zeta[Series[fz,{z,aa,nn}]] /; nn >= 0

Zeta /: Literal[Zeta[s_SeriesData]] :=
    Module[{j, k, f, f0 = faas[s], m = mterms[1,s], z = s[[1]]},
	If[TrueQ[f0==1],
		f = 1/z+EulerGamma+Sum[(-1)^j StieltjesGamma[j] z^j/j!,{j,m}] +
			O[z]^(m+1),
	    (* else *)
		f = Sum[Derivative[k][Zeta][f0]z^k/k!, {k,0,m}] + O[z]^(m+1)
		];
	f[[2]] = f0;
	ComposeSeries[f, s]
	] /; simser[s]

 (* ------------------- Various Integral functions ----------------------- *)

SinIntegral /: Series[SinIntegral[fz_],{z_,aa_,nn_Integer}] :=
	SinIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0

SinIntegral /: Literal[SinIntegral[fz_SeriesData]] :=
   Module[{temp, f, f0 = faas[fz]},
	f = Integrate[Series[Sin[temp+f0]/(temp+f0),{temp,0,mterms[1,fz]}],
		temp];
	SinIntegral[f0] + ComposeSeries[f, fz-f0]
	] /; simser[fz]

CosIntegral /: Series[CosIntegral[fz_],{z_,aa_,nn_Integer}] :=
	CosIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0

CosIntegral /: Literal[CosIntegral[fz_SeriesData]] :=
    Module[{i, temp, f, m = mterms[1,fz], f0 = faas[fz]},
	If[TrueQ[f0 == 0],
		f = Sum[(-1)^i temp^(2i)/(2i(2i)!),{i,1,m/2}] + O[temp]^(m+1);
		ComposeSeries[f, fz] + (Log[fz] /. logrule) + EulerGamma,
	    (* else *)
		f = Integrate[Series[Cos[temp+f0]/(temp+f0),{temp,0,m}], temp];
		CosIntegral[f0] + ComposeSeries[f, fz-f0]
		]
	] /; simser[fz]

SinhIntegral /: Series[SinhIntegral[fz_],{z_,aa_,nn_Integer}] :=
	SinhIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0

SinhIntegral /:
Literal[SinhIntegral[fz_SeriesData]] :=
   Module[{temp, f, f0=faas[fz]},
	f = Integrate[Series[Sinh[temp+f0]/(temp+f0),{temp,0,mterms[1,fz]}],
		temp];
	SinhIntegral[f0] + ComposeSeries[f, fz-f0]
	] /; simser[fz]

CoshIntegral /: Series[CoshIntegral[fz_],{z_,aa_,nn_Integer}] :=
	CoshIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0

CoshIntegral /:
Literal[CoshIntegral[fz_SeriesData]] :=
    Module[{i, temp, f, m = mterms[1,fz], f0 = faas[fz]},
	If[TrueQ[f0 == 0],
		f = Sum[temp^(2i)/(2i(2i)!),{i,1,m/2}] + O[temp]^(m+1);
		ComposeSeries[f, fz] + (Log[fz] /. logrule) + EulerGamma,
	    (* else *)
		f = Integrate[Series[Cosh[temp+f0]/(temp+f0),{temp,0,m}], temp];
		CoshIntegral[f0] + ComposeSeries[f, fz-f0]
		]
	] /; simser[fz]

FresnelC/: Series[FresnelC[fz_],{z_,aa_,nn_Integer}] :=
	FresnelC[Series[fz,{z,aa,nn}]] /; nn >= 0

FresnelC/: Literal[FresnelC[fz_SeriesData]] :=
   Module[{temp, f, f0=faas[fz]},
	f = Integrate[Series[Cos[Pi (temp+f0)^2/2], {temp,0,mterms[1,fz]}],
		temp];
	FresnelC[f0] + ComposeSeries[f, fz-f0]
	] /; simser[fz]

FresnelS/: Series[FresnelS[fz_],{z_,aa_,nn_Integer}] :=
	FresnelS[Series[fz,{z,aa,nn}]] /; nn >= 0

FresnelS/: Literal[FresnelS[fz_SeriesData]] :=
   Module[{temp, f, f0=faas[fz]},
	f = Integrate[Series[Sin[Pi (temp+f0)^2/2], {temp,0,mterms[1,fz]}],
		temp];
	FresnelS[f0] + ComposeSeries[f, fz-f0]
	] /; simser[fz]

 (* ----------------------- Sums and Integrals ------------------------- *)

Series[Literal[Integrate[g_,{t_,f1z_,f2z_}]],{z_,aa_,nn_Integer}] :=
    Module[{sf1 = Series[f1z,{z,aa,nn}],sf2 = Series[f2z,{z,aa,nn}],
		f01,f02,f1,f2,temp},
	(If[Head[sf1] === SeriesData,
		f01 = faas[sf1];
		f1 = Integrate[Series[(g /. t -> temp+f01),
			{temp,0,mterms[1,sf1[[4]],sf1[[5]],sf1[[6]]]}],temp];
		f1 = ComposeSeries[f1,sf1-f01],
	    (* else *)
		f01 = sf1;
		f1 = 0
		];
	If[Head[sf2] === SeriesData,
		f02 = faas[sf2];
		f2 = Integrate[Series[(g /. t -> temp+f02),
			{temp,0,mterms[1,sf2[[4]],sf2[[5]],sf2[[6]]]}],temp];
		f2 = ComposeSeries[f2,sf2-f02],
	    (* else *)
		f02 = sf2;
		f2 = 0
		];
	Integrate[g,{t,f01,f02}] + ((f2-f1) /. temp -> z)) /;
		((Head[sf1] =!= Series) && (Head[sf2] =!= Series))
	] /; (nn >= 0 && FreeQ[g,z] && FreeQ[t,z])

Series[Literal[Sum[g_,iter_List]],{z_,aa_,nn_Integer}] :=
	Sum[Series[g,{z,aa,nn}],iter] /; (nn >= 0 && FreeQ[iter,z])

If[onSDsdatc, On[SeriesData::sdatc]];
Remove[onSDsdatc];

End[] (* "SpecialFunctions`Series`Private`" *)

SetAttributes[AiryAi, ReadProtected];
SetAttributes[AiryAiPrime, ReadProtected];
SetAttributes[AiryBi, ReadProtected];
SetAttributes[AiryBiPrime, ReadProtected];
SetAttributes[BesselI, ReadProtected];
SetAttributes[BesselJ, ReadProtected];
SetAttributes[BesselK, ReadProtected];
SetAttributes[BesselY, ReadProtected];
SetAttributes[Beta, ReadProtected];
SetAttributes[BetaRegularized, ReadProtected];
SetAttributes[ChebyshevT, ReadProtected];
SetAttributes[ChebyshevU, ReadProtected];
SetAttributes[ExpIntegralE, ReadProtected];
SetAttributes[ExpIntegralEi, ReadProtected];
SetAttributes[Factorial, ReadProtected];
SetAttributes[FresnelC, ReadProtected];
SetAttributes[FresnelS, ReadProtected];
SetAttributes[Gamma, ReadProtected];
SetAttributes[GammaRegularized, ReadProtected];
SetAttributes[GegenbauerC, ReadProtected];
SetAttributes[HermiteH, ReadProtected];
SetAttributes[Hypergeometric0F1, ReadProtected];
SetAttributes[Hypergeometric1F1, ReadProtected];
SetAttributes[Hypergeometric2F1, ReadProtected];
SetAttributes[Hypergeometric0F1Regularized, ReadProtected];
SetAttributes[Hypergeometric1F1Regularized, ReadProtected];
SetAttributes[Hypergeometric2F1Regularized, ReadProtected];
SetAttributes[HypergeometricU, ReadProtected];
SetAttributes[LaguerreL, ReadProtected];
SetAttributes[LegendreP, ReadProtected];
SetAttributes[LegendreQ, ReadProtected];
SetAttributes[LerchPhi, ReadProtected];
SetAttributes[LogIntegral, ReadProtected];
SetAttributes[LogGamma, ReadProtected];
SetAttributes[Pochhammer, ReadProtected];
SetAttributes[PolyGamma, ReadProtected];
SetAttributes[PolyLog, ReadProtected];
SetAttributes[SphericalHarmonicY, ReadProtected];
SetAttributes[Zeta, ReadProtected];
SetAttributes[SinIntegral, ReadProtected];
SetAttributes[CosIntegral, ReadProtected];
SetAttributes[SinhIntegral, ReadProtected];
SetAttributes[CoshIntegral, ReadProtected];
SetAttributes[SeriesData, ReadProtected];
SetAttributes[Series, ReadProtected];
SetAttributes[StieltjesGamma, ReadProtected];

Protect[AiryAi,AiryAiPrime,AiryBi,AiryBiPrime,BesselI,BesselJ,BesselK,BesselY,
Beta,BetaRegularized, ChebyshevT,ChebyshevU,ExpIntegralE,ExpIntegralEi,
Factorial,FresnelC,FresnelS,Gamma,GammaRegularized,GegenbauerC,HermiteH,
Hypergeometric0F1,Hypergeometric0F1Regularized,Hypergeometric1F1,
Hypergeometric1F1Regularized,Hypergeometric2F1,Hypergeometric2F1Regularized,
HypergeometricU,LaguerreL,LegendreP,LegendreQ,LerchPhi,LogIntegral,LogGamma,
Pochhammer,PolyGamma,PolyLog,SphericalHarmonicY,Zeta,SinIntegral,CosIntegral,
SinhIntegral,CoshIntegral,SeriesData,Series,StieltjesGamma];

End[]  (* "System`" *)

