@vv For a generalized eigenvalue problem (which is the case considered in this Question), symmetry doesn't guarantee real eigenvalues. For example, consider this case, where the matrices are real symmetric and the computation is exact:

LinearAlgebra:-Eigenvalues( <1,0; 0,-1>, <0,1; 1,0>);

@tomleslie It looks like the MaplePrimes worksheet rendering ignores displayprecision. Actually, it's not even clear whether its value is stored in the worksheet..

@vv Why do you use shape= symmetric? Is the original matrix symmetric?

@Carl Love If the intent is to use the above procedure iteratively, using a list of vectors as the second argument, the overall efficiency can be substantially improved for that case.

Continuing Preben's line of thought, for any matrix that's small enough that I'd enter it manually, I do this (in 1D input (aka Maple Notation)):

M:= < 
   1,   1         ;
   1,   2/3.4
>;

For me, this spacing and visual result lets me totally see it as a mathematical matrix. In addition, the prettyprinted result of the above (explicitly typed without palettes) command is identical regardless of whether you use 1D or 2D Input.

I've taken my suggestion above, and I've written a procedure Winding to more robustly compute the winding numbers. It checks for five error conditions (the details of which are clear in the following code, I hope). This is then incorporated into a slightly modified procedure ContourInt, whose callling protocol remains the same, with a new (optional) keyword argument time which specifies the maximum number of seconds to allow for each numeric integral. All computations for the winding number are done in hardware-float arithmetic to the extent possible. Of course, the final result of ContourInt is still exact symbolic.

Winding:= proc(
   C::algebraic,
   z0::complexcons, 
   trange::(name= range(realcons)), 
   {time::{positive, -1}:= 5} #see ?timelimit
)
description
   "Somewhat robust, yet speedy, exact computation of the winding number of a "
   "curve around a point in the complex plane via low-precision numeric integration"
;
option
   `Author: Carl J. Love <carl.j.love@gmail.com> 22 Feb 2019`,
   hfloat
;
local t:= lhs(trange), w, W;
   Digits:= trunc(evalhf(Digits)); #"environment" variable
   UseHardwareFloats:= 'deduced';  #"environment" variable
   try
      if evalf(abs(`-`(eval~(C, t=~ [op(rhs(trange))])[]))) > 1e-7 then
         error "Ends of contour do not meet"
      fi
   catch:
      error "Cannot numerically evaluate contour endpoints"
   end try;

   try
      w:= timelimit(
         `if`(time=infinity, -1, time), #see ?timelimit
         int(
            diff(C,t)/(C-z0), trange, 
            numeric, 'method'= '_d01akc', 'epsilon'= 1e-3, 'digits'= 4, _rest
         )
      );
      w:= evalhf(w/2/Pi/I)
   catch "time expired":
      error "Couldn't evaluate winding integral in under %1 seconds", time
   end try;

   W:= round(w);
   if W::integer implies abs(w-W) > 1e-2 then
      error "Winding number is not an integer"
   fi;
   W
end proc
:
      
ContourInt:= proc(
   f::algebraic, 
   C::(name= algebraic), 
   trange::(name= range(realcons)), 
   {time::{positive, -1}:= 10} #see ?timelimit
)
local r:= rhs(C), a;
   add(
      residue(f,a)*Winding(r, rhs(a), trange, :-time= time, _rest),
      a= remove(a-> type(rhs(a), infinity), op~({singular(f,lhs(C))}))
   )*(2*Pi*I)
end proc
:

