ggT zweier Polynome

Mathematiker · Verfasst: Fr 01.06.12 11:34

Hallo Delphi-Gemeinde,
ich habe wieder einmal ein Problem, wahrscheinlich eher ein mathematisches, als ein programmtechnisches.
Ich beabsichtige den größten gemeinsamen Teiler zweier Polynome a(x) und b(x) zu berechnen. Dazu nutze ich den Euklidischen Algorithmus.

markieren

Delphi-Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:

procedure ggt;
var  a,b,q,r:array[0..8] of real;
     ende:boolean;
     n,m,i:integer;
begin
    //Felder für Koeffizienten der Polynome
    fillchar(a,sizeof(a),0);
    fillchar(b,sizeof(b),0);
    //Testwerte
    a[3]:=4;    a[2]:=5;   a[1]:=2;
    b[5]:=-16;  b[4]:=4;   b[3]:=26;
    b[2]:=33;   b[1]:=22;  b[0]:=8;

    repeat
      //Quotient auf Null setzen
      fillchar(q,sizeof(q),0);
      //Rest auf Null setzen
      fillchar(r,sizeof(r),0);
      //Grad n von a bestimmen
      n:=8;
      while (a[n]=0) and (n>0) do dec(n);
      //Grad m von b bestimmen
      m:=8;
      while (b[m]=0) and (m>0) do dec(m);

      //Division, wenn Grad n mindestens gleich m ist
      if (m>0) then
      begin
        while (n>=m) do
        begin
          q[n-m]:=a[n]/b[m];
          a[n]:=0;
          for i:=1 to m do a[n-i]:=a[n-i]-q[n-m]*b[m-i];
          dec(n);
        end;
        //in q steht der Quotient, in a der Rest, der nach r kopiert wird
        r:=a;
      end;

      //Test auf Rest = 0
      ende:=true;
      for i:=0 to 8 do if r[i]<>0 then ende:=false;

      //Polynome tauschen
      a:=b;
      b:=r;
    until ende;
//in a stehen am Ende die Koeffizienten des ggT(a,b)
//im Beispiel a[2]=16, a[1]=20, a[0]=8, was allerdings falsch ist,
//richtig wäre a[2]=4, a[1]=5, a[0]=2
end;

Eigentlich müsste es funktionieren, tut es aber nicht. Im konkreten Beispiel (siehe Quelltext) ist das Ergebnis falsch.
Für andere Koeffizienten ergibt sich manchmal das richtige Ergebnis, jedoch leider nicht immer.
Hat jemand eine Idee, wo mein Fehler ist? :?!?:

Beste Grüße
Mathematiker

Gammatester · Beiträge: 328 Erhaltene Danke: 101

Ich sehe mindestens zwei Probleme:

Erstens ist es relativ unklar, woher die Koeffizienten kommen sollen: R, Z, Q, Z/pZ oder was. Selbst wenn man Dein Beispiel in Pari/GP eingibt, erhält man:

markieren

Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:

(12:26) gp > a = 4*x^3 + 5*x^2 +2*x
%1 = 4*x^3 + 5*x^2 + 2*x
(12:26) gp > b = -16*x^5 +4*x^4+26*x^3+33*x^2+22*x+8
%2 = -16*x^5 + 4*x^4 + 26*x^3 + 33*x^2 + 22*x + 8
(12:26) gp > gcd(a,b)
%3 = 4*x^2 + 5*x + 2
(12:26) gp > b = -16.0*x^5 +4*x^4+26*x^3+33*x^2+22*x+8
%4 = -16.00000000000000000000000000*x^5 + 4*x^4 + 26*x^3 + 33*x^2 + 22*x + 8
(12:26) gp > gcd(a,b)
%5 = 16.00000000000000000000000000*x^2 + 20.00000000000000000000000000*x + 8

Wenn die Koeffizienten aus einem Körper sind, wird wohl sinnvollerweise der höchste Koeff. des ggt auf 1 normiert. Selbst für a,b aus Z[x] sollte man mM durch den ggt der Rest-Koeffzienten dividieren (in Deinem Fall 4, dann hättest Du auch das gewünschte Ergebnis); allerdings solltest Du dann auch Integerarithmetik benutzen.

Zweitens: Dein if (m>0) then verhindert, das zB für a=3, b=6 aus Z[x] überhaupt gerechnet wird, und das Ergebnis ggt(a,b)=6 ist doch wohl sehr zweifelhaft.

Mathematiker · Verfasst: Fr 01.06.12 15:15

Hallo Gammatester,
Dein Hinweis ist korrekt: Ich habe vergessen zu erwähnen, dass die Koeffizienten ganze Zahlen sind.

Gammatester hat folgendes geschrieben :

allerdings solltest Du dann auch Integerarithmetik benutzen.

