Karatsuba algorithm for polynomial multiplication

[user07022024]

I found following pseudocode

Code: Text  [Select][+][-]

1. Karatsuba algorithm
2. Input: Polynomials
3.       f = f[0] + f[1]x + f[2]x^2+...+f[m]x^m
4.       g = g[0] + g[1]x + g[2]x^2+...+g[n]x^n
5. Output: Product of these polynomials
6.         h = fg = h[0]+h[1]x+...+h[k]x^k
7. 1. If f = 0 or g = 0 then
8.    return h=0 and exit
9. 2. If both polynomials are constant then
10.    return h = f[0]*g[0] and exit
11. 3. Let d := Ceil(max(deg f , deg g)/2)      
12. 4. Divide both polynomials to lower and upper parts
13.    flow = f[0]+f[1]x+...+f[d-1]x^{d-1} ; fup = f[d]+f[d+1]x+...+f[m]x^{m-d}
14.    glow = g[0]+g[1]x+...+g[d-1]x^{d-1} ; gup = g[d]+g[d+1]x+...+g[n]x^{n-d}
15. 5. Call algorithm recursively and calculate
16.    flow*glow , fup*gup, (flow + fup)(glow + gup)
17. 6. Merge received results
18.    h = flow*glow + ((flow + fup)*(glow + gup) - flow * glow - fup * gup)*x^{d} + fup*gup * x^{2d}
19.  

I wrote unit for polynomials with help of guy from math forum
and I tried to follow this pseudocode but I have got Range check error

Polynomial unit

Code: Pascal  [Select][+][-]