Although you seem satisfied by Acer's Answer, it's not clear to me (and I doubt that it's clear to Acer) whether you want the result to be a new table or simply a list combining the entries of the original tables. And if you do want a new table, then how should duplicate indices be handled?

@C1Ron Yes, {z=p1} is a set, and the output of [singular(f(z), z)] is a list of sets. The command op (short for operand(s)) applied to a set or list will return its underlying sequence of elements. So, op({z=p1}) returns z = p1. Any command and almost all operators can have ~ appended to make them act elementwise, i.e., so that their action will apply to all elements of a container (the relevant containers being sets, lists, tables, and rtables (Vectors, Arrays, Matrices)). So, op~([singular(f(z), z)]) returns [z = p1, z = p2, ....]. Note that the outer container, in this case a list, remains the same.

If you change any add or mul command to seq, the sequence of operands will be returned without being added or multiplied.

@vv Ah, I understand now. Yes, your way is better than mine. Perhaps it'd be worthwhile to issue an error if the integral is significantly different from an integer.

@Carl Love Okay, here's a general procedure for symbolic contour integration (for simple closed contours, counterclockwise, blah, blah). For each point a returned by singular (other than infinity), I compute the numeric contour integral of 1/(z-a) around the contour. If the absolute value of this is significantly^[*1] larger than 0, then a is used.

^[*1]"Significantly" defaults to 10.^(2-Digits), and can be changed by the keyword parameter eps.

ContourInt:= proc(
   f::algebraic, 
   C::(name= algebraic), 
   trange::name= range(realcons), 
   {eps::positive:= 10.^(2-Digits)}
)
local z:= lhs(C), t:= lhs(trange), r:= rhs(C);
   add(
      map2(
         residue,
         f,
         select(
            Z-> abs(evalf(Int(diff(r,t)/(r-eval(z,Z)), trange))) > eps, 
            remove(z-> type(rhs(z), infinity), op~({singular(f,z)}))
         )
      )
   )*(2*Pi*I)
end proc:

Usage for your two cases:

ContourInt(f(sqrt(2),z), z= 2*exp(I*t), t= 0..2*Pi);
ContourInt(f(2,z), z= cos(t)+2*I*sin(t), t= 0..2*Pi);
 

I don't know whether this is more or less complicated or robust than computing the winding number, as suggested by VV above.

@Joe Riel I don't think that Maple V has 2D Input.

@Josolumoh The normalization process that I described won't work for d < 0 because the process requires B >= 0 for all x in the "support" and B > 0 for at least some x. I put "support" in quotes because the technical definition of it is "those x such that B > 0". But if there are some x in the "support" such that B = 0, we can easily work around that, but we can't work around B < 0.

I think that you'll get better search results if you change "cloud" to "cluster".

Do you need all the data to be below the hypothenuse (this seems relatively easy), or just most of it below?

@Josolumoh Yes, I have confirmed that for d = 0 the sum over the support x= 0..infinity is exactly 1 for any nonnegative integer r and any b > 0. Do you wish to proceed with this, using d = 0? 

To proceed for d > 0, the function B would need to be "normalized" for each b by dividing by the sum over the support for that r, d, and b. I'm a bit uncomfortable doing this because I don't understand the significance of the probability function. But it would be a valid probability mass function after the normalization. Perhaps you were thinking of the Beta negative binomial distribution (Wikipedia), also known as the inverse Markov-Polya distribution or the generalized Waring distribution?

@emendes Sorry for putting a mysterious number, 49, into the code. That was my bad. It came from your 50, minus 1. In the code below, I've improved that by renaming your 50 as nspan.

The following code uses a single procedure NestList for hardware floats, software floats, and exact arithmetic:

NestList:= proc(f, x, n::nonnegint, nmlzr:= (), {digits::posint:= trunc(evalhf(Digits))})
uses LT= ListTools;
local 
   k, R:= rtable(0..n, [x]),
   Iterate:= proc(f, R::rtable, n::nonnegint, nmlzr:= (x-> x))
   local k;
      for k to n do R[k]:= nmlzr(f(R[k-1])) od;
      R
   end proc,
   Types:=                [complex(float), radalgnum,    anything],
   nmlzrs:= table(Types=~ [x-> x,          evala@Normal, simplify]),
   NMLZR:= `if`(nmlzr = (), nmlzrs[LT:-SelectFirst((t,x)-> x::t, Types, x)], nmlzr)
;
   Digits:= digits;
   if x::complex(float) and Digits <= trunc(evalhf(Digits)) then
      R:= hfarray(0..n, [x]); 
      try return evalhf(Iterate(f, R, n, x-> x)) catch: print("evalhf failed") end try 
   fi;
   Iterate(f, R, n, NMLZR)
end proc
:
DoPlot:= (aaa::rtable)->
   plot(
      Re~(<<[$-nstep-1..nspan-1]> | aaa^%T >), 
      style= pointline, symbol= solidcircle, symbolsize= 4, 
      thickness= 0, view= [default, 0..1], _rest
   )
:
(nstep, nspan):= (12, 50);
logistic:= x-> 4*x*(1-x):
fp:= [solve((logistic@@5)(x) = x, x)]:
cfp:= allvalues(fp[5]):
p:= nstep+nspan:
plots:-display(<
   DoPlot(NestList(logistic, evalhf(Re(cfp)), p), title= "Using evalhf") |
   DoPlot(
      NestList(logistic, evalf[30](Re(cfp)), p, digits= 30), 
      title= "Using 30 Digits"
   ) |
   DoPlot(NestList(logistic, cfp, p), title= "Using exact arithmetic")
>);