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

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

Earlier answer [answer](https://mathematica.stackexchange.com/q/280401/217) 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](https://github.com/HenrikSchumacher/DPR1Eigensystem) for this.

I'm looking to do this in $O(d)$ time by specializing to the case of convergent matrices. There's a suggestion that this is doable in Section 3 of this [paper](https://arxiv.org/abs/2206.11124) so I'm looking for hints on the best way to implement this (or a cheap alternative) in Mathematica.

------

Example below uses generating function approach to evaluate powers of convergent diagonal matrices. It uses  approximation below, then calls `LaplaceTransform` on integral obtained after `IntegrateChangeVariables`

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

It 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]
```

```
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]];
```