Autor |
Beitrag |
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
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.
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 fillchar(a,sizeof(a),0); fillchar(b,sizeof(b),0); 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 fillchar(q,sizeof(q),0); fillchar(r,sizeof(r),0); n:=8; while (a[n]=0) and (n>0) do dec(n); m:=8; while (b[m]=0) and (m>0) do dec(m);
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; r:=a; end;
ende:=true; for i:=0 to 8 do if r[i]<>0 then ende:=false;
a:=b; b:=r; until ende; 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
|
Verfasst: 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:
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.
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
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
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
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:
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; 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 numberToFGInt(0,null);
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); 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;
repeat for i:=0 to maxgrad do numberToFGInt(0,q[i]); for i:=0 to maxgrad do numberToFGInt(0,r[i]); n:=maxgrad; while (FGIntCompareAbs(a[n],null)=eq) and (n>0) do dec(n); m:=maxgrad; while (FGIntCompareAbs(b[m],null)=eq) and (m>0) do dec(m);
FGIntMod(a[n],b[m],rest);
if FGIntCompareAbs(rest,null)<>eq then begin FGIntGCD(a[n],b[m],ff); FGIntDiv(a[n],ff,f1); FGIntMul(b[m],f1,ff); FGIntDiv(ff,a[n],fa); FGIntDiv(ff,b[m],fb);
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;
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; for i:=0 to maxgrad do r[i]:=a[i]; end;
ende:=true; for i:=0 to maxgrad do if FGIntCompareAbs(r[i],null)<>eq then ende:=false;
a:=b; b:=r; until ende;
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;
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
Beiträge: 1600
Erhaltene Danke: 232
Delphi 2 - RAD-Studio 10.1 Berlin
|
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.
Zuletzt bearbeitet von Delphi-Laie am So 17.06.12 23:37, insgesamt 1-mal bearbeitet
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
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.
Für diesen Beitrag haben gedankt: Delphi-Laie
|
|
Delphi-Laie
Beiträge: 1600
Erhaltene Danke: 232
Delphi 2 - RAD-Studio 10.1 Berlin
|
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.
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
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
Für diesen Beitrag haben gedankt: Delphi-Laie
|
|
|