4programmers.net 4p Forum Delphi i Pascal Algorytm Karacuby do mnożenia wielomianów

Algorytm Karacuby do mnożenia wielomianów

Algorytm Karacuby do mnożenia wielomianów
nowy121105
N1
  • Rejestracja: dni
  • Ostatnio: dni
0

Na kanale Przemysława Koprowskiego widziałem następujący pseudokod algorytmu Karacuby

Kopiuj
Algorytm Karacuby

Wejście : wielomiany
      f = f[0] + f[1]x+f[2]x^2+...+f[m]x^m
      g = g[0] + g[1]x+g[2]x^2+...+g[n]x^n
Wyjście h = f*g = h[0]+h[1]x+...+h[k]x^k

1. Jeżeli f = 0 lub g = 0, to    
   zwróć 0 i zakończ
2. Jeżeli oba wielomiany są stałe
   to zwróć f[0]*g[0] i zakończ
3. Przyjmij d := Ceil(max(deg f,deg g)/2)
4. Rozdziel oba wielomiany na części "górną" i "dolną"
   fl := f[0]+f[1]x+...+f[d-1]x^{d-1},  fu := f[d] + f[d+1]x+...+f[m]x^{m-d}
   gl := g[0]+g[1]x+...+g[d-1]x^{d-1},  gu := g[d] + g[d+1]x+...+g[n]x^{n-d}
5. Wywołaj algorytm rekurencyjnie i policz
   fl*gl , fu*gu , (fl + fu)(gl + gu)
6. Scal otrzymane wyniki
   h = fl*gl + ((fl + fu)(gl + gu) - fl*gl - fu*gu)x^d + fu*gu x^{2d}

Napisałem moduł do obsługi wielomianów zdefiniowanych następująco

Kopiuj
type TPolynomial = record
                     deg:integer;
                     coeff:array of Double;    
                   end;

Mam w nim mnożenie ale takie standardowe , te znane ze szkoły

Źródło tego pseudokodu to


gdzieś w okolicach 32:23

Problem z tym pseudokodem jest taki że wg mnie niepoprawnie dzieli te wielomiany na części dolną i górną
Gdybyśmy chcieli na podstawie tego pseudokodu napisać funkcję realizującą mnożenie wielomianów to mogłaby się zdarzyć sytuacja że
wyszlibyśmy indeksem poza tablicę przechowującą współczynniki

Jak ten pseudokod poprawić ?

...

Udostępnij
Komentuj
lion137
lion137
  • Rejestracja: dni
  • Ostatnio: dni
  • Postów: 5025
0

Możesz podać kontrprzykład dla punktu czwartego pseudokodu?

Udostępnij
Komentuj
nowy121105
N1
  • Rejestracja: dni
  • Ostatnio: dni
0

Dla wielomianów

Kopiuj
9x^4 + 7x^3 + 4x^2 + 8x +  9
8x^7+ 4x^6+ 2x^5 + x^4 + 6x^3 + 5x^2+ 9x+ 4

źle podzieli bo w pewnym momencie d będzie większe od stopnia któregoś z wielomianów f lub g.

Udostępnij
Komentuj
lion137
lion137
  • Rejestracja: dni
  • Ostatnio: dni
  • Postów: 5025
0

To chyba powinieneś uzupełniać zerami, bo i tak operujesz na tablicach.

Udostępnij
Komentuj
nowy121105
N1
  • Rejestracja: dni
  • Ostatnio: dni
0

Tak jak zasugerowałeś próbowałem uzupełniać zerami ale nic z tego nie wyszło

Moduł do obsługi wielomianów

Kopiuj

unit Polynomial;
interface
uses math;

const ZERO = 1e-10;

type TPolynomial = record
                     deg:integer;
                     coeff:array of Double
                   end;

function MakePolynomial(n:integer):TPolynomial;
function Equals(a,b:TPolynomial):boolean;
function AddPolynomials(a,b:TPolynomial):TPolynomial;
function SubtractPolynomials(a,b:TPolynomial):TPolynomial;
function MultiplyPolynomials(a,b:TPolynomial):TPolynomial;
function DividePolynomials(a,b:TPolynomial;var r:TPolynomial):TPolynomial;
function Horner(p:TPolynomial;x:Double):Double;
procedure WritePolynomial(p:TPolynomial;c:char);

implementation



function MakePolynomial(n:integer):TPolynomial;
var p:TPolynomial;
begin
  p.deg := n;
  SetLength(p.coeff,n+1);
  MakePolynomial := p
end;

