Criterion for fundamental solutions to Pell's equation of the form $X^2-DY^2=C$

Question

We are trying to codify in terms of modern algorithm the works of the ancient Indian mathematician Udayadivakara (CE 1073). In his work Sundari, he quotes one Acarya Jayadeva who has given methods to solve Pell's equations. In these methods, one can find the the cyclic Chakravala method to deal with $X^2-DY^2=1$ wrongly attributed to Bhaskara. He also gives the method to solve $X^2-DY^2=C$ for any integer $C$.

1. His algorithm starts off by finding the nearest square integer $>D$ named $P^2$. Then $a=P^2-D$.
2. Now some $b$ is chosen in such a way that $Db^2+Ca$ is some perfect square $Q^2$.
3. Then the $X$ and $Y$ solutions can be found by using $Y=\frac{Q\pm P b}{a}$ and $X=PY \mp b$.
4. This procedure can continue indefinitely to find all the solutions.
5. Coming to the question of the fundamental solution i.e. the solution with which Bhavana has to be performed repeatedly to get other solutions (related to the modern automorphism group of the quadratic form), Prof. K.S. Shukla who first translated the work from Sanskrit to English, in his example says that it should be chosen "appropriately".
6. Our primary question is then what is the criterion to derive this fundamental solution? Is there a way to derive such a criterion?
7. The whole procedure seems to resemble Conway's topograph method which has been posted here several times Does the Pell-like equation $X^2-dY^2=k$ have a simple recursion like $X^2-dY^2=1$? It is quite fascinating to think that some wonderful mind came up with this algorithm about 1000 years ago and the optimality of it is equally amazing.

P.S: If anyone so wishes, we would be happy to provide a version of the original paper written by Prof. Shukla in 1950 dealing with this!

Answer 1

This not answer, only illustration how work pari/gp.

gp-code:

gpell(D,C)=
{
print("\nRoot solutions of Pell equation x^2-",D,"*y^2=",C,"\n");

 Q= iferr(bnfinit('X^2-D), E, 0);

 if(Q,

  U= iferr(Q.fu, E, 0);

  if(U, for(j=1, #U, u= U[j]; print("Q.fu: ",u,"\n");

   N= iferr(bnfisintnorm(Q, C), E, 0); print("bnfisnorm: ",N,"\n");

   if(N, for(k=1, #N, n= N[k];

    for(l=0, 48, \\print("\n",l);

     nu= lift(n*u^l); \\print(nu);

     X= abs(polcoeff(nu, 0)); Y2= (X^2-C)/D;  

     if(X==floor(X)&&Y2==floor(Y2),

      Y= sqrtint(Y2);

      print("X= ",X,"    Y= ",Y,"    l= ",l); break()

     )
    )
   ))
  ))
 )
};

Run code:

? \r gpell.gp
? gpell(97,96)

Root solutions of Pell equation x^2-97*y^2=96

Q.fu: Mod(569*X - 5604, X^2 - 97)

bnfisnorm: [25/2*X + 247/2, -5/2*X + 53/2, 47*X - 463, 6173*X - 60797, 111802*X - 110112
2, 2*X + 22, -2*X + 22, 278*X - 2738, -5035*X + 49589, -661225*X + 6512311, -5/2*X - 53/
2, 25/2*X - 247/2]

X= 463    Y= 47    l= 0
X= 60797    Y= 6173    l= 0
X= 1101122    Y= 111802    l= 0
X= 22    Y= 2    l= 0
X= 22    Y= 2    l= 0
X= 2738    Y= 278    l= 0
X= 49589    Y= 5035    l= 0
X= 6512311    Y= 661225    l= 0
?
?
?
?
? d=1377;gpell(d*(d-8),d*8)

Root solutions of Pell equation x^2-1885113*y^2=11016

Q.fu: Mod(1/333*X - 4, X^2 - 1885113)

bnfisnorm: [-11/74*X + 459/2, -X + 1377, -3/37*X + 153, -45/74*X + 1683/2]

X= 1377    Y= 1    l= 0
X= 74319362381115913874527035825931593    Y= 54129408430647687070089140606679    l= 36
?

Very attention to the parameter "l" in code.

And for contrast work Wolfram.

Answer 2

In the case of $X^2-DY^2=C$, this applies.

1. $C$ is the product of short and special primes, and some $Z^2$ relative to base $D$.
2. $Z=\gcd(X,Y)$ can be any number, but $Z^2\mid C$.
3. Where $D = 1 \pmod{4}$, then $X$ and $Y$ can be integer-halfs (ie both Z+½).

Let's see what this means.

In base $D$, the period of 1/p will divide $(p-1)/2$ an even (short) or odd (long) number of times. For example, in decimal, short primes are of the form $40n\pm 3^x$, so 3, 13, 31, 37 are 'short primes'. This means that the equation has a solution in terms of $C=111$ or $C=1209=3 * 13 * 31 $, for example, since these are the product of short primes relative to 10.

The long primes never occur as solutions unless they divide both $X$ and $Y$, in which case $p^2 \mid C$.

The special primes, are generally taken as the set of 2 and the divisors of $D$. These will necessarily accompany certain primes. This give for example, that $3$ will not occur alone, but requires a $2$ or $5$ to be with it. So there is no solution to $X^2-10Y^2=3$, but there is for $X^2-10Y^2=6$.

General Theory

Algebraically, these equate to values of the form $a \pm b\sqrt{D}$. This is an integer-system, in that every number of this form divides some integer $C=a^2-D.b^2$. Many of these integer systems are not completely factorable in themselves, but rely heavily on 'multiplication-half' systems.

For example, the decimal case $a\pm b\sqrt{10}$ has a half-multiples of the form $a\sqrt 2 \pm b \sqrt 5$. $a$ and $b$ can be integer-halfs when the difference of the numbers under the square roots is divisible by four, ie $4\mid (5-2)$ here.

Note also that there can be a unit in the integer-half system. $D=10$ does not have one, but there is certainly one in $D=42$ as $\sqrt 7 + \sqrt 6$. This system also has multiplication-halfs at $a\sqrt 2 + b\sqrt{21}$, and $a\sqrt 3+b \sqrt{14}$.

Note that any given prime occurs in a particular multiplication-half, and requires its special prime(s) to find a value of $C$. So in decimal, the prime $41$ can occur as a direct solution, whereas $37$ requires an additional special half-factor (2 or 5) to be a solution. An additional multiplication-half will do, so $13*37=481=10.7^2-3^2$ works.

The case of $D=130$ is interesting. It has four half-systems (1,130), (2,65), (5,26) and (10,13), which work on the XOR principle (ie a prime in (2,65) times one in (10,13), requires a prime from (5,12) to make the square (2×10×5 is square). But the unit for this system has no half-form.

I believe that if you take the primitive roots of the primes, and make a table of even and odd indexes (ie $g^i = n\pmod p$), then these rules apply.

1. If the sum of indexes for $D \pmod p$ is even, then $p$ can occur as a single-power divisor of $C$.

2. The actual parity of the individual divisors of $D$ determine which 'half' it falls in. For example, relative to $D=30$, we have the prime 13, as long in binary, short in -3, and long in 5. This gives 101, so the period of 13 in base 30 has a length dividing 6. But the unit for $x+y \sqrt{30}$ has a half falling in $\sqrt 6 + \sqrt 5$. So to make these all even, it's not in the set 000,110,001,111, (which require no helper prime), but in 010, 100, 011, 101, which needs either a $2$ or $3$. So we find for example, that 26 and 39 are possible candidates for $C$, but not 13 itself.