Autor Beitrag
Mathematiker
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Quizmaster
Beiträge: 2615
Erhaltene Danke: 1420

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: 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.
ausblenden volle Höhe 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..8of 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]=0and (n>0do dec(n);
      //Grad m von b bestimmen
      m:=8;
      while (b[m]=0and (m>0do dec(m);

      //Division, wenn Grad n mindestens gleich m ist
      if (m>0then
      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
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Fr 01.06.12 12:47 
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:
ausblenden 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>0then 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.

Für diesen Beitrag haben gedankt: Mathematiker
Mathematiker Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Quizmaster
Beiträge: 2615
Erhaltene Danke: 1420

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: Fr 01.06.12 15:15 
Hallo Gammatester,
Dein Hinweis ist korrekt: Ich habe vergessen zu erwähnen, dass die Koeffizienten ganze Zahlen sind.
user profile iconGammatester hat folgendes geschrieben Zum zitierten Posting springen:
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.
user profile iconGammatester hat folgendes geschrieben Zum zitierten Posting springen:
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 Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Quizmaster
Beiträge: 2615
Erhaltene Danke: 1420

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: 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:
ausblenden volle Höhe 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>0do dec(n);
      //Grad m von b bestimmen
      m:=maxgrad;
      while (FGIntCompareAbs(b[m],null)=eq) and (m>0do 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>0then
      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


Zuletzt bearbeitet von Mathematiker am So 17.06.12 23:22, insgesamt 1-mal bearbeitet
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1586
Erhaltene Danke: 231


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: 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.


Zuletzt bearbeitet von Delphi-Laie am So 17.06.12 23:37, insgesamt 1-mal bearbeitet
Mathematiker Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Quizmaster
Beiträge: 2615
Erhaltene Danke: 1420

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: So 17.06.12 23:18 
Hallo Delphi-Laie,
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
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.
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
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.

Für diesen Beitrag haben gedankt: Delphi-Laie
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1586
Erhaltene Danke: 231


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: 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.

Für diesen Beitrag haben gedankt: Mathematiker
Mathematiker Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Quizmaster
Beiträge: 2615
Erhaltene Danke: 1420

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: Mo 25.06.12 20:21 
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
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

Für diesen Beitrag haben gedankt: Delphi-Laie