procedure FindDeg(var p:TPolynomial);
var i:integer;
begin
  for i := p.deg downto 0 do
    if abs(p.coeff[i]) > ZERO then
    begin
      p.deg := i;
      Exit;
    end;
    p.deg := 0
end;

function AddPolynomials(a,b:TPolynomial):TPolynomial;
var newDeg,range,i:integer;
    n:TPolynomial;
begin
  if a.deg >  b.deg then
     newDeg := a.deg
  else
     newDeg := b.deg;
  range := newDeg - abs(a.deg - b.deg);
  n := MakePolynomial(newDeg);
  for i := range downto 0 do
     n.coeff[i] := a.coeff[i] + b.coeff[i];
  if a.deg > b.deg then
     for i := newDeg downto range+1 do
        n.coeff[i] := a.coeff[i]
  else
     for i := newDeg downto range+1 do
        n.coeff[i] := b.coeff[i];
  findDeg(n);
  AddPolynomials := n
end;

function SubtractPolynomials(a,b:TPolynomial):TPolynomial;
var newDeg,range,i:integer;
    n:TPolynomial;
begin
  if a.deg >  b.deg then
     newDeg := a.deg
  else
     newDeg := b.deg;
  range := newDeg - abs(a.deg - b.deg);
  n := MakePolynomial(newDeg);
  for i := range downto 0 do
     n.coeff[i] := a.coeff[i] - b.coeff[i];
  if a.deg > b.deg then
     for i := newDeg downto range+1 do
        n.coeff[i] := a.coeff[i]
  else
     for i := newDeg downto range+1 do
        n.coeff[i] := -b.coeff[i];
  findDeg(n);
  SubtractPolynomials := n
end;

function MultiplyPolynomials(a,b:TPolynomial):TPolynomial;
var i,j:integer;
    n:TPolynomial;
begin
   n := MakePolynomial(a.deg + b.deg);
   for i := 0 to n.deg do
     n.coeff[i] := 0;
   for i := 0 to a.deg do
      for j := 0 to b.deg do
         n.coeff[i + j] := n.coeff[i + j] + a.coeff[i] * b.coeff[j];
   MultiplyPolynomials := n
end;

function DividePolynomials(a,b:TPolynomial;var r:TPolynomial):TPolynomial;
var q,temp:TPolynomial;
    i,j,qdeg:integer;
begin
  if a.deg < b.deg then
     qdeg := 0
  else
     qdeg := a.deg;
  q := MakePolynomial(qdeg);
  if q.deg = 0 then
  begin
    q.coeff[0] := 0;
    r.deg := a.deg;
    for i :=r.deg downto 0 do
       r.coeff[i] := a.coeff[i];
       DividePolynomials := q;
       Exit
  end;
  temp := MakePolynomial(a.deg);
  for i := 0 to a.deg do
  begin
    temp.coeff[i] := a.coeff[i];
    q.coeff[i] := 0
  end;
  for i := a.deg - b.deg downto 0 do
  begin
    q.coeff[i] := temp.coeff[b.deg + i] /b.coeff[b.deg];
    for j := b.deg + i - 1 downto i do
       temp.coeff[j] := temp.coeff[j] - q.coeff[i] * b.coeff[j - i]
  end;
  for i := 0 to b.deg - 1 do
      r.coeff[i] := temp.coeff[i];
  if r.deg < a.deg then
     for i := b.deg to r.deg do
         r.coeff[i] := 0
     else
        for i := b.deg to a.deg do
           r.coeff[i] := 0;
     findDeg(r);
     findDeg(q);
     DividePolynomials := q
end;

function Equals(a,b:TPolynomial):boolean;
var bool : boolean;
    i : integer;
begin
  bool := a.deg = b.deg;
  i := a.deg;
  while bool and (i >= 0) do
  begin
    bool := bool and (a.coeff[i] = b.coeff[i]);
    Dec(i)
  end;
  Equals := bool
end;

procedure  WritePolynomial(p:TPolynomial;c:char);
var i:integer;
begin
  for i := p.deg downto 0 do
      if p.coeff[i] < 0 then
         write(p.coeff[i]:1:12,'*',c,'^',i)
      else
         write('+',p.coeff[i]:1:12,'*',c,'^',i);
  writeln
end;

function Horner(p:TPolynomial;x:Double):Double;
var v:double;
    i:integer;
begin
  v := p.coeff[p.deg];
  for i := p.deg - 1 downto 0 do
     v := v * x + p.coeff[i];
  Horner := v
end;

begin

end.

Kod programu do testowania funkcji realizującej algorytm Karacuby

Kopiuj



program untitled;

uses crt,math,polynomial;