genau hier liegt ein Problem. Während der Rechnung treten in dem Quotientenfeld q auch gebrochene Zahlen (im Beispiel 1/4) auf; mit dem Ergebnis, dass die Routine hängen bleibt.

Gammatester hat folgendes geschrieben :

Selbst für a,b aus Z[x] sollte man mM durch den ggt der Rest-Koeffzienten dividieren (in Deinem Fall 4, dann hättest Du auch das gewünschte Ergebnis);

Dank für den Tipp. Ich werde testen, ob dies mein Problem löst.
Beste Grüße
Mathematiker

Mathematiker · Verfasst: So 17.06.12 22:15

Hallo,
nach einiger Zeit habe ich das Problem des ggT zweier Polynome gelöst. Vielleicht interessiert es ja jemanden.
Die Schwierigkeit ist, dass während des Euklidischen Algorithmus sehr große Koeffizienten auftreten können. Aus diesem Grund habe ich die Unit FGInt von Walied Othman (siehe www.koders.com/delph...41A8033A3783AD6.aspx) verwendet.

Nachfolgend der Delphi-Text zu Berechnung des ggTs der Polynome a und b:

markieren

Delphi-Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:

uses fgint;
procedure polyggt;
const maxgrad = 8;         //größtmögliche Potenz im Moment auf 8 gesetzt  
var  a,b,q,r:array[0..maxgrad] of TFGInt;  
     ff,fa,fb,rest,f1,f2,null: TFGint;
     n,m,i:integer;
     ende,abbrucha,abbruchb:boolean;
begin
     //Nullpolynom
     numberToFGInt(0,null);

     //Hier ganzzahlige(!) Koeffizienten des ersten Polynoms a[i] und des zweiten b[i] festlegen
     // Absolutglied in a[0], lineares Glied in a[1], ... 
     // 

    //Test auf Nullpolynome
    abbrucha:=true;
    abbruchb:=true;
    for i:=0 to maxgrad do abbrucha:=abbrucha and (FGIntCompareAbs(a[i],null)=eq);
    for i:=0 to maxgrad do abbruchb:=abbruchb and (FGIntCompareAbs(b[i],null)=eq);
    //Abbruch wenn ein Nullpolynom vorliegt
    if abbrucha or abbruchb then
    begin
       for i:=0 to maxgrad do FGIntDestroy(a[i]);
       for i:=0 to maxgrad do FGIntDestroy(b[i]);
       FGIntDestroy(null);
       exit;
    end;

    //ggT-Bestimmung mit Euklidischem Algorithmus
    repeat
      //Quotient auf Null setzen
      for i:=0 to maxgrad do numberToFGInt(0,q[i]);
      //Rest auf Null setzen
      for i:=0 to maxgrad do numberToFGInt(0,r[i]);
      //Grad n von a bestimmen
      n:=maxgrad;
      while (FGIntCompareAbs(a[n],null)=eq) and (n>0) do dec(n);
      //Grad m von b bestimmen
      m:=maxgrad;
      while (FGIntCompareAbs(b[m],null)=eq) and (m>0) do dec(m);

      //Rest der Division der höchsten Koeffizienten von a und b
      FGIntMod(a[n],b[m],rest);

      //wenn nicht Null, dann Koeffizienten der höchsten Potenzen auf gleichen Wert erweitern
      if FGIntCompareAbs(rest,null)<>eq then
      begin
        //kgv suchen
        FGIntGCD(a[n],b[m],ff);
        FGIntDiv(a[n],ff,f1);
        FGIntMul(b[m],f1,ff); //kgV steht jetzt in ff
        //Erweiterungsfaktoren bestimmen 
        FGIntDiv(ff,a[n],fa);
        FGIntDiv(ff,b[m],fb);

        //Erweitern
        for i:=0 to n do begin
          FGINtmul(fa,a[i],f1); a[i]:=f1
        end;
        for i:=0 to m do begin
          FGINtmul(fb,b[i],f1); b[i]:=f1
        end;
      end;

      //Division, wenn Grad n mindestens gleich m ist
      if (m>0) then
      begin
        while (n>=m) do
        begin
          fgintdiv(a[n],b[m],q[n-m]);
          for i:=0 to m do begin
             fgintmul(q[n-m],b[m-i],f1);
             fgintsub(a[n-i],f1,f2);
             a[n-i]:=f2;
          end;
          dec(n);
        end;
        //in q steht der Quotient, in a der Rest
        for i:=0 to maxgrad do r[i]:=a[i];
      end;

      //Test auf Rest = 0
      ende:=true;
      for i:=0 to maxgrad do if FGIntCompareAbs(r[i],null)<>eq then ende:=false;

      //Polynome tauschen
      a:=b;
      b:=r;
    until ende;

    //Koeffizienten um ggT kürzen
    f1:=a[0];
    for i:=1 to maxgrad do begin
      FGIntGCD(a[i],f1,f2);
      f1:=f2;
    end;
    for i:=0 to maxgrad do begin
      FGIntdiv(a[i],f1,f2);
      a[i]:=f2;
    end;

    //zum Abschluss steht das ggT-Polynom im Feld a

    //Variable freigeben
    for i:=0 to maxgrad do FGIntDestroy(a[i]);
    for i:=0 to maxgrad do FGIntDestroy(b[i]);
    for i:=0 to maxgrad do FGIntDestroy(q[i]);
    for i:=0 to maxgrad do FGIntDestroy(r[i]);
    FGIntDestroy(rest);
    FGIntDestroy(ff);
    FGIntDestroy(fa);
    FGIntDestroy(fb);
    FGIntDestroy(f1);
    FGIntDestroy(f2);
    FGIntDestroy(null);