1. unit Polynomial;
2. interface
3. uses math;
4.  
5. const ZERO = 1e-10;
6.  
7. type TPolynomial = record
8.                      deg:integer;
9.                      coeff:array of Double
10.                    end;
11.  
12. function MakePolynomial(n:integer):TPolynomial;
13. function Equals(a,b:TPolynomial):boolean;
14. function AddPolynomials(a,b:TPolynomial):TPolynomial;
15. function SubtractPolynomials(a,b:TPolynomial):TPolynomial;
16. function MultiplyPolynomials(a,b:TPolynomial):TPolynomial;
17. function DividePolynomials(a,b:TPolynomial;var r:TPolynomial):TPolynomial;
18. function Horner(p:TPolynomial;x:Double):Double;
19. procedure WritePolynomial(p:TPolynomial;c:char);
20.  
21. implementation
22.  
23.  
24.  
25. function MakePolynomial(n:integer):TPolynomial;
26. var p:TPolynomial;
27. begin
28.   p.deg := n;
29.   SetLength(p.coeff,n+1);
30.   MakePolynomial := p
31. end;
32.  
33. procedure FindDeg(var p:TPolynomial);
34. var i:integer;
35. begin
36.   for i := p.deg downto 0 do
37.     if abs(p.coeff[i]) > ZERO then
38.     begin
39.       p.deg := i;
40.       Exit;
41.     end;
42.     p.deg := 0
43. end;
44.  
45. function AddPolynomials(a,b:TPolynomial):TPolynomial;
46. var newDeg,range,i:integer;
47.     n:TPolynomial;
48. begin
49.   if a.deg >  b.deg then
50.      newDeg := a.deg
51.   else
52.      newDeg := b.deg;
53.   range := newDeg - abs(a.deg - b.deg);
54.   n := MakePolynomial(newDeg);
55.   for i := range downto 0 do
56.      n.coeff[i] := a.coeff[i] + b.coeff[i];
57.   if a.deg > b.deg then
58.      for i := newDeg downto range+1 do
59.         n.coeff[i] := a.coeff[i]
60.   else
61.      for i := newDeg downto range+1 do
62.         n.coeff[i] := b.coeff[i];
63.   findDeg(n);
64.   AddPolynomials := n
65. end;
66.  
67. function SubtractPolynomials(a,b:TPolynomial):TPolynomial;
68. var newDeg,range,i:integer;
69.     n:TPolynomial;
70. begin
71.   if a.deg >  b.deg then
72.      newDeg := a.deg
73.   else
74.      newDeg := b.deg;
75.   range := newDeg - abs(a.deg - b.deg);
76.   n := MakePolynomial(newDeg);
77.   for i := range downto 0 do
78.      n.coeff[i] := a.coeff[i] - b.coeff[i];
79.   if a.deg > b.deg then
80.      for i := newDeg downto range+1 do
81.         n.coeff[i] := a.coeff[i]
82.   else
83.      for i := newDeg downto range+1 do
84.         n.coeff[i] := -b.coeff[i];
85.   findDeg(n);
86.   SubtractPolynomials := n
87. end;
88.  
89. function MultiplyPolynomials(a,b:TPolynomial):TPolynomial;
90. var i,j:integer;
91.     n:TPolynomial;
92. begin
93.    n := MakePolynomial(a.deg + b.deg);
94.    for i := 0 to n.deg do
95.      n.coeff[i] := 0;
96.    for i := 0 to a.deg do
97.       for j := 0 to b.deg do
98.          n.coeff[i + j] := n.coeff[i + j] + a.coeff[i] * b.coeff[j];
99.    MultiplyPolynomials := n
100. end;
101.  
102. function DividePolynomials(a,b:TPolynomial;var r:TPolynomial):TPolynomial;
103. var q,temp:TPolynomial;
104.     i,j,qdeg:integer;
105. begin
106.   if a.deg < b.deg then
107.      qdeg := 0
108.   else
109.      qdeg := a.deg;
110.   q := MakePolynomial(qdeg);
111.   if q.deg = 0 then
112.   begin
113.     q.coeff[0] := 0;
114.     r.deg := a.deg;
115.     for i :=r.deg downto 0 do
116.        r.coeff[i] := a.coeff[i];
117.        DividePolynomials := q;
118.        Exit
119.   end;
120.   temp := MakePolynomial(a.deg);
121.   for i := 0 to a.deg do
122.   begin
123.     temp.coeff[i] := a.coeff[i];
124.     q.coeff[i] := 0
125.   end;
126.   for i := a.deg - b.deg downto 0 do
127.   begin
128.     q.coeff[i] := temp.coeff[b.deg + i] /b.coeff[b.deg];
129.     for j := b.deg + i - 1 downto i do
130.        temp.coeff[j] := temp.coeff[j] - q.coeff[i] * b.coeff[j - i]
131.   end;
132.   for i := 0 to b.deg - 1 do
133.       r.coeff[i] := temp.coeff[i];
134.   if r.deg < a.deg then
135.      for i := b.deg to r.deg do
136.          r.coeff[i] := 0
137.      else
138.         for i := b.deg to a.deg do
139.            r.coeff[i] := 0;
140.      findDeg(r);
141.      findDeg(q);
142.      DividePolynomials := q
143. end;
144.  
145. function Equals(a,b:TPolynomial):boolean;
146. var bool : boolean;
147.     i : integer;
148. begin
149.   bool := a.deg = b.deg;
150.   i := a.deg;
151.   while bool and (i >= 0) do
152.   begin
153.     bool := bool and (a.coeff[i] = b.coeff[i]);
154.     Dec(i)
155.   end;
156.   Equals := bool
157. end;
158.  
159. procedure  WritePolynomial(p:TPolynomial;c:char);
160. var i:integer;
161. begin
162.   for i := p.deg downto 0 do
163.       if p.coeff[i] < 0 then
164.          write(p.coeff[i]:1:12,'*',c,'^',i)
165.       else
166.          write('+',p.coeff[i]:1:12,'*',c,'^',i);
167.   writeln
168. end;
169.  
170. function Horner(p:TPolynomial;x:Double):Double;
171. var v:double;
172.     i:integer;
173. begin
174.   v := p.coeff[p.deg];
175.   for i := p.deg - 1 downto 0 do
176.      v := v * x + p.coeff[i];
177.   Horner := v
178. end;
179.  
180. begin
181.  
182. end.
183.  
184.  

Program for testing function for Karatsuba multiplication

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

How to use command line gdb to find error

[srvaldez]

what's in the file "input.txt" ?

[user07022024]

First line: degree of polynomial a
Second line: coefficients of polynomial a
Third line: degree of polynomial b
Fourth line: coefficients of polynomial b

Here is example of input file

Code: Text  [Select][+][-]

1. 4
2. 9 7 4 8 9
3. 7
4. 8 4 2 1 6 5 9 4
5.  

[Zvoni]

At a guess: You're reading a String from the file, and assign it to (array of) Double
I would first read the file, parse it (Split the Line for the coefficients), convert everything to its correct Datatype, and then call your functions

[MathMan]

Hello user07022024,

I'm sorry but I never used gdb from the command line, so I cannot help you there.

I quickly looked at your sources and I think you need to change the way you split the polynomes. To me it seems like the high part generates the range error due to access to the f(i), g(i) out of bounds. But I may be wrong here - didn't test.

Generally Karatsuba (or similar divide & conquer) multiplication is used for polynomes of equal order. It might already help if you extend the order of the 'smaller' polynome to the order of the 'larger' polynome by prepending sufficient number of leading 0 coefficients. Also better handle the split size via integer calculation - t.i. instead of

Code: Pascal  [Select][+][-]

1.     //Calculating ceiling of half of max degree  
2.     d := Ceil(0.5 * d);
3.  

do

Code: Pascal  [Select][+][-]

1.     d := ( d+1 ) div 2;
2.  

Some general questions:

- have you tried existing libraries i.e. FLINT - https://flintlib.org ?
- you define the coefficients f(i), g(i) as float. Not impossible but very unusual as it can/will generate rounding errors. Is there a reason?
- what is your objective with this? Do you just want to learn how this works? If not then consider my first question.

Cheers,
MathMan

[440bx]

@user07022024,

Is there a reason you cannot (or do not want to) use FpDebug ?

[mizar]

@MathMan
is there any flintlib bindings for fpc ?

[Fred vS]

Hello Mizar.
Here how to do to use gdb under Linux (for Windows, just use gdb.exe and "\" as directory separator):
All is done via a terminal.
First compile your demo with debug data:

Code: Bash  [Select][+][-]

1. > cd /the_directory_of_project/
2. > fpc -ountitled -l -Mobjfpc -Sh -Fcutf8 -gl -O- -B untitled.pas

Then load the debugger, still in the directory of the project:

Code: Bash  [Select][+][-]

1. > gdb untitled

Once gdb is loaded, press the "r" character (for run) + "enter"

You should get something like this:

(gdb) r
Starting program: /home/fred/Polynomial_test/untitled
Polynomial a
4
+9.000000000000*x^4+7.000000000000*x^3+4.000000000000*x^2+8.000000000000*x^1+9.0
00000000000*x^0
Polynomial b
7
+8.000000000000*x^7+4.000000000000*x^6+2.000000000000*x^5+1.000000000000*x^4+6.0
00000000000*x^3+5.000000000000*x^2+9.000000000000*x^1+4.000000000000*x^0
Polynomial c
11
+72.000000000000*x^11+92.000000000000*x^10+78.000000000000*x^9+103.000000000000*
x^8+173.000000000000*x^7+143.000000000000*x^6+166.000000000000*x^5+176.000000000
000*x^4+158.000000000000*x^3+133.000000000000*x^2+113.000000000000*x^1+36.000000
000000*x^0
[Inferior 1 (process 576738) exited normally]
(gdb)

Note that I did not get any error.
I just changed this from your code:

Code: Pascal  [Select][+][-]

1.  var f:textfile;
2.  Assignfile(f,'input.txt');

 

[MathMan]

@MathMan
is there any flintlib bindings for fpc ?

I don't think there are - the only FPC bindings i know are for GMP.

But again - the question is what do you want to achieve?

[mizar]

@MathMan:

I'm  not the O.P., just read about flintlib (that I didn't know) in your reply and was just curious to know if it was usable from fpc 

[MathMan]

@mizar

Sorry for the confusion on my side. Re FLINT - it has a C interface (not C++!), so it shouldn't be too complex to create a binding for FPC. I searched a bit and it looks like there really is nothing available as of yet.

You can earn a name if you'd create one  But FLINT is really large so you may need some automated thing to start this off. The library is, if not the best, among the top three of its kind - many of the best in computer number theory / computer algebra have contributed to it over time.

[Laksen]

How big do you think integer is?

I would advice using a type that is not dependent on platform, like longint

[TRon]

How big do you think integer is?

I would advice using a type that is not dependent on platform, like longint

🤣

[Handoko]

I just want to add, longint is not always platform/cpu independent at least in Delphi.

In newer Delphi versions, the longint type is platform and CPU dependent.

Source: https://www.freepascal.org/docs-html/ref/refsu4.html

[tetrastes]

There is logical error in your recursion. At some step we have f.deg=0 and g.deg=3, so we are at line 10 in code below with d=2:

Code: Pascal  [Select][+][-]

1.    if f.deg > g.deg then
2.       d := f.deg
3.    else
4.       d := g.deg;
5.    //Calculating ceiling of half of max degree
6.    d := Ceil(0.5 * d);
7.    //Creating temporary polynomial which stores lower part of polynomial f
8.    temp1 := makePolynomial(d-1);
9.    for k := 0 to d - 1 do
10.      temp1.coeff[k] := f.coeff[k];

but there is no f.coef[1] obviously, as f.deg=0, thus we have Range check error.

I changed in this line and to or:

Code: Pascal  [Select][+][-]

1.  //Base case both polynomials are constant
2.  if ((f.deg = 0) or{and} (g.deg = 0)) then

and got the following result without errors, though I don't know if it correct, of course:

C:\work\lazarus\tests\Karatsuba>project1.exe
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
10                                                                                                                                                                                                                 +56.000000000000*x^10+60.000000000000*x^9+103.000000000000*x^8+245.000000000000*x^7+179.000000000000*x^6+184.000000000000*x^5+176.000000000000*x^4+158.000000000000*x^3+133.000000000000*x^2+113.000000000000*x^1+36.000000000000*x^0