function Karatsuba(f,g:TPolynomial):TPolynomial;
var temp1,temp2:TPolynomial;
    h,h1,h2,h3:TPolynomial;
    k,d:integer;
begin
  if ((f.deg = 0)and(f.coeff[0] = 0)) or ((g.deg = 0)and(g.coeff[0] = 0)) then
  begin
    h := makePolynomial(0);
    h.coeff[0] := 0;
    Karatsuba := h;
    Exit
  end;
  if ((f.deg = 0) and (g.deg = 0)) then
  begin
    h := makePolynomial(0);
    h.coeff[0] := f.coeff[0] * g.coeff[0];
    Karatsuba := h;
    Exit
  end;
  begin
    d := max(f.deg,g.deg);
    d := Ceil(0.5 * d);
    temp1 := makePolynomial(d-1);
    for k := 0 to min(d - 1,f.deg) do
        temp1.coeff[k] := f.coeff[k];
    for k := min(d-1,f.deg) + 1 to d - 1 do
        temp1.coeff[k] := 0;
    temp2 := makePolynomial(d-1);
    for k := 0 to min(d - 1,g.deg) do
        temp2.coeff[k] := g.coeff[k];
    for k := min(d - 1,g.deg) + 1 to d - 1 do
        temp2.coeff[k] := 0;
    h1 := Karatsuba(temp1,temp2);
    temp1 := makePolynomial(ord((f.deg - d)>=0)*(f.deg - d));
    for k := 0 to temp1.deg do
       temp1.coeff[k] := f.coeff[k + d];
    temp2 := makePolynomial(ord((g.deg - d)>=0)*(g.deg - d));
    for k := 0 to temp2.deg do
       temp2.coeff[k] := g.coeff[k+d];
    h2 := Karatsuba(temp1,temp2);
    temp1 := makePolynomial(d-1);
    temp2 := makePolynomial(ord((g.deg - d)>=0)*(g.deg - d));
    for k := 0 to min(d - 1,g.deg) do
       temp1.coeff[k] := g.coeff[k];
    for k := min(d - 1,g.deg) + 1 to d-1 do
        temp1.coeff[k] := 0;
    for k := 0 to temp2.deg do
       temp2.coeff[k] := g.coeff[k+d];
    h3 := AddPolynomials(temp1,temp2);
    for k := 0 to min(d - 1,f.deg) do
       temp1.coeff[k] := f.coeff[k];
    for k := min(d - 1,f.deg) + 1 to d - 1 do
       temp1.coeff[k] := 0;
    temp2 := makePolynomial(ord((f.deg - d)>=0)*(f.deg - d));
    for k := 0 to temp2.deg do
       temp2.coeff[k] := f.coeff[k+d];
    temp1 := AddPolynomials(temp1,temp2);
    h3 := Karatsuba(temp1,h3);
    temp1 := SubtractPolynomials(h3,h1);
    temp1 := SubtractPolynomials(temp1,h2);
    temp2 := makePolynomial(temp1.deg+d);
    for k := 0 to temp1.deg do
       temp2.coeff[k + d] := temp1.coeff[k];
    for k := 0 to d - 1 do
       temp2.coeff[k] := 0;
    h := AddPolynomials(h1,temp2);
    temp2 := makePolynomial(h2.deg + 2*d);
    for k := 0 to h2.deg do
        temp2.coeff[k + 2*d] := h2.coeff[k];
    for k := 0 to 2*d - 1 do
        temp2.coeff[k] := 0;
    h := AddPolynomials(h,temp2);
    Karatsuba := h;
  end;
end;

var f:text;
	k,n,m:integer;
	a,b,c:TPolynomial;

BEGIN
	
Assign(f,'input.txt');
Reset(f);
ReadLn(f,n);
a := makePolynomial(n);
k := n;
while not Eoln(f) do 
begin
Read(f,a.coeff[k]);
k:=k-1;
end;
ReadLn(f,m);
b := makePolynomial(m);
k := m;
while not Eoln(f) do 
begin
Read(f,b.coeff[k]);
k:=k-1;
end;
WriteLn('Polynomial a');
WriteLn(a.deg);
writePolynomial(a,'x');
WriteLn('Polynomial b');
WriteLn(b.deg);
writePolynomial(b,'x');
c := Karatsuba(a,b);
WriteLn('Polynomial c');
WriteLn(c.deg);
writePolynomial(c,'x');
Close(f);	
END.

Podobno błąd zakresu jest gdzieś w linii 42

Udostępnij
Komentuj
flowCRANE
flowCRANE
  • Rejestracja: dni
  • Ostatnio: dni
  • Lokalizacja: Tuchów
  • Postów: 12269
