The sum of unitary divisors of a positive integer $n$ can be computed as $$ \sigma^\star(n) = \prod_{i=1}^{\omega(n)}(p_i^{a_{\large i}}+1) $$ where $\prod_{i=1}^{\omega(n)} p_i^{a_{\large i}}$ is the prime factorisation of $n$; assuming that $\sigma^\star(1)=1$, being the empty product.
$\sum_{n=1}^x \sigma^\star(n)$ can either be computed exactly by factorising all numbers up to $x$, or it can be approximated as $\frac{\pi^2 x^2}{12\zeta(3)} + \mathcal{O}(x \thinspace log^{2/3}x)$.
For a prime $p \ll x$, how can we compute the exact modular sum $$ \sum_{\large n=1}^{\large x} \sigma^\star(n) \mod p $$ efficiently, doing fewer than $x$ prime factorisations?