Re: Pell's equations and beyond

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

* hermann@stamm-wilbrandt.de [2023-12-18 19:32]:
> I stumbled over
> https://math.stackexchange.com/a/3341210
> 
> on finding solution for Pell's equation x^2-D*y^2=1 for D=61.
> 
> Then I implemented my pell.gq (bottom) that did the job for any D.
> 
> Then I found this 2008 posting from Karim:
> https://pari.math.u-bordeaux.fr/archives/pari-users-0811/msg00001.html
> 
> 
> 1) How do these commands from Karim's posting reveal x and y?
> 
> ? quadunit(61)
> %15 = 17 + 5*w

For the equation x^2 - D*y^2 = ± 1, you want the fundamental unit in the order
of discriminant 4*D:

? quadunit(61*4)
%1 = 29718 + 3805*w
? norm(%)  \\ the fundamental unit has norm -1, not 1.
%2 = -1
? %1^2 \\ the fundamental unit of positive norm is its square
%3 = 1766319049 + 226153980*w

> ? K = bnfinit(x^2 - 61); K.fu
> %16 = [Mod(-5/2*x + 39/2, x^2 - 61)]

This has better complexity but is more complicated for your case,
because nf / bnf are only work (out of the box) in maximal orders.

? K = bnfinit(x^2 - 61); u = K.fu[1]
%1 = [Mod(-5/2*x + 39/2, x^2 - 61)]

This has a denominator (the maximal order is Z[(1 + sqrt(61)) / 2], not
Z[sqrt(61)]) so it doesn't belong to the order you want. You can
determine naively the power needed for u^i to belong to Z[sqrt(61)]

? for (i=2,oo, if (denominator(u^i, 1) == 1, return(i)))
%2 = 3

or non-naively BUT the following has the wrong complexity in D
(in D^{1/2} instead of D^ε); it has better complexity in the index
aspect though (which is 2 in this case, so this latter point is
irrelevant)

? quadunitindex(61, 2)
%3 = 3

More advanced note: I could modify quadunitindex() trivially for the
important special case where the fundamental unit is already known
(allowing quadunitindex(D, f, u)). So *that* part of the computation
would run in complexity (Df)^ε. But in bad cases, u itself can be as
large as exp(sqrt(D)), so just writing it down suffers from sqrt(D)
complexity ...

In any case

? u^3  \\ unsurprisingly we recover the x,y found previously, up to sign
%4 = Mod(-3805*x + 29718, x^2 - 61)

? norm(u) \\again, u (hence u^3) has negative norm, so take its square
%5 = -1

? u^6
%6 = Mod(-226153980*x + 1766319049, x^2 - 61)

> pell.gp determines x and y given D and verifies result being 1:
> 
> pi@raspberrypi5:~ $ D=61 gp -q < pell.gp
> period=11
> CF=[7, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14, 1, 4, 3, 1, 2, 2, 1, 3, 5]
> [x,y]=[1766319049, 226153980]
> x^2-D*y^2=1
> pi@raspberrypi5:~ $
> 
> 
> 2) Is there a simpler way to compute continued fraction period than calling
> "period()" from contfrac.gp in examples directory?
> 
> $ cat pell.gp
> readvec("pari-2.15.4/examples/contfrac.gp");
> D=eval(getenv("D"));
> p=period(D);
> print("period=",p);
> CF=contfrac(sqrt(D),,2*p);
> print("CF=",CF);
> [x,y]=contfracpnqn(CF,2*p-1)[,2*p-1];
> print("[x,y]=",[x,y]);
> print("x^2-D*y^2=",x^2-D*y^2);
> $

It doesn't really matter: by computing the continued fraction expansion
you get the "bad" complexity D^{1/2}. You might as well use quadunit ...

Cheers,

    K.B.
-- 
Pr. Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/