What is the Mathematica "PowerMod[a, 1/2, m]" equivalent in PARI/GP?

[hermann@stamm-wilbrandt.de]

https://reference.wolfram.com/language/ref/PowerMod.html

I found this mailing list post from 2000:
https://pari.math.u-bordeaux.fr/archives/pari-users-0004/msg00001.html

But its function does not work for non-prime modulus:

hermann@7950x:~$ gp -q
? powermod(a,k,n)=lift(Mod(a,n)^k)
(a,k,n)->lift(Mod(a,n)^k)
? 
powermod(41,1/2,71641520761751435455133616475667090434063332228247871795429)
  ***   at top-level: powermod(41,1/2,716415207617514354551336164756
  ***                 ^----------------------------------------------
  ***   in function powermod: lift(Mod(a,n)^k)
  ***                                      ^---
  *** _^_: not a prime number in gpow: 
71641520761751435455133616475667090434063332228247871795429.
  ***   Break loop: type 'break' to go back to GP prompt
break>


PowerMod[] can deal with non-prime modulus — how can this be computed 
with PARI/GP?

pi@raspberrypi5:~ $ time wolframscript -code 
"PowerMod[41,1/2,71641520761751435455133616475667090434063332228247871795429]"
15567422879070002639383923810785206745982804843948310703484

real	0m59.452s
user	0m0.341s
sys	0m0.071s
pi@raspberrypi5:~ $


Regards,

Hermann.