Does anyone know any efficient algorithm to compute the $p$-adic expansion of $\zeta_p(n)$ with $n\in \mathbb{Z}$ and $n \geq 3$, where $\zeta_p$ is the $p$-adic zeta function.
I also hope to implement the algorithm in Mathematica, does anyone know any Mathematica package for $p$-adic analysis?