rを適当な原子根として、任意の原始根の集合は( ri(modp) : 1≤i≤p−1 , gcd(p−1,i)=1 )で与えられるため、
ans = sum(ri(modp) : 1≤i≤p−1) - sum( ri(modp) : 1≤i≤p−1 , gcd(p−1,i)>1 ) modP
を求めれば良いです。
前者は対数時間で、後者はP−1の素因数を列挙したのちに、包除原理を用いて十分高速に計算できます。
計算量はO(P+2(P−1の素因数の個数)∗log(P)) です。
【余談】テストケースは全てP−1,0,1のケースの場合ですが、これは意図的にそうしたのではなく、他のケースが見つからなかったためです。
少なくともPが1000000以下においてはそのようなケースは確認されませんでした。
これが任意の素数Pについて成り立つかは確認できていませんので、興味のある方はぜひ真偽を確かめてください。