numerics - What happened to SequenceLimit?

In older versions of Mathematica, there was a function called SequenceLimit that allowed taking the limit of a numerical sequence. It is useful for speeding up the convergence of numerical algorithms, and sometimes for getting them unstuck using the last few results. What happened to it? Any improvements available?

Answer

SequenceLimit became NumericalMath`NSequenceLimit. In the past, that change broke some of the examples on NLimit that used SequenceLimit. I reported the breakage, and the examples were fixed to use the same syntax but not actually call the now-non-existent function. The only remaining reference in the documentation that I can find to the use of SequenceLimit as if it were still a function is on the last bullet point of NLimit's detailed description for its Method option: "uses SequenceLimit on constructed sequence".

Recently, I started investigating the behavior of NumericalMath`NSequenceLimit because I need it in my engineering work. I feel NumericalMath`NSequenceLimit has some quite confusing shortcomings in regard to the approximation it returns vs. the number of terms supplied, so I wrote my own version below from some references. I also posted a C++ gist version on GitHub.

A comparison appears below showing the convergence of partial sums of series approximations to Pi and E. Afterward, I use symbolic lists to document several parts of NumericalMath`NSequenceLimit that seem... possibly wrong... to me.

sPi@n_:=Sum[4(-1)^nn/(2nn+1),{nn,0,n}];
Table[sPi@n,{n,0,4-1}]//N
(*{4.,2.66667,3.46667,2.89524}*)

NumericalMath`NSequenceLimit@%
(*3.16667*)

sequenceLimit@%%
(*3.13333*)    
(*my answer is closer to Pi~=3.141, for good reason as shown

  below in the general case with 4 terms {a,b,c,d}*)

sE@n_:=Piecewise[{{(1+1/n)^n ,n!=0}},1] 
Table[sE@n,{n,0,4-1}]//N
(*{1.,2.,2.25,2.37037}*)

NumericalMath`NSequenceLimit@%
(*2.33333*)

sequenceLimit@%%

(*2.48214*)    
(*my answer is closer to E~=2.71828, for good reason as shown
  below in the general case with 4 terms {a,b,c,d}*)

Below, I illustrate the difference between the behavior of the built-in NumericalMath`NSequenceLimit function vs. my sequenceLimit function for symbolic lists, like {a,b,c,d}.

Quiet@NumericalMath`NSequenceLimit[{a,b,c}]
(*NumericalMath`NSequenceLimit[{a,b,c}]*)
(*yes, unevaluated with 3 inputs*)

sequenceLimit[{a,b,c}]//Simplify (*my function*)

(*(-b^2+a*c)/(a-2*b+c)*)
(*same as Mathematica's next result with 4 inputs*)

NumericalMath`NSequenceLimit[{a,b,c,d}]
(*(-b^2+a*c)/(a-2*b+c)*)
(*yes, only references a, b, & c not the later, more converged, d*)

sequenceLimit[{a,b,c,d}]//Simplify (*my function*)
(*(-c^2+b*d)/(b-2*c+d)*)
(*note use of d, which explains why my approximations

  for Pi and E above are better*)

NumericalMath`NSequenceLimit[{a,b,c,d,e}]
(*(-b^2+a*c)/(a-2*b+c)*) 
(*yes, same result as with 4 terms, no d or e in output*)

sequenceLimit[{a,b,c,d,e}]//Simplify (*my function*)
(*(c^3+a*d^2+b^2*e-c*(2*b*d+a*e))/
  (b^2+3*c^2-2*c*d+d^2-2*b*(c+d-e)-c*e-a*(c-2*d+e))*)
(*note use of d and e*)

I also document some error messages that seem spurious.

NumericalMath`NSequenceLimit[{a}]
(*NumericalMath`NSequenceLimit::seqw Sequence of length 1
  is too short for use with Degree -> 1*)
(*NumericalMath`NSequenceLimit::bdmtd WynnEpsilon is not
  a valid specification of a sequence limit extrapolation
  algorithm.*)
(*NumericalMath`NSequenceLimit[{a}]*)


NumericalMath`NSequenceLimit[{a},Degree->1]; (*Nulled output*)
(*NumericalMath`NSequenceLimit:optx Unknown option ° in 
  NumericalMath`NSequenceLimit[{a}]*)

NumericalMath`NSequenceLimit[{a},"Degree"->1]; (*Note quotes*)
(*NumericalMath`NSequenceLimit::optx Unknown option "Degree"
  in NumericalMath`NSequenceLimit[{a}]*)

sequenceLimit[{a}] (*sequence limit is my function*)
(*a*)


Quiet@NumericalMath`NSequenceLimit[{a,b}] (*Note Quiet & no evaluation*)
(*NumericalMath`NSequenceLimit[{a,b}]*)

sequenceLimit[{a,b}] (*my function*)
(*b*)

Everything above is from the current version of Mathematica Online.

$Version
(*11.3.0 for Linux x86 (64-bit) (March 7, 2018)*)

Here is the definition of my function

wynnE[-2,_,_]:=0
wynnE[-1,n_,s_]:=s@n
wynnE[rkp1_,n_,s_]:=(*wynnE[rkp1,n,s]=*)(*optional caching*)
 wynnE[rkp1-2,n+1,s]+1/(wynnE[rkp1-1,n+1,s]-wynnE[rkp1-1,n,s]) 

sequenceLimit[list_?VectorQ]:=
 With[{len=Length@list},
  With[{off=Boole@EvenQ@len},

   wynnE[len-2-off,off,list[[#+1]]&]
  ]/;len>0
 ]