Efficient powers of DPR1 matrices

Question

I'm looking to compute the following quantity for $A$ where $A$ is a convergent positive definite $d\times d$ diagonal + rank1 matrix (DPR1).

$$f(s)=\operatorname{Tr}(A^s)$$

Earlier answer answer by Henrik Schumacher shows how to do this in Mathematica in $O(d^2)$ time for an arbitrary positive definite DPR1 matrix, and provides a nice package for this.

I'm looking for the most elegant way to do this in Mathematica in $O(d)$ time by specializing to the case of convergent matrices. The math for this technique is described here

---

Example below uses LaplaceTransform to find $\operatorname{Tr}(A^s)$ when $A$ is diagonal. It relies on following approximation.

$$\sum_i (1-h_i)^s\approx \int_i \exp(-sh_i)$$

d = 5;
h[i_] = i^-1;
hvals = Table[h[i], {i, 1, d}];
diag[s_] := Tr[MatrixPower[DiagonalMatrix[1 - hvals], s]]/d;
diagInt = Inactive[Integrate][Exp[-s h[i]]/d, {i, 1, d}];

(* changes integration limits (a,b) -> (0,\[Infinity]) *)
expandIntegrationLimits[int_] := 
  Module[{body, limitSpec, var, limits, newLimitSpec, interval},
   (* returns expression that is zero when var is outside of (min,
   var) *)
   interval[var_, {min_, max_}] = 
    HeavisideTheta[var - min] - HeavisideTheta[var - max];
   
   Assert[Head[int] == Inactive[Integrate]];
   {body, limitSpec} = List @@ int;
   {var, limits} = {First[limitSpec], Rest[limitSpec]};
   newLimitSpec = {var, 0, \[Infinity]};
   Inactive[Integrate][body*interval[var, limits], newLimitSpec]
   ];


(* turn into Laplace transform form *)
int = IntegrateChangeVariables[diagInt, hi, hi == h[i]];
int = expandIntegrationLimits[int];
diagApprox[s_] = -LaplaceTransform[First@int/Exp[-hi s], hi, s];

s0 = 5;
Print["Exact value diag: ", diag[s0] // N]
Print["Integral approximation: ", 
 Activate[diagInt /. {Integrate -> NIntegrate, s -> s0}]]

Print["symbolic Laplace: ", N[diagApprox[s0]]];

(* normalize h to ensure DPR1 matrix is convergent *)
hvals = hvals/Total[hvals];
dpr[s_] := 
  Tr[MatrixPower[
     DiagonalMatrix[1 - hvals] + {hvals}\[Transpose] . {hvals}, s]]/
   d;
s0 = 5;
Print["Exact value dpr1: ", dpr[s0] // N];
Print["Approx value dpr1: ", dprApprox[s0]];

should print

Exact value diag: 0.145584
Integral approximation: 0.148296
symbolic Laplace: 0.148296
Exact value dpr1: 0.528644
Approx value dpr1: dprApprox[5]