2

@nowy121105: postaw break point w linii 42 i sprawdź jakie wartości mają zmienne lokalne. Tylko w ten sposób dowiesz się dlaczego dostajesz range check error.

Udostępnij
Komentuj
nowy121105
N1
  • Rejestracja: dni
  • Ostatnio: dni
0
Kopiuj
program untitled;

uses crt,math,polynomial;

function Karatsuba(f,g:TPolynomial):TPolynomial;
var temp1,temp2:TPolynomial;
    h,h1,h2,h3:TPolynomial;
    k,m,n,d:integer;
begin
 //Base case f = 0 or g = 0
  if ((f.deg = 0)and(f.coeff[0] = 0)) or ((g.deg = 0)and(g.coeff[0] = 0)) then
  begin
    h := makePolynomial(0);
    h.coeff[0] := 0;
    Karatsuba := h;
    Exit
  end;
  //Base case both polynomials are constant
  if ((f.deg = 0) and (g.deg = 0)) then
  begin
    h := makePolynomial(0);
    h.coeff[0] := f.coeff[0] * g.coeff[0];
    Karatsuba := h;
    Exit
  end;
  begin
  //Checking which polynomial has maximum degree
    d := max(f.deg,g.deg); 
    //Calculating ceiling of half of max degree  
    d := (d + 1) div 2;
    //Creating temporary polynomial which stores lower part of polynomial f
    m := 1 + min(d - 1,f.deg);
    temp1 := makePolynomial(m-1);
    for k := 0 to m - 1 do
      temp1.coeff[k] := f.coeff[k];
     //Creating temporary polynomial which stores lower part of polynomial g
     n := 1 + min(d - 1,g.deg);
    temp2 := makePolynomial(n-1);
    for k := 0 to n - 1 do
       temp2.coeff[k] := g.coeff[k];
    //Recursive call of function Karatsuba to calculate product flow * glow
    h1 := Karatsuba(temp1,temp2);
    //Creating temporary polynomial which stores upper part of polynomial f
    temp1 := makePolynomial(ord(f.deg - m >= 0)*(f.deg - m));
    temp1.coeff[0] := 0;
    for k := 0 to f.deg - m do
       temp1.coeff[k] := f.coeff[k + m];
    //Creating temporary polynomial which stores upper part of polynomial g
    temp2 := makePolynomial(ord(g.deg - n >= 0)*(g.deg - n));
    temp2.coeff[0] := 0;
    for k := 0 to g.deg - n do
       temp2.coeff[k] := g.coeff[k+n];
    //Recursive call of function Karatsuba to calculate product fup * gup  
    h2 := Karatsuba(temp1,temp2);
    //Creating temporary polynomial to store lower part of polynomial g
    temp1 := makePolynomial(n-1);
    //Creating temporary polynomial to store upper part of polynomial g
    temp2 := makePolynomial(ord(g.deg - n >= 0)*(g.deg - n));
    for k := 0 to n - 1 do
       temp1.coeff[k] := g.coeff[k];
    temp2.coeff[0] := 0;
    for k := 0 to g.deg - n do
       temp2.coeff[k] := g.coeff[k+n];
    //Creating polynomial glow + gup
    h3 := AddPolynomials(temp1,temp2);
    //Set coefficients of temporary polynomial to coefficients of lower part of polynomial f
    for k := 0 to m - 1 do
       temp1.coeff[k] := f.coeff[k];
    //Creating temporary polynomial which stores upper part of polynomial f  
    temp2 := makePolynomial(ord(f.deg - m >= 0)*(f.deg - m));
    temp2.coeff[0] := 0;
    for k := 0 to f.deg - m do
       temp2.coeff[k] := f.coeff[k+m];
    //Creating polynomial flow + fup  
    temp1 := AddPolynomials(temp1,temp2);
    //Recursive call of function Karatsuba to calculate product (flow + fup)*(glow + gup)
    h3 := Karatsuba(temp1,h3);
    //Merging results , calculation middle terms
    temp1 := SubtractPolynomials(h3,h1);
    temp1 := SubtractPolynomials(temp1,h2);
    temp2 := makePolynomial(temp1.deg+d);
    for k := 0 to temp1.deg do
       temp2.coeff[k + d] := temp1.coeff[k];
    for k := 0 to d - 1 do
       temp2.coeff[k] := 0;
    //Adding lower part and middle part
    h := AddPolynomials(h1,temp2);
    //Creating upper part
    temp2 := makePolynomial(h2.deg + 2*d);
    for k := 0 to h2.deg do
        temp2.coeff[k + 2*d] := h2.coeff[k];
    for k := 0 to 2*d - 1 do
        temp2.coeff[k] := 0;
    //Adding upper part to current result
    h := AddPolynomials(h,temp2);
    Karatsuba := h;
  end;
