home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e031 / 3.ddi / MATHZIP2 / STARTUP / SERIES.M < prev    next >
Encoding:
Text File  |  1991-09-23  |  65.7 KB  |  1,965 lines
  1.  
  2. (* :Name: System`Series` *)
  3.  
  4. (* :Title: InverseSeries and Series of Special Functions *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Defines InverseSeries[ ] using Newton's method symbolically.
  10.     Provides series of special functions using special properties
  11.     of those functions.  Also allows series expansions of special
  12.     functions about singularities of certain simple types.
  13. *)
  14.  
  15. (* :Context: System` *)
  16.  
  17. (* :Package Version: 2.0 *)
  18.  
  19. (* :Copyright: Copyright 1990  Wolfram Research, Inc.
  20.  
  21.     Permission is hereby granted to modify and/or make copies of
  22.     this file for any purpose other than direct profit, or as part
  23.     of a commercial product, provided this copyright notice is left
  24.     intact.  Sale, other than for the cost of media, is prohibited.
  25.  
  26.     Permission is hereby granted to reproduce part or all of
  27.     this file, provided that the source is acknowledged.
  28. *)
  29.  
  30. (* :History:
  31.     Updated by Jerry B. Keiper, December 1990.
  32.     Version 1.2 by Jerry B. Keiper, 1988.
  33.     InverseSeries[ ] was originally by Daniel R. Grayson.
  34. *)
  35.  
  36. (* :Keywords:
  37.     Series, InverseSeries, special functions.
  38. *)
  39.  
  40. (* :Source:
  41.     M. Abramowitz & I. Stegun, Handbook of Mathematical Functions With
  42.     Formulas, Graphs, and Mathematical Tables, National Bureau of
  43.     Standards, Applied Mathematics Series, Vol. 55.  Reprinted by
  44.     Dover.
  45. *)
  46.  
  47. (* :Mathematica Version: 2.0 *)
  48.  
  49. (* :Limitation:
  50.     InverseSeries[ ] will not find the inverse series of a series
  51.     expansion about a transcendental singularity.
  52. *)
  53.  
  54.  (* --------------------- reversion of series --------------------------- *)
  55.  
  56. Begin["System`"]
  57.  
  58. Unprotect[ InverseSeries, Series, SeriesData]
  59.  
  60. Map[ Clear, {InverseSeries, Series, SeriesData}]
  61.  
  62. InverseSeries::usage = "InverseSeries[s] takes the series s generated by
  63. Series, and gives a series for the inverse of the function represented by s.
  64. InverseSeries[s, y] gives the inverse series in terms of the variable y."
  65.  
  66. Begin["InverseSeries`Private`"]
  67.  
  68. InverseSeries[f_SeriesData] := InverseSeries[f, f[[1]]];
  69.  
  70. InverseSeries[f_SeriesData, y_Symbol] :=
  71.     Module[{s},
  72.         s /; (s = f; s[[2]] = 0; s = 1/InverseSeries[s, y];
  73.             Head[s] === SeriesData)
  74.     ] /; Head[f[[2]]] === DirectedInfinity
  75.  
  76. InverseSeries[f_SeriesData, y_Symbol] :=
  77.     Which[
  78.         f[[4]] === 0,
  79.         invser[f[[3,1]], 0, 1, f - f[[3,1]], y],
  80.         f[[4]] > 0,
  81.         invser[0, 0, 1, f, y],
  82.         f[[4]] < 0,
  83.         invser[DirectedInfinity[ ], 0, 1, 1/f, y]
  84.     ] /; (Length[f[[3]]] > If[ f[[4]] === 0, 1, 0] && IntegerQ[f[[4]]] &&
  85.         Head[f[[1]]] === Symbol && (FreeQ[f, y] || f[[1]] === y) &&
  86.         Head[f[[2]]] =!= DirectedInfinity)
  87.  
  88. InverseSeries[y0_. + (x_ b_. + mx0_.)^(a_) f_SeriesData, y_Symbol] :=
  89.     invser[y0,a,b,f,y] /;
  90.     (f[[1]] === x && Head[x] === Symbol && TrueQ[f[[2]] b == -mx0] &&
  91.     FreeQ[y0,x] && FreeQ[mx0,x] && FreeQ[a,x] && FreeQ[b,x] &&
  92.     (Length[f[[3]]] > 0) && (x === y ||
  93.         (FreeQ[y0, y] && FreeQ[mx0,y] && FreeQ[a,y] && FreeQ[b,y] &&
  94.             FreeQ[f,y])))
  95.  
  96. InverseSeries[y0_. + b_ (xmx0_)^(a_) f_SeriesData, y_Symbol] :=
  97.     Module[{s = b f}, InverseSeries[y0 + (xmx0)^a s, y]] /; FreeQ[b, x]
  98.  
  99. invser[yc_, a_, b_, f_, y_] :=
  100.     Module[{ba = f[[6]]/(a f[[6]]+f[[4]]), y0, z, s, w, newts},
  101.         y0 = If[Head[yc] === DirectedInfinity, 0, yc];
  102.         s = b z SeriesData[z,0,f[[3]]/f[[3,1]],0,f[[5]]-f[[4]],f[[6]]]^ba;
  103.         s[[3]] = Together[s[[3]]];
  104.         s = Series[Normal[s], {z, 0, Floor[(s[[5]]-1)/s[[6]]]}];
  105.         w = Series[(y-y0)/f[[3,1]], {y, y0, Length[s[[3]]]+1}]^ba;
  106.         newts = newton[s];
  107.         If[Head[w] === SeriesData, 
  108.         w = ComposeSeries[newts, w] + f[[2]];
  109.         If[Head[yc] === DirectedInfinity, w[[2]] = yc],
  110.           (* else *)
  111.         (* Another horrible kludge in Series[ ].  This one is because
  112.             ComposeSeries[ ] doesn't. *)
  113.         w = (newts + f[[2]]) /. z -> ((y-y0)/f[[3,1]])^ba;
  114.         ];
  115.         w] /; Length[f[[3]]] > 0
  116.  
  117. newton[f_SeriesData] :=
  118.     Module[{n=f[[5]], d=f[[6]], g, var=f[[1]], i, ff, gg, hh, kk},
  119.         g = SeriesData[var,0,{1/f[[3,1]]},d,n,d];
  120.         i = 2d;
  121.         While[ i < n,
  122.         i *= 2;
  123.         ff = f + SeriesData[var,0,{},i,i,d];
  124.         gg = g + SeriesData[var,0,{},i,i,d];
  125.         kk = (SeriesData[var,0,{1},d,i,d] - ComposeSeries[ff,gg])/
  126.                 ComposeSeries[D[ff,var],gg];
  127.         kk[[5]] = n;
  128.         g = g + kk
  129.         ];
  130.         g]
  131.  
  132. Protect[InverseSeries];
  133.  
  134. End[]  (* "InverseSeries`Private`" *)
  135.  
  136.  (* --------------------- Special Functions ---------------------------- *)
  137.  
  138. StieltjesGamma::usage =
  139. "StieltjesGamma[n] is the nth Stieltjes constant and is related to the
  140. Laurent series expansion of Zeta[s] about s == 1."
  141.  
  142. Unprotect[AiryAi,AiryAiPrime,AiryBi,AiryBiPrime,BesselI,BesselJ,BesselK,
  143. BesselY,Beta,BetaRegularized,ChebyshevT,ChebyshevU,ExpIntegralE,ExpIntegralEi,
  144. Factorial,FresnelC,FresnelS,Gamma,GammaRegularized,GegenbauerC,HermiteH,
  145. Hypergeometric0F1,Hypergeometric0F1Regularized,Hypergeometric1F1,
  146. Hypergeometric1F1Regularized,Hypergeometric2F1,Hypergeometric2F1Regularized,
  147. HypergeometricU,LaguerreL,LegendreP,LegendreQ,LerchPhi,LogIntegral,LogGamma,
  148. Pochhammer,PolyGamma,PolyLog,SphericalHarmonicY,Zeta,SinIntegral,CosIntegral,
  149. SinhIntegral,CoshIntegral,SeriesData,Series];
  150.  
  151. Begin["SpecialFunctions`Series`Private`"]
  152.  
  153.  
  154. OnQ[mess_] := MemberQ[{MessageName, String}, Head[mess]];
  155.  
  156. onSDsdatc = OnQ[SeriesData::sdatc];
  157. If[onSDsdatc, Off[SeriesData::sdatc]];
  158.  
  159.  (* -----------------  utilities for SeriesData  -------------------------- *)
  160.  
  161.     (* Evaluate the limit of a series if the limit point is
  162.      the same as the expansion point. *)
  163.  
  164. SeriesData/: Literal[Limit[s_SeriesData, x_ -> a_]] :=
  165.     faas[s] /; (simser[s] && s[[1]] === x && s[[2]] == a)
  166.  
  167.     (* Coefficient of a SeriesData object *)
  168.  
  169. SeriesData/: Literal[Coefficient[s_SeriesData, x___]] :=
  170.     Coefficient[Normal[s], x] /; simser[s]
  171.  
  172.     (* CoefficientList of a SeriesData object *)
  173.  
  174. SeriesData/: Literal[CoefficientList[s_SeriesData, x___]] :=
  175.     CoefficientList[Normal[s], x] /; simser[s]
  176.  
  177.  (* --------------------- powers of series ----------------------------- *)
  178.  
  179. (*
  180. Commented out by DJW
  181. SeriesData /: Power[s_SeriesData, v_?NumberQ] :=
  182.     Module[{vv},(s^vv) /. vv -> v]
  183. *)
  184.  
  185. SeriesData /:
  186. Power[SeriesData[z_, z0_, cl_List, nmin_, nmax_, den_], v_] :=
  187.     Module[{cla, a, na, p, exp = v nmin/den},
  188.     If[Length[cl] == 0, a = 1; cla = cl, a = cl[[1]]; cla = cl/a];
  189.     na = N[a];
  190.     If[NumberQ[na] && Re[na] < 0,
  191.         (-a)^v If[Head[z0]===DirectedInfinity, (-z)^(-exp), (z0-z)^exp],
  192.         a^v If[Head[z0]===DirectedInfinity, z^(-exp), (z-z0)^exp]] *
  193.         SeriesData[z, z0, cla, 0, nmax-nmin, den]^v
  194.     ] /; (nmin =!= 0 && !NumberQ[v])
  195.  
  196.  (* -------------------- Log expansion about Infinity ------------------ *)
  197.  
  198.   (* this should also deal with Log[a_. s_SeriesData], but then it couldn't
  199.      be attached to SeriesData (SeriesData would be inside of a Times) *)
  200.  
  201. SeriesData /:
  202. Log[SeriesData[z_, z0_DirectedInfinity, cl_List, nmin_, nmax_, den_]] :=
  203.     Log[cl[[1]] z^(-nmin/den)] +
  204.     ComposeSeries[Series[Log[z], {z, 1, mterms[1, nmin, nmax, den]}],
  205.         SeriesData[z, z0, cl/cl[[1]], 0, nmax-nmin, den]] /;
  206.     Length[cl] > 0 && simser[z, cl]
  207.  
  208.  (* ------- Various simplification rules and support functions --------- *)
  209.  
  210. logrule = ((Log[n_?((NumberQ[#] && #<0)&)] + Log[-1 + y_]) ->
  211.         (Log[-n] + Log[1-y]));
  212.  
  213. logsplit[Plus[logpart_ + rest___]] :=
  214.     {logpart,Plus[rest]} /; (!FreeQ[logpart,Log] && FreeQ[{rest},Log]);
  215.  
  216. logsplit[rest_] := {0,rest} /; FreeQ[{rest},Log];
  217.  
  218. logsplit[logpart_] := {logpart,0} /; !FreeQ[logpart,Log];
  219.  
  220. intQ[m_] := (NumberQ[m] && (m == Round[m]));
  221.  
  222.     (* faa[] takes a coefficient list and the exponent associated with the
  223.     first term in the coefficient list and returns the value of the series
  224.     evaluated at the center of expansion *)
  225.  
  226. faas[s_SeriesData] := faa[s[[1]], s[[2]], s[[3]], s[[4]]] /; Length[s] == 6
  227.  
  228. faa[x_, a_, cl_, nmin_] :=
  229.     If[nmin == 0,
  230.         If[Length[cl] > 0, Limit[cl[[1]], x -> a], Indeterminate],
  231.         (* else *)
  232.         If[nmin < 0, DirectedInfinity[ ], 0]];
  233.  
  234.     (* mterms[ ] takes a number (either 0 or nonzero: presumably the value of
  235.     the derivative of our function f at the expansion point), nmin, nmax,
  236.     and den and returns a guess as to the number of terms needed in the
  237.     expansion of f w.r.t. a temporary variable in order to get the required
  238.     number of terms in the composed series *)
  239.  
  240. mterms[fp_, s_SeriesData] := mterms[fp, s[[4]], s[[5]], s[[6]]]
  241.  
  242. mterms[fp_,nmin_,nmax_,den_] :=
  243.     If[nmin == 0,
  244.         nmax,
  245.         (* else *)
  246.         If[nmin > 0,
  247.             Floor[nmax/nmin],
  248.             (* else *)
  249.             2 - Ceiling[nmax/nmin]
  250.             ]
  251.         ] + If[TrueQ[fp==0], 1, 0];
  252.  
  253.     (* simser[ ] tests whether the series is a simple series.  Series in x
  254.     containing things like Log[x] in their coefficients cannot be handled
  255.     properly by some of the functions below. *)
  256.  
  257. simser[s_SeriesData] :=
  258.     Length[s] == 6 && (Symbol === Head[s[[1]]]) && FreeQ[s[[3]], s[[1]]]
  259.  
  260. simser[x_, cl_] := (Symbol === Head[x]) && FreeQ[cl, x]
  261.  
  262. Attributes[esssQ] = HoldRest
  263.  
  264. esssQ[f0_, s_] := If[Head[f0] === DirectedInfinity,
  265.             Message[Series::esss, HoldForm[s]]; False,
  266.             True]
  267.  
  268.  (* ----------------------- Airy functions ----------------------------- *)
  269.  
  270.     (* special cases: expansions about +/- Infinity *)
  271.  
  272. airyc[k_] := Module[{j}, Product[j, {j, 2k+1, 6k-1, 2}]/(216^k k!)];
  273.  
  274. airyd[k_] := -(6k+1) airyc[k]/(6k-1);
  275.  
  276. airycdsum[h_, zeta_, n_, i_] :=
  277.     Module[{k}, {Sum[i^(k/2) h[k] zeta^(-k), {k, 0, 2n/3 + 2, 2}],
  278.         Sum[i^((k-1)/2) h[k] zeta^(-k), {k, 1, 2n/3 + 2, 2}]}
  279.     ];
  280.  
  281. AiryAi /: Series[AiryAi[z_], {z_, Infinity, n_Integer}] :=
  282.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
  283.     onSesss = OnQ[Series::esss];
  284.     If[onSesss, Off[Series::esss]];
  285.     {s1, s2} = airycdsum[airyc, zeta, n, 1];
  286.     s1 = 1/2 Pi^(-1/2) z^(-1/4) E^(-zeta) Series[s1-s2, {z, Infinity, n}];
  287.     If[onSesss, On[Series::esss]];
  288.     s1
  289.     ] /; n >= 0
  290.  
  291. AiryAiPrime /: Series[AiryAiPrime[z_], {z_, Infinity, n_Integer}] :=
  292.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
  293.     onSesss = OnQ[Series::esss];
  294.     If[onSesss, Off[Series::esss]];
  295.     {s1, s2} = airycdsum[airyd, zeta, n, 1];
  296.     s1 = -1/2 Pi^(-1/2) z^(1/4) E^(-zeta) Series[s1-s2, {z, Infinity, n}];
  297.     If[onSesss, On[Series::esss]];
  298.     s1
  299.     ] /; n >= 0
  300.  
  301. AiryBi /: Series[AiryBi[z_], {z_, Infinity, n_Integer}] :=
  302.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
  303.     onSesss = OnQ[Series::esss];
  304.     If[onSesss, Off[Series::esss]];
  305.     {s1, s2} = airycdsum[airyc, zeta, n, 1];
  306.     s1 = Pi^(-1/2) z^(-1/4) E^zeta Series[s1+s2, {z, Infinity, n}];
  307.     If[onSesss, On[Series::esss]];
  308.     s1
  309.     ] /; n >= 0
  310.  
  311. AiryBiPrime /: Series[AiryBiPrime[z_], {z_, Infinity, n_Integer}] :=
  312.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onSesss},
  313.     onSesss = OnQ[Series::esss];
  314.     If[onSesss, Off[Series::esss]];
  315.     {s1, s2} = airycdsum[airyd, zeta, n, 1];
  316.     s1 = Pi^(-1/2) z^(1/4) E^zeta Series[s1+s2, {z, Infinity, n}];
  317.     If[onSesss, On[Series::esss]];
  318.     s1
  319.     ] /; n >= 0
  320.  
  321. AiryAi /: Series[AiryAi[z_], {z_, -Infinity, n_Integer}] :=
  322.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
  323.     onSesss = OnQ[Series::esss];
  324.     If[onSesss, Off[Series::esss]];
  325.     onGivar = OnQ[General::ivar];
  326.     If[onGivar, Off[General::ivar]];
  327.     {s1, s2} = airycdsum[airyc, zeta, n, -1];
  328.     s1 = Sin[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
  329.     s2 = Cos[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
  330.     s1 = Pi^(-1/2) (-z)^(-1/4) (s1 - s2);
  331.     If[onGivar, On[General::ivar]];
  332.     If[onSesss, On[Series::esss]];
  333.     s1
  334.     ] /; n >= 0
  335.  
  336. AiryAiPrime /: Series[AiryAiPrime[z_], {z_, -Infinity, n_Integer}] :=
  337.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
  338.     onSesss = OnQ[Series::esss];
  339.     If[onSesss, Off[Series::esss]];
  340.     onGivar = OnQ[General::ivar];
  341.     If[onGivar, Off[General::ivar]];
  342.     {s1, s2} = airycdsum[airyd, zeta, n, -1];
  343.     s1 = Cos[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
  344.     s2 = Sin[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
  345.     s1 = -Pi^(-1/2) (-z)^(1/4) (s1 + s2);
  346.     If[onGivar, On[General::ivar]];
  347.     If[onSesss, On[Series::esss]];
  348.     s1
  349.     ] /; n >= 0
  350.  
  351. AiryBi /: Series[AiryBi[z_], {z_, -Infinity, n_Integer}] :=
  352.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
  353.     onSesss = OnQ[Series::esss];
  354.     If[onSesss, Off[Series::esss]];
  355.     onGivar = OnQ[General::ivar];
  356.     If[onGivar, Off[General::ivar]];
  357.     {s1, s2} = airycdsum[airyc, zeta, n, -1];
  358.     s1 = Cos[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
  359.     s2 = Sin[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
  360.     s1 = Pi^(-1/2) (-z)^(-1/4) (s1 + s2);
  361.     If[onGivar, On[General::ivar]];
  362.     If[onSesss, On[Series::esss]];
  363.     s1
  364.     ] /; n >= 0
  365.  
  366. AiryBiPrime /: Series[AiryBiPrime[z_], {z_, -Infinity, n_Integer}] :=
  367.     Module[{s1, s2, zeta = 2 z^(3/2)/3, onGivar, onSesss},
  368.     onSesss = OnQ[Series::esss];
  369.     If[onSesss, Off[Series::esss]];
  370.     onGivar = OnQ[General::ivar];
  371.     If[onGivar, Off[General::ivar]];
  372.     {s1, s2} = airycdsum[airyd, zeta, n, -1];
  373.     s1 = Sin[zeta + Pi/4] Series[s1, {z, Infinity, n}] /. z -> -z;
  374.     s2 = Cos[zeta + Pi/4] Series[s2, {z, Infinity, n}] /. z -> -z;
  375.     s1 = Pi^(-1/2) (-z)^(1/4) (s1 - s2);
  376.     If[onGivar, On[General::ivar]];
  377.     If[onSesss, On[Series::esss]];
  378.     s1
  379.     ] /; n >= 0
  380.  
  381.     (* generic cases *)
  382.  
  383. AiryAi /: Series[AiryAi[g_],{z_,aa_,nn_Integer}] :=
  384.     AiryAi[Series[g,{z,aa,nn}]] /; nn >= 0
  385.  
  386. AiryAiPrime /: Series[AiryAiPrime[g_],{z_,aa_,nn_Integer}] :=
  387.     AiryAiPrime[Series[g,{z,aa,nn}]] /; nn >= 0
  388.  
  389. AiryBi /: Series[AiryBi[g_],{z_,aa_,nn_Integer}] :=
  390.     AiryBi[Series[g,{z,aa,nn}]] /; nn >= 0
  391.  
  392. AiryBiPrime /: Series[AiryBiPrime[g_],{z_,aa_,nn_Integer}] :=
  393.     AiryBiPrime[Series[g,{z,aa,nn}]] /; nn >= 0
  394.  
  395. AiryAi /: AiryAi[s_SeriesData] :=
  396.     Module[{aa = faas[s]},
  397.     diffeqAiry[AiryAi,False,s,aa] /; esssQ[aa, AiryAi[s]]] /; simser[s]
  398.  
  399. AiryAiPrime /: AiryAiPrime[s_SeriesData] :=
  400.     Module[{aa = faas[s]},
  401.     diffeqAiry[AiryAi,True,s,aa] /; esssQ[aa, AiryAiPrime[s]]] /; simser[s]
  402.  
  403. AiryBi /: AiryBi[s_SeriesData] :=
  404.     Module[{aa = faas[s]},
  405.     diffeqAiry[AiryBi,False,s,aa] /; esssQ[aa, AiryBi[s]]] /; simser[s]
  406.  
  407. AiryBiPrime /: AiryBiPrime[s_SeriesData] :=
  408.     Module[{aa = faas[s]},
  409.     diffeqAiry[AiryBi,True,s,aa] /; esssQ[aa, AiryBiPrime[s]]] /; simser[s]
  410.  
  411. diffeqAiry[head_, der_, s_SeriesData, aa_] :=
  412.     Module[{i, m, f0, f1, f, z = s[[1]], w},
  413.     If[SameQ[aa,Indeterminate], Return[SeriesData[z,aa,{},1,1,1]]];
  414.     m = mterms[1, s] + 1;
  415.     f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i,{i,2,m}],{z,aa,m}];
  416.     f = MapAt[Together,
  417.         (f /. Solve[D[f,{z,2}]==z f,Table[w[i],{i,2,m}]])[[1]],
  418.         {3}];
  419.     f = f /. {f0 -> head[aa], f1 -> head'[aa]};
  420.     If[der, f = D[f, z]];
  421.     ComposeSeries[f, s]
  422.     ];
  423.  
  424.  (* ------------------------ Bessel functions ----------------------------- *)
  425.  
  426.     (* special cases: expansion about Infinity *)
  427.  
  428. bessPQ[v_, z_, n_, i_] :=
  429.     Module[{j, k, mu = 4v^2, p, q},
  430.     (* AS 9.2.9, 9.2.10 *)
  431.     p = Sum[i^(k/2) Product[mu-j^2,{j,1,2k-1,2}]/
  432.         (k! 8^k z^k),{k,0,n+1,2}];
  433.     q = Sum[i^((k-1)/2) Product[mu-j^2,{j,1,2k-1,2}]/
  434.         (k! 8^k z^k),{k,1,n+1,2}];
  435.     {p, q}
  436.     ];
  437.  
  438. BesselJ /: Series[BesselJ[v_, z_], {z_, Infinity, n_Integer}] :=
  439.     Module[{p, q, x = z - (v/2 + 1/4)Pi, onSesss},
  440.     (* AS 9.2.5 *)
  441.     onSesss = OnQ[Series::esss];
  442.     If[onSesss, Off[Series::esss]];
  443.     {p, q} = bessPQ[v, z, n, -1];
  444.     p = Series[p, {z, Infinity, n}];
  445.     q = Series[q, {z, Infinity, n}];
  446.     p = Sqrt[2/(Pi z)] (p Cos[x] - q Sin[x]);
  447.     If[onSesss, On[Series::esss]];
  448.     p
  449.     ] /; n >= 0
  450.  
  451. BesselY /: Series[BesselY[v_, z_], {z_, Infinity, n_Integer}] :=
  452.     Module[{p, q, x = z - (v/2 + 1/4)Pi, onSesss},
  453.     (* AS 9.2.6 *)
  454.     onSesss = OnQ[Series::esss];
  455.     If[onSesss, Off[Series::esss]];
  456.     {p, q} = bessPQ[v, z, n, -1];
  457.     p = Series[p, {z, Infinity, n}];
  458.     q = Series[q, {z, Infinity, n}];
  459.     p = Sqrt[2/(Pi z)] (p Sin[x] + q Cos[x]);
  460.     If[onSesss, On[Series::esss]];
  461.     p
  462.     ] /; n >= 0
  463.  
  464. BesselI /: Series[BesselI[v_, z_], {z_, Infinity, n_Integer}] :=
  465.     Module[{p, q, onSesss},
  466.     (* AS 9.7.1, corrected *)
  467.     onSesss = OnQ[Series::esss];
  468.     If[onSesss, Off[Series::esss]];
  469.     {p, q} = bessPQ[v, z, n, 1];
  470.     p = Series[p, {z, Infinity, n}];
  471.     q = Series[q, {z, Infinity, n}];
  472.     p = Sqrt[1/(2 Pi z)] (E^z (p-q) + E^(I (v + 1/2) Pi - z) (p+q));
  473.     If[onSesss, On[Series::esss]];
  474.     p
  475.     ] /; n >= 0
  476.  
  477. BesselK /: Series[BesselK[v_, z_], {z_, Infinity, n_Integer}] :=
  478.     Module[{p, q, onSesss},
  479.     (* AS 9.7.2 *)
  480.     onSesss = OnQ[Series::esss];
  481.     If[onSesss, Off[Series::esss]];
  482.     {p, q} = bessPQ[v, z, n, 1];
  483.     p = Series[p+q, {z, Infinity, n}];
  484.     p = Sqrt[Pi/(2 z)] E^(-z) p;
  485.     If[onSesss, On[Series::esss]];
  486.     p
  487.     ] /; n >= 0
  488.  
  489.     (* generic case *)
  490.  
  491. besser[v_,z_,a_,n_,fa_,fpa_,m_] := 
  492.     Module[{i, f0, f1, df, f},
  493.     f = Series[f0 + f1(z-a) + Sum[w[i] (z-a)^i, {i,2,n}],{z,a,n}];
  494.     df = D[f,z];
  495.     f = MapAt[Expand,
  496.         (f /. Solve[z^2 D[df,z] + z df + (m z^2 - v^2) f == 0,
  497.             Table[w[i],{i,2,n}]])[[1]],
  498.         {3}];
  499.     f /. {f0 -> fa, f1 -> fpa}
  500.     ];
  501.  
  502. BesselI /: Series[BesselI[v_,g_],{z_,aa_,nn_Integer}] :=
  503.     BesselI[v,Series[g,{z,aa,nn}]] /; nn >= 0
  504.  
  505. BesselI /: Literal[BesselI[v_, zser_SeriesData]] :=
  506.     Module[{m, fp, f, aa = faas[zser], z},
  507.     z = zser[[1]];
  508.     If[TrueQ[aa==0],
  509.         m = mterms[1, zser];
  510.         If[NumberQ[v], m -= Floor[Re[v]]];
  511.         f = Series[Hypergeometric0F1Regularized[v+1,z^2/4], {z,aa,m}];
  512.         (zser/2)^v ComposeSeries[f, zser],
  513.         (* else *)
  514.         fp = D[BesselI[v,z],z] /. z->aa;
  515.         m = mterms[fp, zser];
  516.         f = besser[v,z,aa,m,BesselI[v,aa],fp,-1];
  517.         ComposeSeries[f, zser]
  518.         ] /;
  519.     esssQ[aa, BesselI[v, zser]]
  520.     ] /; simser[zser]
  521.  
  522. BesselJ /: Series[BesselJ[v_,g_],{z_,aa_,nn_Integer}] :=
  523.     BesselJ[v,Series[g,{z,aa,nn}]] /; nn >= 0
  524.  
  525. BesselJ /: Literal[BesselJ[v_, zser_SeriesData]] :=
  526.     Module[{m, fp, f, aa = faas[zser], z},
  527.     z = zser[[1]];
  528.     If[TrueQ[aa==0],
  529.         m = mterms[1, zser];
  530.         If[NumberQ[v], m -= Floor[Re[v]]];
  531.         f = Series[Hypergeometric0F1Regularized[v+1,-z^2/4], {z,aa,m}];
  532.         (zser/2)^v ComposeSeries[f, zser],
  533.         (* else *)
  534.         fp = D[BesselJ[v,z],z] /. z->aa;
  535.         m = mterms[fp, zser];
  536.         f = besser[v,z,aa,m,BesselJ[v,aa],fp,1];
  537.         ComposeSeries[f,zser]
  538.         ] /;
  539.     esssQ[aa, BesselJ[v, zser]]
  540.     ] /; simser[zser]
  541.  
  542. BesselK /: Series[BesselK[v_,g_],{z_,aa_,nn_Integer}] :=
  543.     BesselK[v,Series[g,{z,aa,nn}]] /; nn >= 0
  544.  
  545. BesselK /: Literal[BesselK[v_, zser_SeriesData]] :=
  546.     Module[{m, k, nu = v, fp, f, aa = faas[zser], z},
  547.     z = zser[[1]];
  548.     (If[Re[nu] < 0, nu = -nu];
  549.     If[TrueQ[aa==0],
  550.         If[intQ[nu],
  551.         m = mterms[1, zser];
  552.         f = Series[1/2 (2/z)^nu *
  553.             Sum[(nu-k-1)! (-z^2/4)^k/k!, {k,0,nu-1}] -
  554.             (-1)^nu (Log[z/2] BesselI[nu,z] - 1/2(z/2)^nu *
  555.             Sum[(PolyGamma[k+1]+PolyGamma[k+1+nu])(z^2/4)^k/
  556.             (k!(nu+k)!), {k,0,m-Abs[Round[nu]]}]),{z,0,m}];
  557.         ComposeSeries[f,zser],
  558.           (* else (!intQ[nu]) *)
  559.         Pi(BesselI[-nu,zser] - BesselI[nu,zser])/(2 Sin[nu Pi])
  560.         ],
  561.       (* else (aa != 0)*)
  562.         fp = D[BesselK[nu,z],z] /. z->aa;
  563.         m = mterms[fp, zser];
  564.         f = besser[nu,z,aa,m,BesselK[nu,aa],fp,-1];
  565.         ComposeSeries[f,zser]
  566.         ]) /;
  567.     esssQ[aa, BesselK[v, zser]]
  568.     ] /; simser[zser]
  569.  
  570. BesselY /: Series[BesselY[v_,g_],{z_,aa_,nn_Integer}] :=
  571.     BesselY[v,Series[g,{z,aa,nn}]] /; nn >= 0
  572.  
  573. BesselY /: Literal[BesselY[v_, zser_SeriesData]] :=
  574.     Module[{m, k, nu = v, fp, f, aa = faas[zser], sign = 1, z},
  575.     z = zser[[1]];
  576.     (If[intQ[nu/2] && nu < 0, nu = -nu; sign = (-1)^nu];
  577.     If[TrueQ[aa==0],
  578.         If[intQ[nu],
  579.         m = mterms[1, zser];
  580.         f = Series[-1/Pi (2/z)^nu *
  581.             Sum[(nu-k-1)! (z^2/4)^k/k!, {k,0,nu-1}] +
  582.             2/Pi Log[z/2] BesselJ[nu,z] - 1/Pi (z/2)^nu *
  583.             Sum[(PolyGamma[k+1]+PolyGamma[k+1+nu]) *
  584.             (-z^2/4)^k/(k!(nu+k)!), {k,0,m-Abs[Round[nu]]}],
  585.             {z,0,m}];
  586.         ComposeSeries[f,zser],
  587.           (* else (!intQ[nu]) *)
  588.         (BesselJ[v,zser] Cos[v Pi] - BesselJ[-v,zser])/Sin[v Pi]
  589.         ],
  590.       (* else (aa != 0)*)
  591.         fp = D[BesselY[nu,z],z] /. z->aa;
  592.         m = mterms[fp, zser];
  593.         f = besser[nu,z,aa,m,BesselY[nu,aa],fp,1];
  594.         ComposeSeries[f,zser]
  595.         ]) /;
  596.     esssQ[aa, BesselY[v, zser]]
  597.     ] /; simser[zser]
  598.  
  599.  (* ------------------------- Beta functions ----------------------------- *)
  600.  
  601. betaseries[sf1_,sf2_,a_,b_,regflag_,z_,z0_] :=
  602.     Module[{aa=a,bb=b,m1,m2,f1=0,f2=0,f01,f02,t,anint,bnint,logs1,logs2},
  603.     If[intQ[a],aa = Round[aa]];
  604.     anint = IntegerQ[aa] && aa<=0;
  605.     If[intQ[b],bb = Round[bb]];
  606.     bnint = IntegerQ[bb] && bb<=0;
  607.     If[sflag1 = (Head[sf1] === SeriesData),
  608.         f01 = faas[sf1];
  609.         m1 = mterms[1, sf1],
  610.         (* else *)
  611.         f01 = sf1;
  612.         m1 = Infinity;
  613.         ];
  614.     If[sflag2 = (Head[sf2] === SeriesData),
  615.         f02 = faas[sf2];
  616.         m2 = mterms[1,sf2],
  617.         (* else *)
  618.         f02 = sf2;
  619.         m2 = Infinity;
  620.         ];
  621.     logs1 = sflag1 && TrueQ[f01==0];
  622.     logs2 = sflag2 && TrueQ[f02==1];
  623.     If[sflag1, m1 += Max[If[TrueQ[f01==0], -Round[aa],0],
  624.         If[TrueQ[f01==1], -Round[bb],0]]];
  625.     If[sflag2, m2 += Max[If[TrueQ[f02==0], -Round[aa],0],
  626.         If[TrueQ[f02==1], -Round[bb],0]]];
  627.     If[anint && bnint && (logs1 || logs2),
  628.         If[regflag, Return[0]];
  629.         (* this is slow, but I don't know how else to do it *)
  630.         f1 = Integrate[t^(aa-1)(1-t)^(bb-1),t] /. Log[-1+t]->Log[1-t];
  631.         f2 = If[sflag2, ComposeSeries[Series[f1,{t,f02,m2}],sf2],
  632.             (* else *) f1 /. t->f02];
  633.         f1 = If[sflag1, ComposeSeries[Series[f1,{t,f01,m1}],sf1],
  634.             (* else *) f1 /. t->f01];
  635.         Return[f2-f1]
  636.         ];
  637.     If[anint && logs1,
  638.         Return[betahyper[1-f02,1-sf2,m2,1-f01,1-sf1,m1,bb,aa,regflag]]];
  639.     If[bnint && logs2,
  640.         Return[betahyper[f01,sf1,m1,f02,sf2,m2,aa,bb,regflag]]];
  641.     If[sflag1, {f01,f1} = betser0[aa,bb,sf1,f01,m1]];
  642.     If[sflag2, {f02,f2} = betser0[aa,bb,sf2,f02,m2]];
  643.     If[regflag,
  644.         BetaRegularized[f01,f02,a,b] + (f2 - f1)/Beta[a,b],
  645.         Beta[f01,f02,a,b] + f2 - f1
  646.         ]
  647.     ]
  648.  
  649. betahyper[f01_,sf1_,m1_,f02_,sf2_,m2_,a_,b_,regflag_] :=
  650.     Module[{f1,f2,t,beta = Beta[a,b]},
  651.     If[regflag && (Head[beta] === DirectedInfinity), Return[0]];
  652.     (* AS 6.6.8; a is not a non-positive integer! *)
  653.     f2 = ComposeSeries[t^a/a Series[Hypergeometric2F1[
  654.         a,1-b,a+1,t],{t,f02,m2}],sf2];
  655.     f1 = If[!SameQ[Head[sf1],SeriesData],
  656.         f01^a/a Hypergeometric2F1[a,1-b,a+1,f01],
  657.         (* else *)
  658.         ComposeSeries[t^a/a Series[Hypergeometric2F1[
  659.             a,1-b,a+1,t],{t,f01,m1}],sf1]];
  660.     If[regflag, (f2 - f1)/beta, f2 - f1] /. logrule
  661.     ];
  662.  
  663. betser0[a_,b_,sf_,f0_,m_] :=
  664.     Module[{t, f},
  665.     f = Integrate[Series[(t+f0)^(a-1)(1-f0-t)^(b-1),{t,0,m}],t];
  666.     {f0,ComposeSeries[f, sf-f0]}
  667.     ];
  668.  
  669. Beta /: (* complete beta function *)
  670. Series[Beta[f1z_,f2z_],{z_,aa_,nn_Integer}] :=
  671.     Module[{sf1 = Series[f1z,{z,aa,nn}],sf2 = Series[f2z,{z,aa,nn}]},
  672.     Gamma[sf1] Gamma[sf2] / Gamma[sf1+sf2]] /; nn >= 0
  673.  
  674. Beta /: (* complete beta function *)
  675. Literal[Beta[sf1_SeriesData, sf2_SeriesData]]:=
  676.     Gamma[sf1] Gamma[sf2] / Gamma[sf1+sf2] /;
  677.     (simser[sf1] && simser[sf2] &&
  678.     sf1[[1]] === sf2[[1]] && sf1[[2]] == sf2[[2]])
  679.  
  680. Beta /: (* complete beta function *)
  681. Literal[Beta[sf_SeriesData, a_]]:=
  682.     Module[{sa},
  683.     sa = Series[a, {sf[[1]], sf[[2]],
  684.             Round[(sf[[5]] - Min[0, sf[[4]]])/sf[[6]]]}];
  685.     Gamma[sf] Gamma[sa] / Gamma[sf+sa]
  686.     ] /; simser[sf]
  687.  
  688. Beta /: (* complete beta function *)
  689. Literal[Beta[a_, sf_SeriesData]] :=
  690.     Module[{sa},
  691.     sa = Series[a, {sf[[1]], sf[[2]],
  692.             Round[(sf[[5]] - Min[0, sf[[4]]])/sf[[6]]]}];
  693.     Gamma[sf] Gamma[sa] / Gamma[sf+sa]
  694.     ] /; simser[sf]
  695.  
  696. Beta /: (* incomplete beta function *)
  697. Series[Beta[f1z_:0,f2z_,a_,b_],{z_,aa_,nn_Integer}] :=
  698.     Beta[Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}],a,b] /;
  699.         (nn >= 0 && FreeQ[a,z] && FreeQ[b,z])
  700.  
  701. Beta /: (* incomplete beta function *)
  702. Literal[Beta[s_SeriesData, a_, b_]] :=
  703.     betaseries[0, s, a, b, False, s[[1]], s[[2]]] /;
  704.     (simser[s] && FreeQ[a,s[[1]]] && FreeQ[b,s[[1]]])
  705.  
  706. Beta /: (* incomplete beta function *)
  707. Literal[Beta[s1_SeriesData, s2_SeriesData, a_, b_]] :=
  708.     betaseries[s1, s2, a, b, False, s1[[1]], s1[[2]]] /;
  709.     (simser[s1] && simser[s2] && s1[[1]] === s2[[1]] &&
  710.     s1[[2]] == s2[[2]] && FreeQ[a,s1[[1]]] && FreeQ[b,s1[[1]]])
  711.  
  712. Beta /: (* incomplete beta function *)
  713. Literal[Beta[f1_SeriesData,f2_,a_,b_]] :=
  714.     betaseries[f1,
  715.     Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
  716.     a,b,False,f1[[1]],f1[[2]]] /;
  717.     (simser[f1] && FreeQ[a,f1[[1]]] && FreeQ[b,f1[[1]]])
  718.  
  719. Beta /: (* incomplete beta function *)
  720. Literal[Beta[f1_,f2_SeriesData,a_,b_]] :=
  721.     betaseries[
  722.     Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
  723.     f2,a,b,False,f2[[1]],f2[[2]]] /;
  724.     (simser[f2] && FreeQ[a,f2[[1]]] && FreeQ[b,f2[[1]]])
  725.  
  726. BetaRegularized /: (* incomplete beta function *)
  727. Series[BetaRegularized[f1z_:0,f2z_,a_,b_],{z_,aa_,nn_Integer}] :=
  728.     BetaRegularized[Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}],a,b] /;
  729.     (nn >= 0 && FreeQ[a,z] && FreeQ[b,z])
  730.  
  731. BetaRegularized /: (* incomplete beta function *)
  732. Literal[BetaRegularized[s2_SeriesData, a_, b_]] :=
  733.     betaseries[0, s2, a, b, True, s2[[1]], s2[[2]]] /;
  734.     (simser[s2] && FreeQ[a,s2[[1]]] && FreeQ[b,s2[[1]]])
  735.  
  736. BetaRegularized /: (* incomplete beta function *)
  737. Literal[BetaRegularized[s1_SeriesData, s2_SeriesData, a_, b_]] :=
  738.     betaseries[s1, s2, a, b, True, s1[[1]], s1[[2]]] /;
  739.     (simser[s1] && simser[s2] && s1[[1]] === s2[[1]] &&
  740.     s1[[2]] == s2[[2]] && FreeQ[a,s1[[1]]] && FreeQ[b,s1[[1]]])
  741.  
  742. BetaRegularized /: (* incomplete beta function *)
  743. Literal[BetaRegularized[f1_SeriesData,f2_,a_,b_]] :=
  744.     betaseries[f1,
  745.     Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
  746.     a,b,True,f1[[1]],f1[[2]]] /;
  747.     (simser[f1] && FreeQ[a,f1[[1]]] && FreeQ[b,f1[[1]]])
  748.  
  749. BetaRegularized /: (* incomplete beta function *)
  750. Literal[BetaRegularized[f1_,f2_SeriesData,a_,b_]] :=
  751.     betaseries[
  752.     Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
  753.     f2,a,b,True,f2[[1]],f2[[2]]] /;
  754.     (simser[f2] && FreeQ[a,f2[[1]]] && FreeQ[b,f2[[1]]])
  755.  
  756.  (* --------------------- Binomial function ---------------------------- *)
  757.  
  758. Series[Binomial[x_,y_],{z_,aa_,nn_Integer}] :=
  759.     Module[{ys = Series[y, {z, aa, nn}]},
  760.     1/(ys Beta[ys, 1+Series[x, {z, aa, nn}] - ys])] /; (nn >= 0)
  761.  
  762. SeriesData /:
  763. Literal[Binomial[xs_SeriesData, y_]] :=
  764.     Module[{ys = Series[y, {xs[[1]], xs[[2]],
  765.         Round[(xs[[5]] - Min[0, xs[[4]]])/xs[[6]]]}]},
  766.     1/(ys Beta[ys, xs + (1-ys)])
  767.     ] /; simser[xs]
  768.  
  769. SeriesData /:
  770. Literal[Binomial[x_, ys_SeriesData]] :=
  771.     Module[{xs = Series[x, {ys[[1]], ys[[2]],
  772.         Round[(ys[[5]] - Min[0, ys[[4]]])/ys[[6]]]}]},
  773.     1/(ys Beta[ys, 1+xs-ys])
  774.     ] /; simser[ys]
  775.  
  776. SeriesData /:
  777. Literal[Binomial[xs_SeriesData, ys_SeriesData]] :=
  778.     1/(ys Beta[ys, 1+xs-ys]) /; simser[xs] && simser[ys]
  779.  
  780.  (* -------------------- Chebyshev functions --------------------------- *)
  781.  
  782. ChebyshevT /: Series[ChebyshevT[n_,fz_],{z_,aa_,nn_Integer}] :=
  783.     ChebyshevT[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])
  784.  
  785. ChebyshevT /: ChebyshevT[n_,s_SeriesData] := Cos[n ArcCos[s]] /; simser[s]
  786.  
  787. ChebyshevU /: Series[ChebyshevU[n_,fz_],{z_,aa_,nn_Integer}] :=
  788.     ChebyshevU[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])
  789.  
  790. ChebyshevU /: ChebyshevU[n_,s_SeriesData] :=
  791.     Module[{theta = ArcCos[s]}, Sin[(n+1) theta]/Sin[theta]] /; simser[s]
  792.  
  793.  (* ------------------- Exponential integrals -------------------------- *)
  794.  
  795.     (* special cases: expansions about Infinity *)
  796.  
  797. ExpIntegralE /: Series[ExpIntegralE[n_, z_], {z_, Infinity, nn_Integer}] :=
  798.     Module[{j, onSesss},
  799.     onSesss = OnQ[Series::esss];
  800.     If[onSesss, Off[Series::esss]];
  801.     j = E^(-z) Series[Sum[(-1)^j Pochhammer[n, j] z^(-j-1), {j, 0, nn}],
  802.             {z, Infinity, nn}];
  803.     If[onSesss, On[Series::esss]];
  804.     j
  805.     ] /; (nn >= 0 && FreeQ[n, z])
  806.  
  807. ExpIntegralEi /: Series[ExpIntegralEi[z_], {z_, Infinity, nn_Integer}] :=
  808.     Module[{j, onSesss},
  809.     onSesss = OnQ[Series::esss];
  810.     If[onSesss, Off[Series::esss]];
  811.     j = E^z Series[Sum[j! z^(-j-1), {j, 0, nn}], {z, Infinity, nn}];
  812.     If[onSesss, On[Series::esss]];
  813.     j
  814.     ] /; nn >= 0
  815.  
  816.     (* generic cases *)
  817.  
  818. ExpIntegralE /: Series[ExpIntegralE[n_,fz_],{z_,aa_,nn_Integer}] :=
  819.     ExpIntegralE[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])
  820.  
  821. ExpIntegralE /:
  822. Literal[ExpIntegralE[n_, fz_SeriesData]] :=
  823.     fz^(n-1) Gamma[1-n,fz] /; (simser[fz] && FreeQ[n,fz[[1]]])
  824.  
  825. ExpIntegralEi /: Series[ExpIntegralEi[fz_],{z_,aa_,nn_Integer}] :=
  826.     ExpIntegralEi[Series[fz,{z,aa,nn}]] /; nn >= 0
  827.  
  828. ExpIntegralEi /:
  829. Literal[ExpIntegralEi[s_SeriesData]] :=
  830.     Module[{f0 = faas[s], m = mterms[1, s], temp, fs},
  831.     (fs = If[TrueQ[f0==0],
  832.         EulerGamma+Log[temp]+Sum[temp^k/(k k!),{k,m}] + O[temp]^(m+1),
  833.         (* else *)
  834.         Integrate[Series[Exp[temp+f0]/(temp+f0),{temp,0,m}],temp] +
  835.             + ExpIntegralEi[f0]
  836.         ];
  837.     ComposeSeries[fs, s-f0]) /; esssQ[f0, ExpIntegralEi[s]]
  838.     ] /; simser[s]
  839.  
  840.  (* ----------------- Factorial and Gamma  functions ---------------------- *)
  841.  
  842.     (* special cases: expansions about Infinity *)
  843.  
  844. Gamma /: Series[Gamma[z_],{z_, Infinity, nn_Integer}] :=
  845.     Module[{s, i, onSesss},
  846.     onSesss = OnQ[Series::esss];
  847.     If[onSesss, Off[Series::esss]];
  848.     s = Sum[BernoulliB[i]/(i(i-1) z^(i-1)), {i, 2, nn+2, 2}];
  849.     s = z^z E^(-z) Sqrt[2 Pi/z] E^Series[s, {z, Infinity, nn}];
  850.     If[onSesss, On[Series::esss]];
  851.     s
  852.     ] /; nn >= 0
  853.  
  854. Factorial/: Series[Factorial[z_],{z_, Infinity, nn_Integer}] :=
  855.     Module[{s, i, onSesss},
  856.     onSesss = OnQ[Series::esss];
  857.     If[onSesss, Off[Series::esss]];
  858.     s = Sum[BernoulliB[i]/(i(i-1) (z+1)^(i-1)), {i, 2, nn+2, 2}];
  859.     s = (z+1)^z E^(-z-1) Sqrt[2 Pi (z+1)] E^Series[s, {z, Infinity, nn}];
  860.     If[onSesss, On[Series::esss]];
  861.     s
  862.     ] /; nn >= 0
  863.  
  864.     (* generic cases *)
  865.  
  866. gammaseries[a_,sf1_,sf2_,regflag_,z_,z0_] :=
  867.     Module[{f1,f2,f01,f02},
  868.     If[Head[sf1] === SeriesData,
  869.         {f01,f1} = gamser0[a,sf1], f01 = sf1; f1 = 0];
  870.     If[Head[sf2] === SeriesData,
  871.         {f02,f2} = gamser0[a,sf2], f02 = sf2; f2 = 0];
  872.     If[f01 === $Failed || f02 === $Failed,
  873.         HoldForm[Head[a, sf1, sf2]] /.
  874.             Head -> If[regflag, GammaRegularized, Gamma],
  875.         (* else *)
  876.         If[regflag, GammaRegularized[a,f01,f02] + (f2 - f1)/Gamma[a],
  877.             Gamma[a,f01,f02] + f2 - f1]
  878.         ]
  879.     ];
  880.  
  881. gamser0[a_, sf_] :=
  882.     Module[{t, f, ma, i, m = mterms[1,sf], f0 = faas[sf]},
  883.     If[Head[f0] === DirectedInfinity,
  884.         If[TrueQ[f0 == Infinity],
  885.             Return[{f0, 0}],
  886.             (* else *)
  887.             Message[Series::esss, HoldForm[Gamma[a, sf]]];
  888.             Return[{$Failed, $Failed}]
  889.             ]
  890.         ];
  891.     If[TrueQ[f0==0] && intQ[a] && a <= 0,
  892.         ma = -Round[a];
  893.         m += ma;
  894.         f = Sum[If[i==ma,0,(-1)^i t^(i-ma)/((i-ma)i!)],{i,0,m}] +
  895.                 O[t]^(m+1);
  896.         f = ComposeSeries[f,sf] + (-1)^ma((Log[sf] /. logrule) -
  897.             PolyGamma[ma+1])/ma!;
  898.         {Infinity,f},
  899.         (* else *)
  900.         f = Integrate[Series[Exp[-(t+f0)] (t+f0)^(a-1), {t,0,m}],t];
  901.         {f0,ComposeSeries[f, sf-f0]}
  902.         ]
  903.     ];
  904.  
  905. Series[Factorial[fz_],{z_,aa_,nn_Integer}] :=
  906.     Gamma[Series[fz+1,{z,aa,nn}]] /; nn >= 0
  907.  
  908. SeriesData /: Literal[Factorial[s_SeriesData]] := Gamma[1+s] /; simser[s]
  909.  
  910. Series[Gamma[fz_],{z_,aa_,nn_Integer}] :=
  911.     Gamma[Series[fz,{z,aa,nn}]] /; nn >= 0
  912.  
  913. SeriesData /:  (* complete gamma function *)
  914. Literal[Gamma[fz_SeriesData]] :=
  915.     Module[{f0 = faas[fz], i},
  916.     If[NumberQ[f0] && f0 < 1/2,
  917.         Pi/(Sin[Pi fz] Gamma[1-fz]),
  918.         (* else *)
  919.         Sum[Derivative[i][Gamma][f0](fz-f0)^i/(i!),{i,0,mterms[1,fz]}]
  920.         ] /;
  921.     esssQ[f0, Gamma[fz]]
  922.     ] /; simser[fz]
  923.  
  924. Series[Gamma[a_,f1z_,f2z_:Infinity],{z_,aa_,nn_Integer}] :=
  925.     Gamma[a,Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}]] /;
  926.         (nn >= 0 && FreeQ[a,z])
  927.  
  928. SeriesData /:  (* incomplete gamma function *)
  929. Literal[Gamma[a_, s_SeriesData]]:=
  930.     gammaseries[a,s,Infinity,False,s[[1]],s[[2]]] /;
  931.     (simser[s] && FreeQ[a,s[[1]]])
  932.  
  933. SeriesData /:  (* incomplete gamma function *)
  934. Literal[Gamma[a_, s1_SeriesData, s2_SeriesData]]:=
  935.     gammaseries[a, s1, s2, False, s1[[1]], s1[[2]]] /;
  936.     (FreeQ[a,z] && s1[[1]] === s1[[1]] && s1[[2]] == s2[[2]] &&
  937.     simser[z, cl1] && simser[z, cl2])
  938.  
  939. SeriesData /:  (* incomplete gamma function *)
  940. Literal[Gamma[a_,f1_,f2_SeriesData]] :=
  941.     gammaseries[a,
  942.     Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
  943.     f2, False, f2[[1]], f2[[2]]] /;
  944.     (simser[f2] && FreeQ[a,f2[[1]]])
  945.  
  946. SeriesData /:  (* incomplete gamma function *)
  947. Literal[Gamma[a_, f1_SeriesData, f2_:Infinity]] :=
  948.     gammaseries[a, f1,
  949.     Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
  950.     False, f1[[1]], f1[[2]]] /;
  951.     (simser[f1] && FreeQ[a,f1[[1]]])
  952.  
  953. GammaRegularized /: (* incomplete gamma function *)
  954. Series[GammaRegularized[a_,f1z_,f2z_:Infinity],{z_,aa_,nn_Integer}] :=
  955.     GammaRegularized[a,Series[f1z,{z,aa,nn}],Series[f2z,{z,aa,nn}]] /;
  956.         (nn >= 0 && FreeQ[a,z])
  957.  
  958. GammaRegularized /: (* incomplete gamma function *)
  959. Literal[GammaRegularized[a_, s_SeriesData]]:=
  960.     gammaseries[a,s,Infinity,True,s[[1]],s[[2]]] /;
  961.     (simser[s] && FreeQ[a,s[[1]]])
  962.  
  963. GammaRegularized /: (* incomplete gamma function *)
  964. Literal[GammaRegularized[a_, s1_SeriesData, s2_SeriesData]]:=
  965.     gammaseries[a,s1,s2,True,s1[[1]],s1[[2]]] /;
  966.     (FreeQ[a,z] && s1[[1]] === s2[[1]] && s1[[2]] == s2[[2]] &&
  967.     simser[z,cl1] && simser[z,cl2])
  968.  
  969. GammaRegularized /: (* incomplete gamma function *)
  970. Literal[GammaRegularized[a_, f1_, f2_SeriesData]]:=
  971.     gammaseries[a,
  972.     Series[f1,{f2[[1]],f2[[2]],Round[(f2[[5]]-Min[0,f2[[4]]])/f2[[6]]]}],
  973.     f2, True, f2[[1]], f2[[2]]] /;
  974.     (simser[f2] && FreeQ[a,f2[[1]]])
  975.  
  976. GammaRegularized /: (* incomplete gamma function *)
  977. Literal[GammaRegularized[a_, f1_SeriesData, f2_:Infinity]] :=
  978.     gammaseries[a, f1,
  979.     Series[f2,{f1[[1]],f1[[2]],Round[(f1[[5]]-Min[0,f1[[4]]])/f1[[6]]]}],
  980.     True, f1[[1]], f1[[2]]] /;
  981.     (simser[f1] && FreeQ[a,f1[[1]]])
  982.  
  983.  (* ----------------- Gegenbauer polynomials --------------------------- *)
  984.  
  985. GegenbauerC /: Series[GegenbauerC[n_,m_,fz_],{z_,aa_,nn_Integer}] :=
  986.     GegenbauerC[n,m,Series[fz,{z,aa,nn}]] /;
  987.         (nn >= 0 && FreeQ[n,z] && FreeQ[m,z])
  988.  
  989. GegenbauerC /:
  990. Literal[GegenbauerC[n_, 0, s_SeriesData]] :=
  991.     2/n ChebyshevT[n,s] /; (simser[s] && FreeQ[n,s[[1]]])
  992.  
  993. GegenbauerC /:
  994. Literal[GegenbauerC[n_, m_, s_SeriesData]] :=
  995.     Gamma[n+2m] Sqrt[Pi]/(n! 2^(2m-1) Gamma[m]) *
  996.         Hypergeometric2F1Regularized[-n, n+2m, m+1/2, (1-s)/2] /;
  997.     (simser[s] && FreeQ[n,s[[1]]] && FreeQ[m,s[[1]]])
  998.  
  999.  (* --------------------- Hermite polynomials --------------------------- *)
  1000.  
  1001. HermiteH /: Series[HermiteH[n_,fz_],{z_,aa_,nn_Integer}] :=
  1002.     HermiteH[n,Series[fz,{z,aa,nn}]] /; (nn >= 0 && FreeQ[n,z])
  1003.  
  1004. HermiteH /: Literal[HermiteH[n_, fz_SeriesData]] :=
  1005.     2^n fz HypergeometricU[(1-n)/2,3/2,fz^2] /;
  1006.     (simser[fz] && FreeQ[n,fz[[1]]])
  1007.  
  1008.  (* ------------------- Hypergeometric functions ------------------------- *)
  1009.  
  1010.     (* utility for expansions about Infinity *)
  1011.  
  1012. hyper2F0r[a_, b_, z_, nn_] :=
  1013.     Module[{n}, Sum[Pochhammer[a,n] Pochhammer[b,n] z^(-n)/n!, {n, 0, nn}]];
  1014.  
  1015.     (* utility for distinction between regularized and unregularized *)
  1016.  
  1017. regden[a_, i_, regflag_] := If[regflag, Gamma[a+i], Pochhammer[a, i]]
  1018.  
  1019.  (* --- Hypergeometric 0F1 --- *)
  1020.  
  1021.     (* special case: expansions about Infinity *)
  1022.  
  1023. hyper0F1inf[a_, z_, nn_, regflag_] :=
  1024.     Module[{s1, s2, zz = 4 Sqrt[z], aa = a - 1/2, onSesss},
  1025.     (* based on AS 13.5.1, with 0F1[a,z] == E^(-zz/2) 1F1[aa, 2aa, zz] *)
  1026.     onSesss = OnQ[Series::esss];
  1027.     If[onSesss, Off[Series::esss]];
  1028.     s1 = Series[hyper2F0r[aa, 1-aa, -zz, nn], {z, Infinity, nn}];
  1029.     s2 = Series[hyper2F0r[aa, 1-aa, zz, nn], {z, Infinity, nn}];
  1030.     s1 = (E^(I Pi aa-zz/2)s1+E^(zz/2)s2) zz^(-aa) 4^(a-1)
  1031.                 If[regflag, 1, Gamma[a]]/Sqrt[Pi];
  1032.     If[onSesss, On[Series::esss]];
  1033.     s1
  1034.     ] /; (nn >= 0 && FreeQ[a, z])
  1035.  
  1036. Hypergeometric0F1 /:
  1037. Series[Hypergeometric0F1[a_, z_], {z_, Infinity, nn_Integer}] :=
  1038.     hyper0F1inf[a, z, nn, False (* not regularized *)] /;
  1039.     (nn >= 0 && FreeQ[a, z])
  1040.  
  1041. Hypergeometric0F1Regularized /:
  1042. Series[Hypergeometric0F1Regularized[a_, z_], {z_, Infinity, nn_Integer}] :=
  1043.     hyper0F1inf[a, z, nn, True (* regularized *)] /; (nn >= 0 && FreeQ[a, z])
  1044.  
  1045.     (* generic case *)
  1046.  
  1047. diffeq0F1[a_,fz_,aa_,z_,m_,head_] :=
  1048.     Module[{i, f0, f1, df, f},
  1049.     f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i, {i,2,m}],{z,aa,m}];
  1050.     df = D[f,z];
  1051.     f = MapAt[Expand,
  1052.         (f /. Solve[z D[df,z]+a df-f == 0,Table[w[i],{i,2,m}]])[[1]],
  1053.         {3}];
  1054.     ComposeSeries[f /. {f0 -> head[a,aa],
  1055.                 f1 -> Derivative[0,1][head][a,aa]},
  1056.             fz]
  1057.     ];
  1058.  
  1059. hyper0F1[a_, fz_SeriesData, regflag_] :=
  1060.     Module[{i, m, f0 = faas[fz], tmp, z=fz[[1]]},
  1061.     If[intQ[a] && a <= 0,
  1062.         If[regflag,
  1063.         m = 1-Round[a];
  1064.         fz^m hyper0F1[m+1,fz,True],
  1065.           (* else *)
  1066.         ComplexInfinity
  1067.         ],
  1068.       (* else *)
  1069.         If[Head[f0] === DirectedInfinity,
  1070.         tmp = 2 Sqrt[fz];
  1071.         If[regflag,
  1072.             Hypergeometric1F1Regularized[a-1/2,2a-1,2tmp]*
  1073.                 4^(a-1) Gamma[a-1/2]/Sqrt[Pi],
  1074.           (* else *)
  1075.             Hypergeometric1F1[a-1/2,2a-1,2tmp]
  1076.             ] Exp[-tmp],
  1077.           (* else *)
  1078.         m = mterms[1, fz];
  1079.         If[TrueQ[f0==0],
  1080.             ComposeSeries[
  1081.             Sum[z^i/(i! regden[a,i,regflag]), {i,0,m}] + O[z]^(m+1),
  1082.             fz],
  1083.           (* else *)
  1084.             diffeq0F1[a,fz,f0,z,m, If[regflag,
  1085.             Hypergeometric0F1Regularized, Hypergeometric0F1]]
  1086.             ]
  1087.         ]
  1088.         ]
  1089.     ];
  1090.  
  1091. Hypergeometric0F1 /: Series[Hypergeometric0F1[a_,fz_],{z_,aa_,nn_Integer}] :=
  1092.     hyper0F1[a, Series[fz,{z,aa,nn}], False (* regularized *)] /;
  1093.     (nn >= 0 && FreeQ[a,z])
  1094.  
  1095. Hypergeometric0F1 /: Hypergeometric0F1[a_,s_SeriesData] :=
  1096.     hyper0F1[a, s, False (* regularized *)] /;
  1097.     (simser[s] && FreeQ[a,s[[1]]])
  1098.  
  1099. Hypergeometric0F1Regularized /:
  1100. Series[Hypergeometric0F1Regularized[a_,fz_],{z_,aa_,nn_Integer}] :=
  1101.     hyper0F1[a, Series[fz,{z,aa,nn}], True (* regularized *)] /;
  1102.     (nn >= 0 && FreeQ[a,z])
  1103.  
  1104. Hypergeometric0F1Regularized /:
  1105. Literal[Hypergeometric0F1Regularized[a_, s_SeriesData]] :=
  1106.     hyper0F1[a, s, True (* regularized *)] /; (simser[s] && FreeQ[a,s[[1]]])
  1107.  
  1108.  (* --- Hypergeometric 1F1 --- *)
  1109.  
  1110.     (* special cases: expansion about Infinity *)
  1111.  
  1112. Hypergeometric1F1 /:
  1113. Literal[ Series[Hypergeometric1F1[a_, b_, z_], {z_, Infinity, nn_Integer}]] :=
  1114.     Series[Hypergeometric1F1Regularized[a, b, z], {z, Infinity, nn}] Gamma[b] /;
  1115.     (nn >= 0 && FreeQ[a, z] && FreeQ[b, z])
  1116.  
  1117. Hypergeometric1F1Regularized /:
  1118. Series[Hypergeometric1F1Regularized[a_, b_, z_], {z_, Infinity, nn_Integer}] :=
  1119.     Module[{s1, s2, onSesss},
  1120.     (* AS 13.5.1 *)
  1121.     onSesss = OnQ[Series::esss];
  1122.     If[onSesss, Off[Series::esss]];
  1123.     s1 = Series[hyper2F0r[a, 1+a-b, -z, nn], {z, Infinity, nn}];
  1124.     s2 = Series[hyper2F0r[b-a, 1-a, z, nn], {z, Infinity, nn}];
  1125.     s1 = E^(I Pi a) z^(-a) s1/Gamma[b-a] + E^z z^(a-b) s2/Gamma[a];
  1126.     If[onSesss, On[Series::esss]];
  1127.     s1
  1128.     ] /; (nn >= 0 && FreeQ[a, z] && FreeQ[b, z])
  1129.  
  1130.     (* generic case *)
  1131.  
  1132. hyper1F1[a_,b_,fz_,regflag_,Uflag_] :=
  1133.     Module[{i, n = 1+a-b, m = mterms[1,fz], f0 = faas[fz], z = fz[[1]]},
  1134.     If[Head[f0] === DirectedInfinity,
  1135.         If[Uflag,
  1136.         (* Abramowitz & Stegun 13.5.2 *)
  1137.         fz^(-a)Sum[Pochhammer[a,i]Pochhammer[n,i]/(i! (-fz)^i),{i,0,m}],
  1138.           (* else *)
  1139.         (* Abramowitz & Stegun 13.5.1 *)
  1140.         (Exp[I Pi a] HypergeometricU[a,b,fz]/Gamma[b-a] +
  1141.                 Exp[fz]/Gamma[a] HypergeometricU[b-a,b,-fz]) *
  1142.                 If[regflag, 1, Gamma[b]]
  1143.         ],
  1144.       (* else *)
  1145.         If[Uflag && TrueQ[f0==0],
  1146.         (* Abramowitz & Stegun 13.1.3 *)
  1147.         Pi/Sin[Pi b] (
  1148.             diffeq1F1[a,b,fz,f0,z,m,True,False]/Gamma[n] -
  1149.             diffeq1F1[n,2-b,fz,f0,z,m,True,False]/(Gamma[a] fz^(b-1))
  1150.             ),
  1151.           (* else *)
  1152.         diffeq1F1[a,b,fz,f0,z,m,regflag,Uflag]
  1153.         ]
  1154.         ]
  1155.     ];
  1156.  
  1157. diffeq1F1[a_,b_,fz_,aa_,z_,m_,regflag_,Uflag_] :=
  1158.     Module[{i, f0, f1, df, f},
  1159.     If[TrueQ[aa==0] && !Uflag,
  1160.         f = Sum[Pochhammer[a,i] z^i/(i! regden[b,i,regflag]), {i,0,m}];
  1161.         Return[ComposeSeries[f + O[z]^(m+1), fz]]
  1162.         ];
  1163.     f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i, {i,2,m}],{z,aa,m}];
  1164.     df = D[f,z];
  1165.     f = MapAt[Expand, (f /. Solve[z D[df,z] + (b-z) df - a f == 0,
  1166.                     Table[w[i],{i,2,m}]])[[1]], {3}];
  1167.     If[regflag,
  1168.         f = f /. {f0 -> Hypergeometric1F1Regularized[a,b,aa],
  1169.             f1 -> a Hypergeometric1F1Regularized[a+1,b+1,aa]},
  1170.         (* else *)
  1171.         f = f /. If[Uflag,
  1172.             {f0 -> HypergeometricU[a,b,aa],
  1173.             f1 -> -a HypergeometricU[a+1,b+1,aa]},
  1174.             (* else *)
  1175.             {f0 -> Hypergeometric1F1[a,b,aa],
  1176.             f1 -> a Hypergeometric1F1[a+1,b+1,aa]/b}]
  1177.         ];
  1178.     ComposeSeries[f, fz]
  1179.     ];
  1180.  
  1181. Hypergeometric1F1 /:
  1182.  Literal[ Series[Hypergeometric1F1[a_,b_,fz_],{z_,aa_,nn_Integer}]] :=
  1183.     Hypergeometric1F1[a,b,Series[fz,{z,aa,nn}]] /;
  1184.         (nn >= 0 && FreeQ[a,z] && FreeQ[b,z])
  1185.  
  1186. Hypergeometric1F1 /:
  1187. Literal[Hypergeometric1F1[a_,b_,s_SeriesData]] :=
  1188.     hyper1F1[a, b, s, False (* reg *), False (* U *)] /;
  1189.     (simser[s] && FreeQ[a,s[[1]]] && FreeQ[b,s[[1]]])
  1190.  
  1191. Hypergeometric1F1Regularized /:
  1192. Series[Hypergeometric1F1Regularized[a_,b_,fz_],{z_,aa_,nn_Integer}] :=
  1193.     Hypergeometric1F1Regularized[a,b,Series[fz,{z,aa,nn}]] /;
  1194.     (nn >= 0 && FreeQ[a,z] && FreeQ[b,z])
  1195.  
  1196. Hypergeometric1F1Regularized /:
  1197. Literal[Hypergeometric1F1Regularized[a_,b_,s_SeriesData]] :=
  1198.     Module[{m=1-Round[b]},
  1199.     Pochhammer[a,m] s^m * hyper1F1[a+m,m+1,s,True (* reg *),False (* U *)]
  1200.     ]/; (simser[s] && FreeQ[a,s[[1]]] && intQ[b] && b <= 0)
  1201.  
  1202. Hypergeometric1F1Regularized /:
  1203. Literal[Hypergeometric1F1Regularized[a_,b_,s_SeriesData]] :=
  1204.     hyper1F1[a, b, s, True (* reg *), False (* U *)] /;
  1205.     (simser[s] && FreeQ[a,s[[1]]] && FreeQ[b,s[[1]]])
  1206.  
  1207.  (* --- Hypergeometric 2F1 --- *)
  1208.  
  1209. hyper2F1[a_,b_,c_,fz_,regflag_] :=
  1210.     Module[{m=b-a, ab, l, f0 = faas[fz], nn = mterms[1,fz], z = fz[[1]]},
  1211.     Which[
  1212.         Head[f0] === DirectedInfinity,
  1213.         If[intQ[m],
  1214.             m = Round[m];
  1215.             ab = If[m > 0, b, a];
  1216.             m = Abs[m];
  1217.             l = c - a - m;
  1218.             If[intQ[l], l = Round[l]];
  1219.             If[IntegerQ[l] && l>=0,
  1220.             (* Erdelyi 2.1.4(19) p63 *)
  1221.             erdelyi21419[ab,m,l,fz,z,nn],
  1222.               (* else *)
  1223.             (* Abramowitz & Stegun 15.3.14 *)
  1224.             as15314[ab,m,c,fz,z,nn]
  1225.             ],
  1226.           (* else *)
  1227.             (* Abramowitz & Stegun 15.3.7 *)
  1228.             as1537[a,b,c,fz,z]
  1229.             ] * If[regflag, 1, Gamma[c]],
  1230.         TrueQ[f0==1],
  1231.         m=c-a-b;
  1232.         If[intQ[m],
  1233.             m = Round[m];
  1234.             If[m >= 0,
  1235.             (* Abramowitz & Stegun 15.3.11 *)
  1236.             as15311[a,b,m,fz,z,nn],
  1237.               (* else *)
  1238.             (* Abramowitz & Stegun 15.3.12 *)
  1239.             as15312[a,b,-m,fz,z,nn]
  1240.             ],
  1241.           (* else *)
  1242.             (* Abramowitz & Stegun 15.3.6 *)
  1243.             as1536[a,b,c,fz,z]
  1244.             ] * If[regflag, 1, Gamma[c]],
  1245.         TrueQ[f0==0],
  1246.         ComposeSeries[Sum[Pochhammer[a,l] Pochhammer[b,l]
  1247.             z^l/(l! regden[c,l,regflag]), {l,0,nn}] + O[z]^(nn+1),fz],
  1248.         True,
  1249.         diffeq2F1[a,b,c,fz,f0,z,nn,regflag]
  1250.         ]
  1251.     ];
  1252.  
  1253. diffeq2F1[a_,b_,c_,fz_,aa_,z_,m_,regflag_] :=
  1254.     Module[{i, f0, f1, df, f},
  1255.     f = Series[f0+f1(z-aa)+Sum[w[i](z-aa)^i, {i,2,m}],{z,aa,m}];
  1256.     df = D[f,z];
  1257.     f = MapAt[Expand, (f /. Solve[z(1-z)D[df,z] + (c-(a+b+1)z) df
  1258.             - a b f == 0, Table[w[i],{i,2,m}]])[[1]], {3}];
  1259.     If[regflag,
  1260.         f = f /. {f0 -> Hypergeometric2F1Regularized[a,b,c,aa],
  1261.             f1 -> a b Hypergeometric2F1Regularized[a+1,b+1,c+1,aa]},
  1262.       (* else *)
  1263.         f = f /. {f0 -> Hypergeometric2F1[a,b,c,aa],
  1264.             f1 -> a b Hypergeometric2F1[a+1,b+1,c+1,aa]/c}
  1265.         ];
  1266.     ComposeSeries[f, fz]
  1267.     ]
  1268.  
  1269. as1537[a_,b_,c_,fz_,z_] :=
  1270.     Gamma[b-a]/(Gamma[b] Gamma[c-a]) (-fz)^-a *
  1271.             Hypergeometric2F1[a,a+1-c,a+1-b,1/fz] +
  1272.         Gamma[a-b]/(Gamma[a] Gamma[c-b]) (-fz)^-b *
  1273.             Hypergeometric2F1[b,b+1-c,b+1-a,1/fz];
  1274.  
  1275. as15314[a_,m_,c_,fz_,z_,nn_] :=
  1276.     Module[{n},
  1277.     1/Gamma[a+m] *
  1278.     (1/Gamma[c-a] (-fz)^(-a-m) *
  1279.         Sum[Pochhammer[a,n+m] Pochhammer[1-c+a,n+m]/(n! (n+m)! fz^n)
  1280.             ((Log[-fz] /. logrule)+PolyGamma[1+m+n]+PolyGamma[1+n]-
  1281.             PolyGamma[a+m+n]-PolyGamma[c-a-m-n]),{n,0,nn-Re[a]-m}]+
  1282.         (-fz)^(-a) Sum[Gamma[m-n] Pochhammer[a,n]/
  1283.             (n! Gamma[c-a-n] fz^n), {n,0,Min[m-1,nn-Re[a]]}])
  1284.         ];
  1285.  
  1286. erdelyi21419[a_,m_,l_,fz_,z_,nn_] :=
  1287.     Module[{n},
  1288.     1/Gamma[a+m] * ((-fz)^(-a)((-fz)^(-m) *
  1289.         ((-1)^(m+l) Sum[Pochhammer[a,n+m](n-l)!/
  1290.                 (n! (n+m)! fz^n),{n,l,nn-Re[a]-m}] +
  1291.         1/(l+m-1)! Sum[Pochhammer[n+m] Pochhammer[1-m-l,n+m]/
  1292.             (n! (n+m)! fz^n)((Log[-fz]/.logrule)+PolyGamma[1+m+n]+
  1293.             PolyGamma[1+n] - PolyGamma[a+m+n] - PolyGamma[l-n]),
  1294.             {n,0,Min[l-1,nn-Re[a]-m]}]) +
  1295.         Sum[(m-n-1)! Pochhammer[a,n]/(n! (m+l-n-1)! fz^n),
  1296.             {n,0,Min[m-1,nn-Re[a]]}]))
  1297.         ];
  1298.  
  1299. as1536[a_,b_,c_,fz_,z_] :=
  1300.     Gamma[c-a-b]/(Gamma[c-a] Gamma[c-b]) *
  1301.             Hypergeometric2F1[a,b,a+b+1-c,1-fz] + 
  1302.         Gamma[a+b-c]/(Gamma[a] Gamma[b]) (1-fz)^(c-a-b) *
  1303.             Hypergeometric2F1[c-a,c-b,c+1-a-b,1-fz];
  1304.  
  1305. as15311[a_,b_,m_,fz_,z_,nn_] :=
  1306.     Module[{n,hyper=0},
  1307.     If[m!=0, hyper = Gamma[m]/(Gamma[a+m] Gamma[b+m]) *
  1308.             Sum[Pochhammer[a,n] Pochhammer[b,n] (1-fz)^n/
  1309.             (n! Pochhammer[1-m,n]), {n,0,Min[m-1,nn]}]];
  1310.     hyper - (fz-1)^m/(Gamma[a] Gamma[b]) *
  1311.             Sum[Pochhammer[a+m,n] Pochhammer[b+m,n]/(n! (n+m)!)
  1312.                 (1-fz)^n ((Log[1-fz] /. logrule)-PolyGamma[n+1]-
  1313.                 PolyGamma[n+m+1] + PolyGamma[a+n+m] +
  1314.                 PolyGamma[b+n+m]), {n,0,nn}]
  1315.         ];
  1316.  
  1317. as15312[a_,b_,m_,fz_,z_,nn_] :=
  1318.     Module[{n,hyper=0},
  1319.     If[m!=0, hyper = Gamma[m] (1-fz)^(-m)/(Gamma[a] Gamma[b]) *
  1320.             Sum[Pochhammer[a-m,n] Pochhammer[b-m,n] (1-fz)^n/
  1321.                 (n! Pochhammer[1-m,n]), {n,0,Min[m-1,nn]}]];
  1322.     hyper - (-1)^m/(Gamma[a-m] Gamma[b-m]) *
  1323.             Sum[Pochhammer[a,n] Pochhammer[b,n]/(n! (n+m)!)
  1324.                 (1-fz)^n ((Log[1-fz]/.logrule)-PolyGamma[n+1]-
  1325.                 PolyGamma[n+m+1] + PolyGamma[a+n] +
  1326.                 PolyGamma[b+n]), {n,0,nn}]
  1327.         ];
  1328.  
  1329. Hypergeometric2F1 /:
  1330. Literal[ Series[Hypergeometric2F1[a_,b_,c_,fz_],{z_,aa_,nn_Integer}]] :=
  1331.     Hypergeometric2F1[a,b,c,Series[fz,{z,aa,nn}]] /;
  1332.         (nn >= 0 && FreeQ[a,z] && FreeQ[b,z] && FreeQ[c,z])
  1333.  
  1334. Hypergeometric2F1 /:
  1335. Literal[Hypergeometric2F1[a_,b_,c_,fz_SeriesData]] :=
  1336.     hyper2F1[a,b,c,fz,False] /;
  1337.     (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]] && FreeQ[c,fz[[1]]])
  1338.  
  1339. Hypergeometric2F1Regularized /:
  1340. Series[Hypergeometric2F1Regularized[a_,b_,c_,fz_],{z_,aa_,nn_Integer}] :=
  1341.     Hypergeometric2F1Regularized[a,b,c,Series[fz,{z,aa,nn}]] /;
  1342.     (nn >= 0 && FreeQ[a,z] && FreeQ[b,z] && FreeQ[c,z])
  1343.  
  1344. Hypergeometric2F1Regularized /:
  1345. Literal[Hypergeometric2F1Regularized[a_,b_,c_,fz_SeriesData]] :=
  1346.     Module[{m=1-Round[c]},
  1347.     Pochhammer[a,m] Pochhammer[b,m] fz^m hyper2F1[a+m,b+m,m+1,fz,True]
  1348.     ]/; (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]] &&
  1349.         intQ[c] && c <= 0)
  1350.  
  1351. Hypergeometric2F1Regularized /:
  1352. Literal[Hypergeometric2F1Regularized[a_,b_,c_,fz_SeriesData]] :=
  1353.     hyper2F1[a,b,c,fz,True] /;
  1354.     (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]] && FreeQ[c,fz[[1]]])
  1355.  
  1356.  (* --- Hypergeometric U --- *)
  1357.  
  1358.     (* special case: expansion about Infinity *)
  1359.  
  1360. HypergeometricU /:
  1361. Series[HypergeometricU[a_,b_,z_],{z_, Infinity, nn_Integer}] :=
  1362.     Module[{n, onSesss}, (* AS 13.5.2 *)
  1363.     onSesss = OnQ[Series::esss];
  1364.     If[onSesss, Off[Series::esss]];
  1365.     n = z^(-a) Series[hyper2F0r[a, 1+a-b, -z, nn], {z, Infinity, nn}];
  1366.     If[onSesss, On[Series::esss]];
  1367.     n
  1368.     ] /; (nn >= 0 && FreeQ[a, z] && FreeQ[b, z])
  1369.  
  1370.     (* generic cases *)
  1371.  
  1372. as1316[a_,b_,fz_,z_,nn_] :=
  1373.     Module[{k, n = Round[b]-1},
  1374.     (n-1)!/(Gamma[a] fz^n) *
  1375.         Sum[Pochhammer[a-n,k] fz^k/(k! Pochhammer[1-n,k]),
  1376.             {k,0,Min[n-1,nn+n]}] -
  1377.         (-1)^n/Gamma[a-n] *
  1378.         (Hypergeometric1F1Regularized[a,n+1,fz] (Log[fz]/.logrule) +
  1379.             Sum[Pochhammer[a,k] fz^k/(k! (n+k)!) (PolyGamma[a+k] -
  1380.                 PolyGamma[1+k]-PolyGamma[1+n+k]),{k,0,nn}])
  1381.         ];
  1382.  
  1383. HypergeometricU /: Series[HypergeometricU[a_,b_,fz_],{z_,aa_,nn_Integer}] :=
  1384.     HypergeometricU[a,b,Series[fz,{z,aa,nn}]] /;
  1385.         (nn >= 0 && FreeQ[a,z] && FreeQ[b,z])
  1386.  
  1387. HypergeometricU /:
  1388. Literal[HypergeometricU[a_, b_, fz_SeriesData]] :=
  1389.     Module[{k, n = a+1-b, nn = mterms[1,fz], f0 = faas[fz], z = fz[[1]]},
  1390.     If[intQ[b] && (f0 == 0),
  1391.         If[b < 0,
  1392.             (* Abramowitz & Stegun 13.1.29 *)
  1393.             fz^(1-b) as1316[n,2-b,fz,z,nn],
  1394.             (* else *)
  1395.             (* Abramowitz & Stegun 13.1.6 *)
  1396.             as1316[a,b,fz,z,nn]
  1397.             ],
  1398.         (* else *)
  1399.         hyper1F1[a, b, fz, False (* reg *), True (* U *)]
  1400.         ] /; (* exclude simple closed forms *)
  1401.             !((intQ[a] && a < .5) || (intQ[n] && n < .5)) &&
  1402.             esssQ[f0, HypergeometricU[a, b, fz]]
  1403.     ]/; (simser[fz] && FreeQ[a,fz[[1]]] && FreeQ[b,fz[[1]]])
  1404.  
  1405.  (* -------------------- Laguerre polynomials ------------------------ *)
  1406.  
  1407. LaguerreL /: Series[LaguerreL[n_,a_:0,fz_],{z_,aa_,nn_Integer}] :=
  1408.     LaguerreL[n,a,Series[fz,{z,aa,nn}]] /;
  1409.         (nn >= 0 && FreeQ[n,z] && FreeQ[a,z])
  1410.  
  1411. LaguerreL /: LaguerreL[n_,a_:0,s_SeriesData] :=
  1412.     Pochhammer[n+1,a] Hypergeometric1F1Regularized[-n,a+1,s] /;
  1413.     (simser[s] && FreeQ[n,s[[1]]] && FreeQ[a,s[[1]]])
  1414.  
  1415.  (* -------------- Legendre polynomials and functions ------------------ *)
  1416.  
  1417. legendreP[n_,m_,fz_,opt___Rule] :=
  1418.     If[SameQ[LegendreType /. {opt} /. Options[LegendreP], Real],
  1419.         (* Gradshteyn & Ryzhik 8.704(b) *)
  1420.         ((1+fz)/(1-fz))^(m/2) *
  1421.             Hypergeometric2F1Regularized[-n,n+1,1-m,(1-fz)/2],
  1422.         (* else *)
  1423.         (* Gradshteyn & Ryzhik 8.702 *)
  1424.         ((fz+1)/(fz-1))^(m/2) *
  1425.             Hypergeometric2F1Regularized[-n,n+1,1-m,(1-fz)/2]
  1426.         ];
  1427.  
  1428. legendreQ[n_,m_,fz_,opt___Rule] :=
  1429.     If[SameQ[LegendreType /. {opt} /. Options[LegendreQ], Real],
  1430.         If[TrueQ[m==Round[m]],
  1431.             (* Gradshteyn & Ryzhik 8.737.2, Erdelyi 3.4(14) *)
  1432.             Pi/(2 Sin[(m+n)Pi]) * (
  1433.                Cos[(m+n)Pi] legendreP[n,m,fz,LegendreType->Real] -
  1434.                 legendreP[n,m,-fz,LegendreType->Real]),
  1435.             (* else *)
  1436.             (* Gradshteyn & Ryzhik 8.705(b), Erdelyi 3.4(17) *)
  1437.             Pi/(2 Sin[m Pi]) * (
  1438.                 Cos[m Pi] legendreP[n,m,fz,LegendreType->Real] -
  1439.                 Gamma[m+n+1]/Gamma[n-m+1]
  1440.                     legendreP[n,-m,fz,LegendreType->Real])
  1441.             ],
  1442.         (* else *)
  1443.         (* AS 8.1.3, Gradshteyn & Ryzhik 8.703 *)
  1444.         Exp[m Pi I] Gamma[m+n+1] Gamma[1/2] 2^-(n+1) (fz^2-1)^(m/2) *
  1445.         Hypergeometric2F1Regularized[(m+n+2)/2,(m+n+1)/2,n+3/2,fz^-2]/
  1446.         fz^(m+n+1)
  1447.         ];
  1448.  
  1449. LegendreP /:
  1450. Series[Literal[LegendreP[n_,m_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
  1451.     LegendreP[n,m,Series[fz,{z,aa,nn}],opt] /;
  1452.         (nn >= 0 && FreeQ[n,z] && FreeQ[m,z])
  1453.  
  1454. LegendreP /: LegendreP[n_,m_,s_SeriesData,opt___Rule] :=
  1455.     legendreP[n,m,s,opt] /; (simser[s] && FreeQ[n,s[[1]]] && FreeQ[m,s[[1]]])
  1456.  
  1457. LegendreQ /:
  1458. Series[Literal[LegendreQ[n_,m_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
  1459.     LegendreQ[n,m,Series[fz,{z,aa,nn}],opt] /;
  1460.         (nn >= 0 && FreeQ[n,z] && FreeQ[m,z])
  1461.  
  1462. LegendreQ /: LegendreQ[n_,m_,s_SeriesData,opt___Rule] :=
  1463.     legendreQ[n,m,s,opt] /; (simser[s] && FreeQ[n,s[[1]]] && FreeQ[m,s[[1]]])
  1464.  
  1465. (* the following Legendre definitions are only needed because the pattern
  1466. matcher doesn't work right. the "m"'s above need to be ":0"'ed *)
  1467.  
  1468. LegendreP /:
  1469. Series[Literal[LegendreP[n_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
  1470.     LegendreP[n,0,Series[fz,{z,aa,nn}],opt] /;
  1471.     (nn >= 0 && FreeQ[n,z])
  1472.  
  1473. LegendreP /: LegendreP[n_,s_SeriesData,opt___Rule] :=
  1474.     legendreP[n,0,s,opt] /; (simser[s] && FreeQ[n,s[[1]]])
  1475.  
  1476. LegendreQ /:
  1477. Series[Literal[LegendreQ[n_,fz_,opt___Rule]],{z_,aa_,nn_Integer}] :=
  1478.     LegendreQ[n,0,Series[fz,{z,aa,nn}],opt] /;
  1479.     (nn >= 0 && FreeQ[n,z])
  1480.  
  1481. LegendreQ /: LegendreQ[n_,s_SeriesData,opt___Rule] :=
  1482.     legendreQ[n,0,s,opt] /; (simser[s] && FreeQ[n,s[[1]]])
  1483.  
  1484.  (* ----------------------- Lerch's transcendant ---------------------- *)
  1485.  
  1486. LerchPhi /:
  1487. Series[Literal[LerchPhi[fz_,s_,a_,opt___Rule]],{z_,aa_,nn_Integer}] :=
  1488.     LerchPhi[Series[fz,{z,aa,nn}],s,a,opt] /;
  1489.     (nn >= 0 && FreeQ[s,z] && FreeQ[a,z])
  1490.  
  1491. LerchPhi /:
  1492. Literal[LerchPhi[fz_SeriesData,s_,a_,opt___Rule]] :=
  1493.     Module[{k},
  1494.     Sum[If[TrueQ[a+k==0], 0, fz^k/((a+k)^2)^(s/2)],
  1495.         {k,0,mterms[1,fz]}] /;
  1496.     !(DoublyInfinite /. {opt} /. Options[LerchPhi]) && (faas[fz] == 0)
  1497.     ] /; (simser[fz] && FreeQ[s,fz[[1]]] && FreeQ[a,fz[[1]]])
  1498.  
  1499.  (* ------------------ Logarithmic integral ----------------------- *)
  1500.  
  1501.     (* special case: expansion about Infinity *)
  1502.  
  1503. LogIntegral /: Series[LogIntegral[z_], {z_, Infinity, nn_Integer}] :=
  1504.     Module[{j, s, onSesss},
  1505.     onSesss = OnQ[Series::esss];
  1506.     If[onSesss, Off[Series::esss]];
  1507.     s = Series[Sum[j! z^(-j-1), {j, 0, nn}], {z, Infinity, nn}];
  1508.     s = z (s /. z -> Log[z]);
  1509.     If[onSesss, On[Series::esss]];
  1510.     s
  1511.     ] /; nn >= 0
  1512.  
  1513.     (* generic case *)
  1514.  
  1515. LogIntegral /: Series[LogIntegral[fz_],{z_,aa_,nn_Integer}] :=
  1516.     LogIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0
  1517.  
  1518. LogIntegral /:
  1519. Literal[LogIntegral[s_SeriesData]] :=
  1520.     Module[{f0 =faas[s], i},
  1521.     ExpIntegralEi[Log[s]] /;
  1522.         !TrueQ[f0 == 0] && esssQ[f0, LogIntegral[s]]
  1523.     ] /; simser[s]
  1524.  
  1525.  (* ------------------------- Pochhammer --------------------------------- *)
  1526.  
  1527. Pochhammer/:
  1528. Series[Pochhammer[f1z_,f2z_],{z_,aa_,nn_Integer}] :=
  1529.     Module[{sf1 = Series[f1z,{z,aa,nn}],sf2 = Series[f2z,{z,aa,nn}]},
  1530.     Gamma[sf1+sf2] / Gamma[sf1]] /; nn >= 0
  1531.  
  1532. Pochhammer/:
  1533. Literal[Pochhammer[sf1_SeriesData, sf2_SeriesData]]:=
  1534.     Gamma[sf1+sf2] / Gamma[sf1] /; simser[sf1] && simser[sf1]
  1535.  
  1536. Pochhammer/:
  1537. Literal[Pochhammer[sf_SeriesData,a_]]:=
  1538.     Gamma[sf +
  1539.     Series[a,{sf[[1]],sf[[2]],Round[(sf[[5]]-Min[0,sf[[4]]])/sf[[6]]]}]] /
  1540.     Gamma[sf] /; simser[sf]
  1541.  
  1542. Pochhammer/:
  1543. Literal[Pochhammer[a_,sf_SeriesData]] :=
  1544.     Module[{sa},
  1545.     sa = Series[a, {sf[[1]], sf[[2]],
  1546.         Round[(sf[[5]] - Min[0, sf[[4]]])/sf[[6]]]}];
  1547.     Gamma[sf+sa] / Gamma[sa]] /; simser[sf]
  1548.  
  1549.  (* -------------------- PolyGamma functions -------------------------- *)
  1550.  
  1551.     (* special cases: expansion about Infinity *)
  1552.  
  1553. LogGamma /: Series[LogGamma[z_],{z_, Infinity, nn_Integer}] :=
  1554.     Module[{s, i, onSesss},
  1555.     onSesss = OnQ[Series::esss];
  1556.     If[onSesss, Off[Series::esss]];
  1557.     s = Sum[BernoulliB[i]/(i(i-1) z^(i-1)), {i, 2, nn, 2}];
  1558.     s = (z - 1/2) Log[z] - z + Log[2 Pi]/2 + Series[s, {z, Infinity, nn}];
  1559.     If[onSesss, On[Series::esss]];
  1560.     s
  1561.     ] /; nn >= 0
  1562.  
  1563. PolyGamma /: Series[PolyGamma[n_Integer:0, z_],{z_, Infinity, nn_Integer}] :=
  1564.     Module[{s, i},
  1565.     onSesss = OnQ[Series::esss];
  1566.     If[onSesss, Off[Series::esss]];
  1567.     s = Sum[BernoulliB[i]/(i z^i), {i, 2, nn, 2}];
  1568.     s = Log[z] - 1/(2z) - Series[s, {z, Infinity, nn}];
  1569.     s = D[s, {z, n}];
  1570.     If[onSesss, On[Series::esss]];
  1571.     s
  1572.     ] /; n >= 0 && nn >= 0
  1573.  
  1574.     (* generic cases *)
  1575.  
  1576. LogGamma /: Series[LogGamma[fz_],{z_,aa_,nn_Integer}] :=
  1577.     LogGamma[Series[fz,{z,aa,nn}]] /; nn >= 0
  1578.  
  1579. LogGamma /: LogGamma[s_SeriesData] :=
  1580.     Module[{f0 = faas[s]},
  1581.     (f0 = If[intQ[f0] && f0 < .5,
  1582.         (Round[f0]-1) Pi I,
  1583.         LogGamma[f0] - Log[Gamma[f0]]
  1584.         ];
  1585.     Log[Gamma[s]] + f0) /; esssQ[f0, LogGamma[s]]
  1586.     ] /; simser[s]
  1587.  
  1588. PolyGamma /: Series[PolyGamma[n_Integer:0,fz_],{z_,aa_,nn_Integer}] :=
  1589.     PolyGamma[n,Series[fz,{z,aa,nn+n}]] /; nn >= n
  1590.  
  1591. PolyGamma /: PolyGamma[-1,s_SeriesData] :=
  1592.     LogGamma[s] /; simser[s] && esssQ[faas[s], LogGamma[s]]
  1593.  
  1594. PolyGamma /: PolyGamma[n_Integer:0, s_SeriesData] :=
  1595.     Module[{a, pgs, f0 = faas[s], z = s[[1]]},
  1596.     (pgs = Gamma[Series[z, {z, f0, mterms[1, s]}]];
  1597.     pgs = D[D[pgs, z]/pgs, {z, n}];
  1598.     pgs = ComposeSeries[pgs, s] /. (a_ :> (Expand[a] /; FreeQ[a,z]))) /;
  1599.     Head[f0] =!= DirectedInfinity
  1600.     ] /; n >= 0 && simser[s]
  1601.  
  1602.  (* ------------------- PolyLogarithm functions ------------------------- *)
  1603.  
  1604. PolyLog /: Series[PolyLog[n_,fz_],{z_,aa_,nn_Integer}] :=
  1605.     PolyLog[n,Series[fz,{z,aa,nn}]] /; (nn >= 0  && FreeQ[n,z])
  1606.  
  1607. polylog1[n_,fz_,m_] :=
  1608.     Module[{k, fpow=1, plder, logpart=0, t, sum=0, i},
  1609.     If[TrueQ[n==1], Return[-Log[1-fz] //. logrule]];
  1610.     plder = PolyLog[n,t];
  1611.     Do[sum += (plder /. t-> 1)fpow/i!;
  1612.         If[logpart =!= 0,
  1613.             logpart = logpart /. t -> 1-t;
  1614.             logpart = Normal[Series[logpart,{t,0,m-i}]];
  1615.             Do[logpart = -Integrate[logpart,t],{i}];
  1616.             sum += (logpart /. t -> 1-fz)
  1617.             ];
  1618.         fpow *= fz-1;
  1619.         {logpart,plder} = logsplit[D[plder,t]],
  1620.         {i,0,m}];
  1621.     sum //. logrule
  1622.     ]
  1623.  
  1624. PolyLog /: Literal[PolyLog[n_,fz_SeriesData]] :=
  1625.     Module[{k, f0=faas[fz], m = mterms[1, fz]},
  1626.     Which[TrueQ[f0==0], Sum[fz^k/k^n,{k,m}],
  1627.         TrueQ[f0==1] && IntegerQ[n] && n>0, polylog1[n,fz,m],
  1628.         True, Sum[Derivative[0,k][PolyLog][n,f0](fz-f0)^k/(k!),{k,0,m}]
  1629.         ]
  1630.     ] /; simser[fz] && FreeQ[n, fz[[1]]]
  1631.  
  1632.  (* ------------------- Spherical Harmonic functions --------------------- *)
  1633.  
  1634. SphericalHarmonicY /:
  1635. Series[SphericalHarmonicY[l_,m_,ftheta_,phi_],{theta_,aa_,nn_Integer}] :=
  1636.     SphericalHarmonicY[l,m,Series[ftheta,{theta,aa,nn}],phi] /;
  1637.     (nn >= 0 && FreeQ[l,theta] && FreeQ[m,theta] && FreeQ[phi,theta])
  1638.  
  1639. SphericalHarmonicY /:
  1640. Literal[SphericalHarmonicY[l_,m_,s_SeriesData,phi_]] :=
  1641.     Exp[m I phi] Sqrt[(2l+1)/(4 Pi Pochhammer[l-m+1,2m])] *
  1642.         LegendreP[l,m,Cos[s],LegendreType->Real] /;
  1643.     (simser[s] && FreeQ[l,s[[1]]] && FreeQ[m,s[[1]]] && FreeQ[phi,s[[1]]])
  1644.  
  1645.  (* ----------------------- Zeta function ---------------------------- *)
  1646.  
  1647. $stltjesgammas;
  1648. $NStltjesGamma;
  1649. $PrecNStltjesGamma;
  1650.  
  1651. (* values for debugging:
  1652. $NStltjesGamma[1] = -0.0728158454836767248605863758749013191377363383343
  1653. $NStltjesGamma[2] = -0.00969036319287231848453038604
  1654. $NStltjesGamma[3] = 0.0020538344203033458661600465
  1655. $NStltjesGamma[4] = 0.0023253700654673000574681702
  1656. $NStltjesGamma[5] = 0.000793323817301062701753335
  1657. $NStltjesGamma[6] = -0.000238769345430199609872422
  1658. *)
  1659.  
  1660. StieltjesGamma[0] = EulerGamma;
  1661.  
  1662. N[StieltjesGamma[n_Integer]] ^:= N[StieltjesGamma[n], $MachinePrecision]
  1663.  
  1664. N[StieltjesGamma[n_Integer], m_Integer] ^:=
  1665.     Module[{eglist, i},
  1666.     If[NumberQ[$NStltjesGamma[n]] && m <= $PrecNStltjesGamma[n],
  1667.         N[$NStltjesGamma[n], m],
  1668.       (* else *)
  1669.         eglist = EGList[n, m];
  1670.         AbortProtect[
  1671.         Do[ If[!NumberQ[$NStltjesGamma[i]] || m > $PrecNStltjesGamma[i],
  1672.             {$NStltjesGamma[i], $PrecNStltjesGamma[i]} =
  1673.             {eglist[[i]], m}], {i, n}]];
  1674.         $NStltjesGamma[n]]] /; (n > 0 && m > 0)
  1675.  
  1676. EGList[n_, m_] :=
  1677.     Module[{zser, i},
  1678.     zser = ($stltjesgammas - 1) N[EulerGamma, m] -
  1679.         Sum[$NSigmaZeta[i, m] (1 - $stltjesgammas)^i, {i, n+1, 2, -1}];
  1680.     zser = Exp[Series[zser, {$stltjesgammas, 1, n+1}]];
  1681.     zser = Drop[zser[[3]], 2];
  1682.     zser Table[(-1)^i i!, {i, n}]];
  1683.  
  1684. $NSigmaZeta;
  1685. $PrecNSigmaZeta;
  1686.  
  1687. (* values for debugging:
  1688. $NSigmaZeta[2] = 0.0937731164201826122986016923027207941919722317905
  1689. $NSigmaZeta[3] = 0.017229544011064297934002741028
  1690. $NSigmaZeta[4] = 0.00368791470636343601614505920
  1691. $NSigmaZeta[5] = 0.000904895577699075748249223220
  1692. $NSigmaZeta[6] = 0.000241132534087530523369496737
  1693. $NSigmaZeta[7] = 0.000067363439740772150048502904
  1694. *)
  1695.  
  1696. $NSigmaZeta[n_, m_] :=
  1697.     Module[{szlist, i},
  1698.     (* coefficients in the series expansion of Log[(s-1) Zeta[s]]
  1699.         about s == 1 *)
  1700.     If[NumberQ[$NSigmaZeta[n]] && m <= $PrecNSigmaZeta[n],
  1701.         Return[N[$NSigmaZeta[n], m]]];
  1702.     szlist = SZList[n, m];
  1703.     AbortProtect[
  1704.         Do[ If[!NumberQ[$NSigmaZeta[i]] || m > $PrecNSigmaZeta[i],
  1705.         {$NSigmaZeta[i], $PrecNSigmaZeta[i]} = {szlist[[i-1]], m}],
  1706.         {i, 2, n}]];
  1707.     $NSigmaZeta[n]]
  1708.  
  1709. $Nzeta32;
  1710. $PrecNzeta32;
  1711.  
  1712. (* values for debugging:
  1713. $Nzeta32[2] = 0.1168502750680849136771556874922594459571062129525
  1714. $Nzeta32[3] = 0.017266596754881666574923630441
  1715. $Nzeta32[4] = 0.00366950790104801363656336638
  1716. $Nzeta32[5] = 0.000904752559027923226702063001
  1717. $Nzeta32[6] = 0.000241179440157020317746431190
  1718. $Nzeta32[7] = 0.000067364093196650679301661642
  1719. *)
  1720.  
  1721. $Nzeta32[i_, m_] :=
  1722.     Module[{zeta},
  1723.     If[NumberQ[$Nzeta32[j]] && m <= $PrecNzeta32[j],
  1724.         Return[N[$Nzeta32[j], m]]];
  1725.     zeta = N[Zeta[i, 3/2], m]/(i 2^i);
  1726.     AbortProtect[{$Nzeta32[i], $PrecNzeta32[i]} = {zeta, m}];
  1727.     zeta];
  1728.  
  1729. SZList[n_, m_] :=
  1730.     Module[{xiser, i},
  1731.     (* series expansion of 2*xi(s) about s == 1 *)
  1732.     xiser = 1 + Sum[($NXiBeta[i-2, m] + $NXiBeta[i-1, m]) *
  1733.         ($stltjesgammas - 1)^i, {i, n}];
  1734.     (* series expansion of log(2*xi(s)) about s == 1 *)
  1735.     xiser = Log[Series[xiser, {$stltjesgammas, 1, n}]];
  1736.     xiser = Rest[xiser[[3]]] Table[(-1)^i, {i, n-1}];
  1737.     (* coefficients in the series expansion of Log[(s-1) Zeta[s]]
  1738.         about s == 1 *)
  1739.     Table[$Nzeta32[i, m], {i, 2, n}] + xiser];
  1740.  
  1741. $NXiBeta;
  1742. $PrecNXiBeta;
  1743.  
  1744. $NXiBeta[-1, _] := 0;
  1745. $NXiBeta[0, m_] := N[1 + EulerGamma/2 - Log[2 Sqrt[Pi]], m];
  1746.  
  1747. (* values for debugging:
  1748. $NXiBeta[1] = 0.000248155568105149320572108979315189423336753288730593
  1749. $NXiBeta[2] = 0.0002498282818177993517790627950742989298
  1750. $NXiBeta[3] = 3.3534484987276532770582890420745483*10^-6
  1751. $NXiBeta[4] = 1.69680629349152089252694618293863650*10^-6
  1752. $NXiBeta[5] = 2.418074837001468525320879096545763*10^-8
  1753. $NXiBeta[6] = 8.197666248796084350271868231276427*10^-9
  1754. *)
  1755.  
  1756. betatheta[t_, m_] :=
  1757.     Module[{sum, term, cterm, dterm, eps = 0.1^m},
  1758.     (* evaluate the theta function Sum[E^(-n^2 Pi t), {n, 1, Infinity}] *)
  1759.     sum = term = Exp[-N[Pi, m] t];
  1760.     cterm = term;
  1761.     dterm = cterm^2;
  1762.     While[term > sum eps,
  1763.         cterm *= dterm;
  1764.         term *= cterm;
  1765.         sum += term
  1766.         ];
  1767.     sum
  1768.     ];
  1769.  
  1770. $NXiBeta[j_, m_] :=
  1771.     Module[{prec = m+10, t, beta, onNIslwcon},
  1772.     If[NumberQ[$NXiBeta[j]] && m <= $PrecNXiBeta[j],
  1773.         Return[N[$NXiBeta[j], m]]];
  1774.     onNIslwcon = OnQ[NIntegrate::slwcon];
  1775.     If[onNIslwcon, Off[NIntegrate::slwcon]];
  1776.     beta = NIntegrate[betatheta[t, prec] (Log[t]/2)^j (Sqrt[t] + (-1)^j)/t,
  1777.         {t, 1, Infinity}, Method -> DoubleExponential,
  1778.         WorkingPrecision -> prec, PrecisionGoal -> m,
  1779.         MaxRecursion -> 20]/j!;
  1780.     If[onNIslwcon, On[NIntegrate::slwcon]];
  1781.     AbortProtect[{$NXiBeta[j], $PrecNXiBeta[j]} = {beta, m}];
  1782.     beta];
  1783.  
  1784.  
  1785. Zeta /: Series[Zeta[fz_],{z_,aa_,nn_Integer}] :=
  1786.     Zeta[Series[fz,{z,aa,nn}]] /; nn >= 0
  1787.  
  1788. Zeta /: Literal[Zeta[s_SeriesData]] :=
  1789.     Module[{j, k, f, f0 = faas[s], m = mterms[1,s], z = s[[1]]},
  1790.     If[TrueQ[f0==1],
  1791.         f = 1/z+EulerGamma+Sum[(-1)^j StieltjesGamma[j] z^j/j!,{j,m}] +
  1792.             O[z]^(m+1),
  1793.         (* else *)
  1794.         f = Sum[Derivative[k][Zeta][f0]z^k/k!, {k,0,m}] + O[z]^(m+1)
  1795.         ];
  1796.     f[[2]] = f0;
  1797.     ComposeSeries[f, s]
  1798.     ] /; simser[s]
  1799.  
  1800.  (* ------------------- Various Integral functions ----------------------- *)
  1801.  
  1802. SinIntegral /: Series[SinIntegral[fz_],{z_,aa_,nn_Integer}] :=
  1803.     SinIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0
  1804.  
  1805. SinIntegral /: Literal[SinIntegral[fz_SeriesData]] :=
  1806.    Module[{temp, f, f0 = faas[fz]},
  1807.     f = Integrate[Series[Sin[temp+f0]/(temp+f0),{temp,0,mterms[1,fz]}],
  1808.         temp];
  1809.     SinIntegral[f0] + ComposeSeries[f, fz-f0]
  1810.     ] /; simser[fz]
  1811.  
  1812. CosIntegral /: Series[CosIntegral[fz_],{z_,aa_,nn_Integer}] :=
  1813.     CosIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0
  1814.  
  1815. CosIntegral /: Literal[CosIntegral[fz_SeriesData]] :=
  1816.     Module[{i, temp, f, m = mterms[1,fz], f0 = faas[fz]},
  1817.     If[TrueQ[f0 == 0],
  1818.         f = Sum[(-1)^i temp^(2i)/(2i(2i)!),{i,1,m/2}] + O[temp]^(m+1);
  1819.         ComposeSeries[f, fz] + (Log[fz] /. logrule) + EulerGamma,
  1820.         (* else *)
  1821.         f = Integrate[Series[Cos[temp+f0]/(temp+f0),{temp,0,m}], temp];
  1822.         CosIntegral[f0] + ComposeSeries[f, fz-f0]
  1823.         ]
  1824.     ] /; simser[fz]
  1825.  
  1826. SinhIntegral /: Series[SinhIntegral[fz_],{z_,aa_,nn_Integer}] :=
  1827.     SinhIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0
  1828.  
  1829. SinhIntegral /:
  1830. Literal[SinhIntegral[fz_SeriesData]] :=
  1831.    Module[{temp, f, f0=faas[fz]},
  1832.     f = Integrate[Series[Sinh[temp+f0]/(temp+f0),{temp,0,mterms[1,fz]}],
  1833.         temp];
  1834.     SinhIntegral[f0] + ComposeSeries[f, fz-f0]
  1835.     ] /; simser[fz]
  1836.  
  1837. CoshIntegral /: Series[CoshIntegral[fz_],{z_,aa_,nn_Integer}] :=
  1838.     CoshIntegral[Series[fz,{z,aa,nn}]] /; nn >= 0
  1839.  
  1840. CoshIntegral /:
  1841. Literal[CoshIntegral[fz_SeriesData]] :=
  1842.     Module[{i, temp, f, m = mterms[1,fz], f0 = faas[fz]},
  1843.     If[TrueQ[f0 == 0],
  1844.         f = Sum[temp^(2i)/(2i(2i)!),{i,1,m/2}] + O[temp]^(m+1);
  1845.         ComposeSeries[f, fz] + (Log[fz] /. logrule) + EulerGamma,
  1846.         (* else *)
  1847.         f = Integrate[Series[Cosh[temp+f0]/(temp+f0),{temp,0,m}], temp];
  1848.         CoshIntegral[f0] + ComposeSeries[f, fz-f0]
  1849.         ]
  1850.     ] /; simser[fz]
  1851.  
  1852. FresnelC/: Series[FresnelC[fz_],{z_,aa_,nn_Integer}] :=
  1853.     FresnelC[Series[fz,{z,aa,nn}]] /; nn >= 0
  1854.  
  1855. FresnelC/: Literal[FresnelC[fz_SeriesData]] :=
  1856.    Module[{temp, f, f0=faas[fz]},
  1857.     f = Integrate[Series[Cos[Pi (temp+f0)^2/2], {temp,0,mterms[1,fz]}],
  1858.         temp];
  1859.     FresnelC[f0] + ComposeSeries[f, fz-f0]
  1860.     ] /; simser[fz]
  1861.  
  1862. FresnelS/: Series[FresnelS[fz_],{z_,aa_,nn_Integer}] :=
  1863.     FresnelS[Series[fz,{z,aa,nn}]] /; nn >= 0
  1864.  
  1865. FresnelS/: Literal[FresnelS[fz_SeriesData]] :=
  1866.    Module[{temp, f, f0=faas[fz]},
  1867.     f = Integrate[Series[Sin[Pi (temp+f0)^2/2], {temp,0,mterms[1,fz]}],
  1868.         temp];
  1869.     FresnelS[f0] + ComposeSeries[f, fz-f0]
  1870.     ] /; simser[fz]
  1871.  
  1872.  (* ----------------------- Sums and Integrals ------------------------- *)
  1873.  
  1874. Series[Literal[Integrate[g_,{t_,f1z_,f2z_}]],{z_,aa_,nn_Integer}] :=
  1875.     Module[{sf1 = Series[f1z,{z,aa,nn}],sf2 = Series[f2z,{z,aa,nn}],
  1876.         f01,f02,f1,f2,temp},
  1877.     (If[Head[sf1] === SeriesData,
  1878.         f01 = faas[sf1];
  1879.         f1 = Integrate[Series[(g /. t -> temp+f01),
  1880.             {temp,0,mterms[1,sf1[[4]],sf1[[5]],sf1[[6]]]}],temp];
  1881.         f1 = ComposeSeries[f1,sf1-f01],
  1882.         (* else *)
  1883.         f01 = sf1;
  1884.         f1 = 0
  1885.         ];
  1886.     If[Head[sf2] === SeriesData,
  1887.         f02 = faas[sf2];
  1888.         f2 = Integrate[Series[(g /. t -> temp+f02),
  1889.             {temp,0,mterms[1,sf2[[4]],sf2[[5]],sf2[[6]]]}],temp];
  1890.         f2 = ComposeSeries[f2,sf2-f02],
  1891.         (* else *)
  1892.         f02 = sf2;
  1893.         f2 = 0
  1894.         ];
  1895.     Integrate[g,{t,f01,f02}] + ((f2-f1) /. temp -> z)) /;
  1896.         ((Head[sf1] =!= Series) && (Head[sf2] =!= Series))
  1897.     ] /; (nn >= 0 && FreeQ[g,z] && FreeQ[t,z])
  1898.  
  1899. Series[Literal[Sum[g_,iter_List]],{z_,aa_,nn_Integer}] :=
  1900.     Sum[Series[g,{z,aa,nn}],iter] /; (nn >= 0 && FreeQ[iter,z])
  1901.  
  1902. If[onSDsdatc, On[SeriesData::sdatc]];
  1903. Remove[onSDsdatc];
  1904.  
  1905. End[] (* "SpecialFunctions`Series`Private`" *)
  1906.  
  1907. SetAttributes[AiryAi, ReadProtected];
  1908. SetAttributes[AiryAiPrime, ReadProtected];
  1909. SetAttributes[AiryBi, ReadProtected];
  1910. SetAttributes[AiryBiPrime, ReadProtected];
  1911. SetAttributes[BesselI, ReadProtected];
  1912. SetAttributes[BesselJ, ReadProtected];
  1913. SetAttributes[BesselK, ReadProtected];
  1914. SetAttributes[BesselY, ReadProtected];
  1915. SetAttributes[Beta, ReadProtected];
  1916. SetAttributes[BetaRegularized, ReadProtected];
  1917. SetAttributes[ChebyshevT, ReadProtected];
  1918. SetAttributes[ChebyshevU, ReadProtected];
  1919. SetAttributes[ExpIntegralE, ReadProtected];
  1920. SetAttributes[ExpIntegralEi, ReadProtected];
  1921. SetAttributes[Factorial, ReadProtected];
  1922. SetAttributes[FresnelC, ReadProtected];
  1923. SetAttributes[FresnelS, ReadProtected];
  1924. SetAttributes[Gamma, ReadProtected];
  1925. SetAttributes[GammaRegularized, ReadProtected];
  1926. SetAttributes[GegenbauerC, ReadProtected];
  1927. SetAttributes[HermiteH, ReadProtected];
  1928. SetAttributes[Hypergeometric0F1, ReadProtected];
  1929. SetAttributes[Hypergeometric1F1, ReadProtected];
  1930. SetAttributes[Hypergeometric2F1, ReadProtected];
  1931. SetAttributes[Hypergeometric0F1Regularized, ReadProtected];
  1932. SetAttributes[Hypergeometric1F1Regularized, ReadProtected];
  1933. SetAttributes[Hypergeometric2F1Regularized, ReadProtected];
  1934. SetAttributes[HypergeometricU, ReadProtected];
  1935. SetAttributes[LaguerreL, ReadProtected];
  1936. SetAttributes[LegendreP, ReadProtected];
  1937. SetAttributes[LegendreQ, ReadProtected];
  1938. SetAttributes[LerchPhi, ReadProtected];
  1939. SetAttributes[LogIntegral, ReadProtected];
  1940. SetAttributes[LogGamma, ReadProtected];
  1941. SetAttributes[Pochhammer, ReadProtected];
  1942. SetAttributes[PolyGamma, ReadProtected];
  1943. SetAttributes[PolyLog, ReadProtected];
  1944. SetAttributes[SphericalHarmonicY, ReadProtected];
  1945. SetAttributes[Zeta, ReadProtected];
  1946. SetAttributes[SinIntegral, ReadProtected];
  1947. SetAttributes[CosIntegral, ReadProtected];
  1948. SetAttributes[SinhIntegral, ReadProtected];
  1949. SetAttributes[CoshIntegral, ReadProtected];
  1950. SetAttributes[SeriesData, ReadProtected];
  1951. SetAttributes[Series, ReadProtected];
  1952. SetAttributes[StieltjesGamma, ReadProtected];
  1953.  
  1954. Protect[AiryAi,AiryAiPrime,AiryBi,AiryBiPrime,BesselI,BesselJ,BesselK,BesselY,
  1955. Beta,BetaRegularized, ChebyshevT,ChebyshevU,ExpIntegralE,ExpIntegralEi,
  1956. Factorial,FresnelC,FresnelS,Gamma,GammaRegularized,GegenbauerC,HermiteH,
  1957. Hypergeometric0F1,Hypergeometric0F1Regularized,Hypergeometric1F1,
  1958. Hypergeometric1F1Regularized,Hypergeometric2F1,Hypergeometric2F1Regularized,
  1959. HypergeometricU,LaguerreL,LegendreP,LegendreQ,LerchPhi,LogIntegral,LogGamma,
  1960. Pochhammer,PolyGamma,PolyLog,SphericalHarmonicY,Zeta,SinIntegral,CosIntegral,
  1961. SinhIntegral,CoshIntegral,SeriesData,Series,StieltjesGamma];
  1962.  
  1963. End[]  (* "System`" *)
  1964.  
  1965.