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

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Tue, Oct 01, 2024 at 01:55:57AM +0200, hermann@stamm-wilbrandt.de wrote:
> 
> 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?

This is not well defined for non-integral exponents.
In your case there are 4 possible square roots.

You can use issquare/ispower
? issquare(Mod(41,71641520761751435455133616475667090434063332228247871795429),&z)
%3 = 1
? z
%4 = Mod(19404231649676781033243368819091487975254979923473930822898,71641520761751435455133616475667090434063332228247871795429)
? z^2
%5 = Mod(41,71641520761751435455133616475667090434063332228247871795429)

Internally, we have a function Zn_quad_roots that compute all the solution of x^2+b*x+c mod N 
for composite N.
Maybe we could add it to GP if we find a GP interface to it.

Cheers,
Bill