end;

var f:text;
	k,n,m:integer;
	a,b,c:TPolynomial;

BEGIN
	
Assign(f,'input.txt');
Reset(f);
ReadLn(f,n);
a := makePolynomial(n);
k := n;
while not Eoln(f) do 
begin
Read(f,a.coeff[k]);
k:=k-1;
end;
ReadLn(f,m);
b := makePolynomial(m);
k := m;
while not Eoln(f) do 
begin
Read(f,b.coeff[k]);
k:=k-1;
end;
WriteLn('Polynomial a');
WriteLn(a.deg);
writePolynomial(a,'x');
WriteLn('Polynomial b');
WriteLn(b.deg);
writePolynomial(b,'x');
c := Karatsuba(a,b);
WriteLn('Polynomial c');
WriteLn(c.deg);
writePolynomial(c,'x');
Close(f);	
END.

Dla pliku tekstowego z danymi

Kopiuj
4
9 7 4 8 9
7
8 4 2 1 6 5 9 4

gdb po załadowaniu i uruchomieniu programu wyświetlił następujący komunikat

Kopiuj

(gdb) r
Starting program: F:\Karatsuba/testKaratsuba
[New Thread 16736.0x5bd0]
[New Thread 16736.0x23c8]
[New Thread 16736.0x6c4]
[New Thread 16736.0x59f0]
Polynomial a
4
+9.000000000000*x^4+7.000000000000*x^3+4.000000000000*x^2+8.000000000000*x^1+9.000000000000*x^0
Polynomial b
7
+8.000000000000*x^7+4.000000000000*x^6+2.000000000000*x^5+1.000000000000*x^4+6.000000000000*x^3+5.000000000000*x^2+9.000000000000*x^1+4.000000000000*x^0
Polynomial c
12
+20.000000000000*x^12+82.000000000000*x^11+92.000000000000*x^10+78.000000000000*x^9+83.000000000000*x^8+163.000000000000*x^7+143.000000000000*x^6+166.000000000000*x^5+176.000000000000*x^4+158.000000000000*x^3+133.000000000000*x^2+113.000000000000*x^1+36.000000000000*x^0
[Inferior 1 (process 16736) exited normally]
(gdb)

Czyli wygląda na to że uruchamia się bez błędów ale zwraca błędny wynik

Udostępnij
Komentuj
yarel
YA
  • Rejestracja: dni
  • Ostatnio: dni
  • Postów: 2384
0

Nie pamiętam już pascala, ale:

  • czy w MakePolynomail nie powinieneś jawnie zerować współczynników wielomianu? Nie wiem, czy w Pascalu jest szansa, że ta tablica będzie zawierać jakieś śmieciowe rzeczy i możesz dostać undefined behaviour.
  • w Karatsuba
Kopiuj
  if ((f.deg = 0)and(f.coeff[0] = 0)) or ((g.deg = 0)and(g.coeff[0] = 0)) then

coeff[0] jest double, więc może powinieneś porównywać numerycznie i mieć funkcję IsZero, która sprawdza czy wartość bezwzględna jest mniejsza niż jakaś ustalona delta (const ZERO = 1e-10;).

  • jakimś cudem wyszedł Ci stopień wielomianu x^12, przy x^7 * x^4 (czyli powinno być x^11) - jeśli nie chcesz bawić się debuggerem, to dodaj sobie najprostsze logowanie i zobacz, w którym miejscu tworzony jest stopień x^12.
Udostępnij
Komentuj
flowCRANE
flowCRANE
  • Rejestracja: dni
  • Ostatnio: dni
  • Lokalizacja: Tuchów
  • Postów: 12269
1
yarel napisał(a):

Nie wiem, czy w Pascalu jest szansa, że ta tablica będzie zawierać jakieś śmieciowe rzeczy i możesz dostać undefined behaviour.

W Pascalach, zmienne lokalne nie są inicjalizowane automatycznie (a więc zawierają śmieci), natomiast wszystkie zmienne globalne (deklarowane bezpośrednio w module, bez względu na sekcję) oraz wszystkie pola klas są inicjalizowane (zmienne globalne w trakcie inicjalizacji modułu, a pola klas w konstruktorach).

Udostępnij
Komentuj

Zarejestruj się i dołącz do największej społeczności programistów w Polsce.

Otrzymaj wsparcie, dziel się wiedzą i rozwijaj swoje umiejętności z najlepszymi.

Utwórz konto