Re: Pari/GP "halfgcd()" implementation

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

On Thu, Jun 08, 2023 at 10:56:19PM +0200, hermann@stamm-wilbrandt.de wrote:
> I really like Pari/GP "halfgcd()" function:
> https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.15.1/users.pdf#page=204
> 
> It allows very fast computation even on 388342-digit prime:
> Allows to determine  sqrt(-1) (mod p)  given sum of squares  p = x^2 + y^2
> (in 123ms on i7-11850H):
> https://forums.raspberrypi.com/viewtopic.php?p=2112166#p2112166
> 
> I wanted to see whether it is possible to make "halfgcd()" available for
> Python.
> Just did analyze "halfgcdii()" implementation in master branch
> "pari/Olinux-aarch64/mpker.c".
> The huge C callgraph made me step away from that plan, I will simply use
> Pari/GP instead.

Well, if you want a slow halfgcd, it is not harder than extended GCD:

myhalfgcd(a,b)=
{
  my(M=matid(2));
  my(n = max(a,b));
  while (b^2>=n,
    my([q,r]=divrem(a,b));
    [a,b] = [b,r];
    M = [0,1;1,-q]*M;
  );
  [M,[a,b]~];
}

(the matrix product can be inlined.)

halfgcdii implements both the Lehmer-Jebelean acceleration and the
subquadratic algorithm (following Thull-Yap).

Cheers,
Bill