end;

Viel Spaß beim Testen.

Beste Grüße
Mathematiker

Delphi-Laie · Verfasst: So 17.06.12 23:13

Interessieren tut es mich durchaus, und auch meine Anerkennung für dieses Ergebnis.

Nur, was mich von Anfang an ein wenig irritierte: Ich nehme an, daß Dein Pseudonym Dein (nicht nur Delphi-)Programm ist. Insofern verwundert es mich, daß das für Dich eine Herausforderung sein soll.

Die Polynomdivision kenne ich auch, obwohl es schon 20 Jahre her ist, immer noch. Wenistens bei einer Variablen kann ich sie auch noch (mehr lernte ich als Nichtmathematiker ohnehin nie kennen).

Die Beschränkung auf den 8. Grad halte ich für ein wenig zu restriktiv. Ist bei heutigen Computerleistungen nicht ein Grad von "fast unendlich" möglich?

Edit: Die verlinkte Großzahlenunit ist sehr interessant, und sie ist sogar schon mit Delphi 4 compilierbar. Zur Zeit gefällt mir für solche Zwecke allerdings BenBes BigNum2-Unit am besten.

Mathematiker · Verfasst: So 17.06.12 23:18

Hallo Delphi-Laie,

Delphi-Laie hat folgendes geschrieben :

Ich nehme an, daß Dein Pseudonym Dein (nicht nur Delphi-)Programm ist. Insofern verwundert es mich, daß das für Dich eine Herausforderung sein soll.

Doch ist es, leider.

Delphi-Laie hat folgendes geschrieben :

Die Beschränkung auf den 8. Grad halte ich für ein wenig zu restriktiv.

Ist kein Problem. Es kann leicht abgeändert werden. Allerdings hätte ich auch eine Konstante einführen können, was weniger "Stress" beim Ändern macht. Sorry!
Beste Grüße
Mathematiker

Nachtrag: Ich habe die 8 auf maxgrad geändert.

Delphi-Laie · Verfasst: Mo 25.06.12 19:54

Hallo Mathematiker, meine Beschäftigung mit den Polynomen wird nunmehr - nach 20 Jahren - wieder konkreter. Zum Glück liegen mir alle meine damaligen (Turbo-)Pascal-Programme, wenn auch nur als Ausdrucke, noch vor.

Zunächst versuche ich, Deinen Prozedur zum Laufen zu bringen, doch es scheitert schon an diesem Aufruf:

numberToFGInt(0,null).

Die Unit "FGInt" deklariert diese Prozedur nicht. Ist das eine "Eigenbastelei" Deinerseits?

Überhaupt, diese Unit macht es (mir) ungleich schwieriger als BenBes BigNum-Units, sich (mich) in sie einzuarbeiten. Welches sind die "Schnittstellenprozeduren", also String <-> TFGInt sowie Integer (oder Int64) <-> TFGInt? (Nur so in den Raum geworfen.)

Vorausdank und Gruß

Delphi-Laie

Edit: Die letzte, eher rhetorische Frage hat sich wahrscheinlich erledigt, es scheinen die beiden Prozeduren

Base10StringToFGInt(Base10 : String; Var FGInt : TFGInt);
FGIntToBase10String(Const FGInt : TFGInt; Var Base10 : String);

zu sein.

Mathematiker · Verfasst: Mo 25.06.12 20:21

Delphi-Laie hat folgendes geschrieben :

Zunächst versuche ich, Deinen Prozedur zum Laufen zu bringen, doch es scheitert schon an diesem Aufruf:
numberToFGInt(0,null).
Die Unit "FGInt" deklariert diese Prozedur nicht.

Ich vermutete es schon fast, habe aber nicht nochmal nachgesehen. Sorry.
In der von mir verwendeten, "alten" FGInt gab es die Funktion noch. In der gegenwärtig im Netz verfügbaren Version hat der Autor ein paar Änderungen vorgenommen, u.a. alle longword gegen int64 getauscht. numberToFGInt(0,null) gibt es nicht mehr, kann aber durch Base10StringToFGInt ersetzt werden, wie Du selbst schon gefunden hast.
Beste Grüße
Mathematiker

_________________
Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme