Recent

Author Topic: Karatsuba algorithm for polynomial multiplication  (Read 2482 times)

MathMan

  • Sr. Member
  • ****
  • Posts: 376
Re: Karatsuba algorithm for polynomial multiplication
« Reply #30 on: October 27, 2024, 10:32:02 pm »
My school age passed away nearly quarter century ago
so it is not a homework

I then apologise for my assumption and hope you were not offended by it.

I found this pseudocode on youtube and tried to follow it or at least stay as close as possible to it

If you had stated that when I initially asked for the purpose/intention of it, I'd have acted a bit different. So I understand now that this was a learning experience.

I hope that nonetheless I could help you a bit ...

Cheers,
MathMan

user07022024

  • New Member
  • *
  • Posts: 13
Re: Karatsuba algorithm for polynomial multiplication
« Reply #31 on: October 28, 2024, 10:43:05 am »
Your function works only in objfpc mode
and if i declare Result variable in local variable section then add Karatsuba := Result
then  gdb prints Range check error

In objfpc mode it is better to use classes
but I had only procedural Pascal and many years ago

Your function would be good if we have generic class for Polynomials
(generic to be able to choose type of coefficients)
or even better if we apply improvements which you mentioned like changing base case etc

Why , when I remove objfpc mode and replace Result variable with other variable I get range check error
« Last Edit: October 28, 2024, 11:15:16 am by user07022024 »

MathMan

  • Sr. Member
  • ****
  • Posts: 376
Re: Karatsuba algorithm for polynomial multiplication
« Reply #32 on: October 28, 2024, 12:33:28 pm »
Your function works only in objfpc mode
and if i declare Result variable in local variable section then add Karatsuba := Result
then  gdb prints Range check error

In objfpc mode it is better to use classes
but I had only procedural Pascal and many years ago

Your function would be good if we have generic class for Polynomials
(generic to be able to choose type of coefficients)
or even better if we apply improvements which you mentioned like changing base case etc

Why , when I remove objfpc mode and replace Result variable with other variable I get range check error

What exactly did you do regarding 'ObjFpc' mode? Did you replace it with a different mode, and if so which? Could it be that you are using 'Turbo Pasccal' mode?

If you want to remove the 'Result' variable, then simply replace 'Result' with 'Karatsuba' in the function. No need to add a 'Result' variable to the var section of the function and assign 'Karatsuba := Result' at the end.

Why GDB throws a range error depends on the actual mode you use for compilation.

PS - I have to admit that there still was an error in my function - I forgot to include the special case that one of the polys is '0'. Should have tested more thoroughly  :-[ Here is the corrected version

Code: Pascal  [Select][+][-]
  1. function Karatsuba(
  2.  f: TPolynomial;
  3.  g: TPolynomial
  4. ):TPolynomial;
  5.  
  6. var
  7.   Count: NativeInt; // generic loop control
  8.   Split: NativeInt; // degree of polynome split
  9.  
  10.   pl: tPolynomial;  // larger degree polynomial
  11.   ps: tPolynomial;  // smaller degree polynomial
  12.  
  13.   plHi, plLo: tPolynomial; // high & low part of larger polynome
  14.   psHi, psLo: tPolynomial; // same for smaller polynome
  15.  
  16.   pHi, pMid, pLo: tPolynomial; // partial results of Karatsuba
  17.  
  18. begin
  19.   // special case:
  20.   // - at least one of the polynomes is '0'
  21.   // - the result polynome then is '0' too
  22.   if( ( ( f.Deg=0 ) and ( f.Coeff[ 0 ]=0 ) ) or
  23.       ( ( g.Deg=0 ) and ( g.Coeff[ 0 ]=0 ) ) ) then begin
  24.     Result := MakePolynomial( 0 );
  25.     Result.Coeff[ 0 ] := 0;
  26.     Exit;
  27.   end;
  28.  
  29.   // base case:
  30.   // - both polinomes are constant
  31.   // - the polynome product is the product of constants
  32.   if( ( f.Deg=0 ) and ( g.Deg=0 ) ) then begin
  33.     Result := MakePolynomial( 0 );
  34.     Result.Coeff[ 0 ] := f.Coeff[ 0 ] * g.Coeff[ 0 ];
  35.     Exit;
  36.   end;
  37.  
  38.   // sort the polynomes by degree to simplify the following calculations
  39.   if( f.Deg>=g.Deg ) then begin
  40.     pl := f;
  41.     ps := g;
  42.   end else begin
  43.     pl := g;
  44.     ps := f;
  45.   end;
  46.  
  47.   // get the split degree
  48.   Split := pl.Deg div 2;
  49.  
  50.   // now split the polynomials into high & low parts
  51.   // the split of the larger polynome is trivial, as high & low part exist
  52.   plLo := MakePolynomial( Split );
  53.   for Count := 0 to Split do plLo.Coeff[ Count ] := pl.Coeff[ Count ];
  54.   plHi := MakePolynomial( pl.Deg-Split-1 );
  55.   for Count := Split+1 to pl.Deg do plHi.Coeff[ Count-Split-1 ] := pl.Coeff[ Count ];
  56.  
  57.   // for the smaller polynome we have to check, if the high part exists
  58.   psLo := MakePolynomial( Min( Split, ps.Deg ) );
  59.   for Count := 0 to Min( Split, ps.Deg ) do psLo.Coeff[ Count ] := ps.Coeff[ Count ];
  60.   if( Split>=ps.Deg ) then begin
  61.     psHi := MakePolynomial( 0 );
  62.     psHi.Coeff[ 0 ] := 0;
  63.   end else begin
  64.     psHi := MakePolynomial( ps.Deg-Split-1 );
  65.     for Count := Split+1 to ps.Deg do psHi.Coeff[ Count-Split-1 ] := ps.Coeff[ Count ];
  66.   end;
  67.  
  68.   // calculate low & high product
  69.   pLo := Karatsuba( plLo, psLo );
  70.   pHi := Karatsuba( plHi, psHi );
  71.  
  72.   // calculate middle product ( plHi+plLo )*( psHi+psLo )
  73.   // we can re-use low parts for the sums
  74.   plLo := AddPolynomials( plLo, plHi );
  75.   psLo := AddPolynomials( psLo, psHi );
  76.   pMid := Karatsuba( plLo, psLo );
  77.  
  78.   Result := MakePolynomial( f.Deg+g.Deg );
  79.  
  80.   // finally sum all parts up, positionally correct
  81.   for Count := 0 to pLo.Deg do
  82.     Result.Coeff[ Count ] := pLo.Coeff[ Count ];
  83.   for Count := 0 to pHi.Deg do
  84.     Result.Coeff[ Result.Deg-pHi.Deg+Count ] := pHi.Coeff[ Count ];
  85.   for Count := 0 to pMid.Deg do
  86.     Result.Coeff[ Count+Split+1 ] := Result.Coeff[ Count+Split+1 ] + pMid.Coeff[ Count ];
  87.   for Count := 0 to pLo.Deg do
  88.     Result.Coeff[ Count+Split+1 ] := Result.Coeff[ Count+Split+1 ] - pLo.Coeff[ Count ];
  89.   for Count := 0 to pHi.Deg do
  90.     Result.Coeff[ Count+Split+1 ] := Result.Coeff[ Count+Split+1 ] - pHi.Coeff[ Count ];
  91. end;
  92.  
« Last Edit: October 28, 2024, 01:04:05 pm by MathMan »

MathMan

  • Sr. Member
  • ****
  • Posts: 376
Re: Karatsuba algorithm for polynomial multiplication
« Reply #33 on: October 30, 2024, 01:48:30 pm »
@user07022024

Out of curiosity I did some benchmarking of the version from my last post against the 'schoolbook' variant MultiplyPolynomials to get an insight on the efficiency. I initially said that the approach is not that efficient, but I was astonished how inefficient it actually is.

The break even between the 'schoolbook' and the Karatsuba variant is around polynomial degree 65,000 and 500,000 - depending on compiler settings.

So I did another implementation, the way I would do it (see below). Here the break even is around 64 to 512 - again depending on compiler settings.

Code: Pascal  [Select][+][-]
  1. procedure KaraRec(
  2.   Prod: pDouble;
  3.   Poly1: pDouble;
  4.   Poly2: pDouble;
  5.   Size: NativeInt;
  6.   Heap: pDouble );
  7.  
  8. var
  9.   Count: NativeInt;
  10.  
  11.   Hi: NativeInt;
  12.   Lo: NativeInt;
  13.  
  14.   Sum1: pDouble;
  15.   Sum2: pDouble;
  16.  
  17. begin
  18.   // check if split is required
  19.   if( Size>1 ) then begin
  20.     // split with lower part as shorter - doesn't matter anyway
  21.     Lo := Size shr 1;
  22.     Hi := Size - Lo;  // Hi is in range [ Lo..Lo+1 ]
  23.  
  24.     // recurse for products of lower / higher parts and put on Heap
  25.     KaraRec( Heap, Poly1, Poly2, Lo, Heap+2*Lo );
  26.     KaraRec( Heap+2*Lo, Poly1+Lo, Poly2+Lo, Hi, Heap+2*Size );
  27.  
  28.     // calculate the sums of lower / higher part for both polynomials
  29.     Sum1 := Prod;
  30.     Sum2 := Prod + Hi;
  31.  
  32.     for Count := 0 to Lo-1 do
  33.     begin
  34.       Sum1[ Count ] := Poly1[ Count ] + Poly1[ Lo+Count ];
  35.       Sum2[ Count ] := Poly2[ Count ] + Poly2[ Lo+Count ];
  36.     end;
  37.     if( Hi<>Lo ) then
  38.     begin
  39.       Sum1[ Lo ] := Poly1[ Size-1 ];
  40.       Sum2[ Lo ] := Poly2[ Size-1 ];
  41.     end;
  42.  
  43.     // recurse for product of sums
  44.     KaraRec( Heap+2*Size, Sum1, Sum2, Hi, Heap+2*( Size+Hi ) );
  45.  
  46.     // re-assemble, positionally correct
  47.     for Count := 0 to 2*Lo-1 do
  48.       Prod[ Count ] := Heap[ Count ];
  49.     for Count := 0 to 2*Hi-1 do
  50.       Prod[ 2*Lo+Count ] := Heap[ 2*Lo+Count ];
  51.     for Count := 0 to 2*Hi-1 do begin
  52.       Prod[ Lo+Count ] := Prod[ Lo+Count ] + Heap[ 2*Size+Count ];
  53.     end;
  54.     for Count := 0 to 2*Lo-1 do
  55.       Prod[ Lo+Count ] := Prod[ Lo+Count ] - Heap[ Count ];
  56.     for Count := 0 to 2*Hi-1 do
  57.       Prod[ Lo+Count ] := Prod[ Lo+Count ] - Heap[ 2*Lo+Count ];
  58.   end else begin
  59.     // both polynomials are constant
  60.     // calculate the product and virtually extend with 0*x^1, size regime!
  61.     Prod[ 0 ] := Poly1[ 0 ] * Poly2[ 0 ];
  62.     Prod[ 1 ] := 0;
  63.   end;
  64. end;
  65.  
  66. procedure Karatsuba(
  67.   p: tPolynomial;
  68.   f: tPolynomial;
  69.   g: tPolynomial );
  70.  
  71. var
  72.   Count: NativeInt;
  73.   Size: NativeInt;
  74.  
  75.   Prod: pDouble;
  76.   Poly1: pDouble;
  77.   Poly2: pDouble;
  78.   Heap: pDouble;
  79.  
  80. begin
  81.   // size of bare Karatsuba must be one larger than the tPolynomial degree
  82.   // as the recursive KaraRec works on floating point arrays!
  83.   Size := Max( f.Deg, g.Deg ) + 1;
  84.  
  85.   // organise buffers
  86.   GetMem( Poly1, SizeOf( Double )*Size   );
  87.   GetMem( Poly2, SizeOf( Double )*Size   );
  88.   GetMem( Prod,  SizeOf( Double )*2*Size );
  89.   GetMem( Heap,  SizeOf( Double )*8*Size );
  90.  
  91.   // transfer polynomial coefficients to arrays, equalising sizes
  92.   for Count := 0 to f.Deg do Poly1[ Count ] := f.Coeff[ Count ];
  93.   for Count := f.Deg+2 to Size do Poly1[ Count ] := 0;
  94.   for Count := 0 to g.Deg do Poly2[ Count ] := g.Coeff[ Count ];
  95.   for Count := g.Deg+2 to Size do Poly2[ Count ] := 0;
  96.  
  97.   // calculate the product
  98.   KaraRec( Prod, Poly1, Poly2, Size, Heap );
  99.  
  100.   // transfer back from array to polynomial type
  101.   // Attn:
  102.   //   the resulting product polynomial p must be initialised a priori to a
  103.   //   sufficiently large degree - f.Deg+g.Deg at least. The polynomial could
  104.   //   also be instantiated here, but this would cause memory fragmentation, as
  105.   //   the buffers can only be released after the copy!
  106.   for Count := 0 to f.Deg+g.Deg do p.Coeff[ Count ] := Prod[ Count ];
  107.  
  108.   // release buffers
  109.   FreeMem( Heap,  SizeOf( Double )*8*Size );
  110.   FreeMem( Prod,  SizeOf( Double )*2*Size );
  111.   FreeMem( Poly2, SizeOf( Double )*Size   );
  112.   FreeMem( Poly1, SizeOf( Double )*Size   );
  113. end;
  114.  

user07022024

  • New Member
  • *
  • Posts: 13
Re: Karatsuba algorithm for polynomial multiplication
« Reply #34 on: October 30, 2024, 02:58:41 pm »
I changed base case as user tetrastes suggested and
at first look seemed to work but it still is someting wrong
Base case is a little bit more complicated but it should be less recursive calls

In my opinion the original base case should also work but with more recursive calls


Code: Pascal  [Select][+][-]
  1.  
  2.  
  3. program untitled;
  4.  
  5. uses crt,math,polynomial;
  6.  
  7. function Karatsuba(f,g:TPolynomial):TPolynomial;
  8. var temp1,temp2:TPolynomial;
  9.     h,h1,h2,h3:TPolynomial;
  10.     k,m,n,d:integer;
  11. begin
  12.   if f.deg < g.deg then
  13.   begin
  14.     temp1 := f;
  15.     f := g;
  16.     g := temp1;
  17.   end;
  18.  //Base case f = 0 or g = 0
  19.   if ((f.deg = 0)and(f.coeff[0] = 0)) or ((g.deg = 0)and(g.coeff[0] = 0)) then
  20.   begin
  21.     h := makePolynomial(0);
  22.     h.coeff[0] := 0;
  23.     Karatsuba := h;
  24.     Exit
  25.   end;
  26.   //Base case both polynomials are constant
  27.   if g.deg = 0 then
  28.   begin
  29.     h := makePolynomial(f.deg);
  30.     for k := 0 to h.deg do
  31.         h.coeff[k] := g.coeff[0] * f.coeff[k];
  32.     Karatsuba := h;
  33.     Exit
  34.   end;
  35.   begin
  36.   //Checking which polynomial has maximum degree
  37.     d := max(f.deg,g.deg);
  38.     //Calculating ceiling of half of max degree  
  39.     d := (d + 1) div 2;
  40.     //Creating temporary polynomial which stores lower part of polynomial f
  41.     m := min(d ,f.deg + 1);
  42.     temp1 := makePolynomial(m-1);
  43.     for k := 0 to m - 1 do
  44.       temp1.coeff[k] := f.coeff[k];
  45.      //Creating temporary polynomial which stores lower part of polynomial g
  46.      n := min(d,g.deg + 1);
  47.     temp2 := makePolynomial(n-1);
  48.     for k := 0 to n - 1 do
  49.        temp2.coeff[k] := g.coeff[k];
  50.     //Recursive call of function Karatsuba to calculate product flow * glow
  51.     h1 := Karatsuba(temp1,temp2);
  52.     //Creating temporary polynomial which stores upper part of polynomial f
  53.     temp1 := makePolynomial(ord(f.deg - m >= 0)*(f.deg  - m));
  54.     temp1.coeff[0] := 0;
  55.     for k := 0 to f.deg - m do
  56.        temp1.coeff[k] := f.coeff[k + m];
  57.     //Creating temporary polynomial which stores upper part of polynomial g
  58.     temp2 := makePolynomial(ord(g.deg - n >= 0)*(g.deg  - n));
  59.     temp2.coeff[0] := 0;
  60.     for k := 0 to g.deg - n do
  61.        temp2.coeff[k] := g.coeff[k+n];
  62.     //Recursive call of function Karatsuba to calculate product fup * gup  
  63.     h2 := Karatsuba(temp1,temp2);
  64.     //Creating temporary polynomial to store lower part of polynomial g
  65.     temp1 := makePolynomial(n-1);
  66.     for k:= 0 to n - 1 do
  67.         temp1.coeff[k] := g.coeff[k];
  68.     //Creating temporary polynomial to store upper part of polynomial g
  69.     temp2 := makePolynomial(ord(g.deg - n >= 0)*(g.deg  - n));
  70.     for k := 0 to n - 1 do
  71.        temp1.coeff[k] := g.coeff[k];
  72.     temp2.coeff[0] := 0;
  73.     for k := 0 to g.deg - n do
  74.        temp2.coeff[k] := g.coeff[k+n];
  75.     //Creating polynomial glow + gup
  76.     h3 := AddPolynomials(temp1,temp2);
  77.     //Set coefficients of temporary polynomial to coefficients of lower part of polynomial f
  78.     for k := 0 to m - 1 do
  79.        temp1.coeff[k] := f.coeff[k];
  80.     //Creating temporary polynomial which stores upper part of polynomial f  
  81.     temp2 := makePolynomial(ord(f.deg - m >= 0)*(f.deg  - m));
  82.     temp2.coeff[0] := 0;
  83.     for k := 0 to f.deg - m do
  84.        temp2.coeff[k] := f.coeff[k+m];
  85.     //Creating polynomial flow + fup  
  86.     temp1 := AddPolynomials(temp1,temp2);
  87.     //Recursive call of function Karatsuba to calculate product (flow + fup)*(glow + gup)
  88.     h3 := Karatsuba(temp1,h3);
  89.     //Merging results , calculation middle terms
  90.     d := min(m,n);
  91.     temp1 := SubtractPolynomials(h3,h1);
  92.     temp1 := SubtractPolynomials(temp1,h2);
  93.     temp2 := makePolynomial(temp1.deg+d);
  94.     for k := 0 to temp1.deg do
  95.        temp2.coeff[k + d] := temp1.coeff[k];
  96.     for k := 0 to d - 1 do
  97.        temp2.coeff[k] := 0;
  98.     //Adding lower part and middle part
  99.     h := AddPolynomials(h1,temp2);
  100.     //Creating upper part
  101.     temp2 := makePolynomial(h2.deg + 2*d);
  102.     for k := 0 to h2.deg do
  103.         temp2.coeff[k + 2*d] := h2.coeff[k];
  104.     for k := 0 to 2*d - 1 do
  105.         temp2.coeff[k] := 0;
  106.     //Adding upper part to current result
  107.     h := AddPolynomials(h,temp2);
  108.     Karatsuba := h;
  109.   end;
  110. end;
  111.  
  112. var f:text;
  113.         k,n,m:integer;
  114.         a,b,c:TPolynomial;
  115.  
  116. BEGIN
  117.        
  118. Assign(f,'input.txt');
  119. Reset(f);
  120. ReadLn(f,n);
  121. a := makePolynomial(n);
  122. k := n;
  123. while not Eoln(f) do
  124. begin
  125. Read(f,a.coeff[k]);
  126. k:=k-1;
  127. end;
  128. ReadLn(f,m);
  129. b := makePolynomial(m);
  130. k := m;
  131. while not Eoln(f) do
  132. begin
  133. Read(f,b.coeff[k]);
  134. k:=k-1;
  135. end;
  136. WriteLn('Polynomial a');
  137. WriteLn(a.deg);
  138. writePolynomial(a,'x');
  139. WriteLn('Polynomial b');
  140. WriteLn(b.deg);
  141. writePolynomial(b,'x');
  142. c := Karatsuba(a,b);
  143. WriteLn('Polynomial c');
  144. WriteLn(c.deg);
  145. writePolynomial(c,'x');
  146. Close(f);      
  147. END.
  148.  

Still is something wrong in this code
« Last Edit: October 30, 2024, 04:32:33 pm by user07022024 »

 

TinyPortal © 2005-2018