Re: Pari/GP "halfgcd()" implementation

[hermann@stamm-wilbrandt.de]

> 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).
I tried it, and it is really slow!

tail of "git diff 388342.halfgcd.gp", with your "myhalfgcd()" added at 
top:
...
 ##
-assert(V[2]^2 + M[2,1]^2 == p, "V[2],M[2,1]", " not sum of squares");
+[x,y] = [V[2], M[2,1]];
+assert(x^2 + y^2 == p, "x,y", " not sum of squares");
+
+[M,V]=myhalfgcd(sqrtm1,p);
+##
+[x,y] = [V[2], M[2,1]];
+assert(x^2 + y^2 == p, "x,y", " not sum of squares");
+
 print("done");


Output on i7-11850H, that is a slowdown of more than 13000x ...

$ python -c "print((13*60000+1518)/60)"
13025.3
$

$ gp-2.16 < 388342.halfgcd.gp
...
parisizemax = 2000003072, primelimit = 1048576, nbthreads = 8
foobar
  ***   last result: cpu time 60 ms, real time 60 ms.
  ***   last result: cpu time 13min, 1,518 ms, real time 13min, 4,658 
ms.
done
Goodbye!
$


Regards,

Hermann Stamm-Wilbrandt.