Given that you only have two phases (0 and pi) in your input, if p is fixed you can easily precompute that and just multiply the sign with the result. Then you multiply with the p-norm divided by the vector size to 1/p (or absorb that factor into the phase).
Best regards
Thomas
Edited P.S.: Looking at the paper authors’ implementation, they use min, max (limits as p approaches -/+inf,  mean and 3rd-power mean, so it seems simple enough and no complex numbers involved (if you take power of 1/3 to be the inverse of power of 3).