Pellsche Gleichung

Beispiel für die Lösung mittels Kettenbruchentwicklung

Wir wollen die Gleichung

lösen, d. h. für .

All erstes benötigen wir die Kettenbruchentwicklung für . Diese berechnet sich durch rekursives Anwenden folgender Schritte:

• Berechnung von : Dies ist der Wert, den wir zur genauen Darstellung von durch einen hier abgebrochenen Kettenbruch benötigen würden.

Im ersten Schritt ist der Ausgangsterm , in den folgenden ist der Restterm .

Der Nenner der Form wird durch Multiplikation mit der konjugierten Wurzel ganzzahlig gemacht, anschließend wird, wenn möglich, gekürzt.

• Berechnung von : Dies ist der ganzzahlige Teil (Vorkommastellen) von . Durch Einsetzen von für kann man das Abrunden von unter Umgehung von irrationalen Zahlen durchführen.
• Berechnung von : Dies ist der Rest (Nachkommastellen) von und berechnet sich zu .

Erhalten wir hierbei ein schon mal vorgekommenes , können wir die Entwicklung hier abbrechen, die Kettenbruchentwicklung ist ab hier periodisch.

• Die Kettenbruchentwicklung entsteht aus den Werten für und lautet .
• Die Rechnung verläuft dann so:

Da ist, wiederholen sich ab hier die Werte periodisch und wir können die Entwicklung abbrechen.

Die Kettenbruchentwicklung für ist daher

und hat die Periode .

Die Näherungsbrüche der Kettenbruchentwicklung erhalten wir durch Einsetzen der Kettenbruchentwicklung in die Bildungsvorschrift. Die Näherungsbrüche bei Abbruch nach Stellen lauten:

Nach Umwandlung in gewöhnliche Brüche interessieren jetzt die Werte nach Abbruch nach einem ganzzahligen Vielfachen von :

und findet an den Stellen , , und

						die Triviallösung von
						die erste Lösung der negativen Pellschen Gleichung
						die erste (nicht triviale) Lösung von
						die zweite Lösung der negativen Pellschen Gleichung von

Weiter stellt man fest, dass für jedes Element der abgebrochenen Kettenbruchentwicklung der Länge eine Lösung einer Pellschen Gleichung mit rechter Seite ist.

Die Näherungsbrüche dazwischen stellen einige (aber nicht alle!) Lösungen der verallgemeinerten Pellschen Gleichung mit rechter Seite und dar.

#include <stdio.h>
#include <stdarg.h>
#include <cmath>
#include <vector>
#include <cstdint>
#include <stdexcept>
#include <utility>

using namespace std;
using i64 = int64_t;
using u64 = uint64_t;

static void log (const char* format, ...)
{
    if (1)
    {
        va_list  args;
        va_start(args, format);
        vprintf (format, args);
        va_end  (args);
    }
}

vector<u64> const Fahrradkettenbruch (u64 const n)
{
    log ("Kettenbruch-Entwicklung von sqrt %llu\n\n", n);

    u64 const   b0 = (u64)sqrt(n);
    u64         b  = b0;
    u64         m  = 0;
    u64         d  = 1;
    vector<u64> cfs;

    cfs.push_back (b);

    log ("  b%-2d =%3lld = (sqrt %lld +%4lld) / %lld\n", 0, b, n, m, d);

    for (int i = 1; ; i++)
    {
        m = d * b - m;
        d = (n - m * m) / d;
        if (d == 0)
            break;           // Abbruch, falls d = 0 (n ist ein Quadrat)
        b = (b0 + m) / d;
        cfs.push_back (b);

        log ("  b%-2d =%3lld = (sqrt %lld +%4lld) / %lld\n", i, b, n, m, d);

        if (b == 2*b0) break;  // Periodenende erkennen
    }
    return cfs;
}

pair <u64, u64> SolvePell (u64 const n, i64 const d)  // Funktion zur Berechnung der Lösung der Pellschen Gleichung
{
    auto cfs = Fahrradkettenbruch (n);
    u64  p0  = 1     , q0 = 0;
    u64  p1  = cfs[0], q1 = 1;

    log ("\nBerechnung der Naeherungsbrueche\n");
    log (                "  p%-2d = %lld\n  q%-2d = %lld\n\n",            0, p0, 0, q0);
    log ("  b%-2d = %lld\n  p%-2d = %lld\n  q%-2d = %lld\n\n", 0, cfs[0], 1, p1, 1, q1);

    for (int i = 1; ; i++)
    {
        u64 const b = cfs [(i-1) % (cfs.size()-1) + 1];
        u64 const p = b*p1 + p0;
        u64 const q = b*q1 + q0;
        if (1.0*b*p1 + p0 > 0xFFFFFFFFFFFFFFFF || 1.0*b*q1 + q0 > 0xFFFFFFFFFFFFFFFF)
        {
            log ("\n\nOverflow u64\n");
            return { 0, 0 };
        }

        log ("  b%-2d = %lld\n  p%-2d = %lld*%lld + %lld = %lld\n  q%-2d = %lld*%lld + %lld = %lld\n", i, b, i, b, p1, p0, p, i, b, q1, q0, q);

        p0 = p1;
        q0 = q1;
        p1 = p;
        q1 = q;

        i64 const dd = p1 * p1 - n * q1 * q1;
        log ("    -->%+4lld = %llu^2 - %llu * %llu^2\n", dd, p1, n, q1);

        if (dd == d) return { p1, q1 };
    }
}

void SolveAndPrint (u64 const n, i64 const d)
{
    log ("\n=== Berechnung fuer x^2 - %llu y^2 = %lld ===\n", n, d);
    auto const [x, y] = SolvePell (n, d);
    if (d == x*x - n * y*y)
        printf ("\nDie kleinste Loesung der Pellschen Gleichung x^2 - %lld y^2 = %+lld ist\n"
                "  x = %llu\n  y = %llu\n", n, d, x, y);
}

int main (int argc, char** argv)
{
    SolveAndPrint (argc < 2 ? 13 : ::atoll (argv[1]),   // n, default ist 13
                   argc < 3 ? +1 : ::atoll (argv[2]));  // d, default ist +1
    return 0;
}