1/*
2 * toomcook - implementation of Toom-Cook(3,4) multiplication algorithm
3 *
4 * Copyright (C) 2013,2021 Christoph Zurnieden
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	2013/08/11 01:31:28
21 * File existed as early as:	2013
22 */
23
24/*
25 * hide internal function from resource debugging
26 */
27static resource_debug_level;
28resource_debug_level = config("resource_debug", 0);
29
30
31/* */
32define toomcook3(a,b){
33  local alen blen  a0 a1 a2  b0 b1 b2 m ret sign mask;
34  local S0 S1 S2 S3 S4 T1 T2;
35
36  if(!isint(a) || !isint(b))
37      return newerror("toomcook3(a,b): a and/or b is not an integer");
38
39  alen = digits(a,2);
40  blen = digits(b,2);
41
42  sign = sgn(a) * sgn(b);
43  /* sgn(x) returns 0 if x = 0 */
44  if(sign == 0) return 0;
45
46  m = min(alen,blen)//3;
47  mask = ~-(1<<m);
48
49  /*
50    Cut-off at about 4,000 dec. digits
51    TODO: check
52  */
53  if(isdefined("test8900")){
54    if(m < 20) return a*b;
55  }
56  else{
57    if(m < 4096 ) return a*b;
58  }
59  a = abs(a);
60  b = abs(b);
61
62  a0 = a & mask;
63  a1 = (a>>m) & mask;
64  a2 = (a>>(2*m));
65
66  b0 =  b         & mask;
67  b1 = (b>>m)     & mask;
68  b2 = (b>>(2*m));
69
70  /*
71    Zimmermann
72  */
73
74  S0 = toomcook3(a0 , b0);
75  S1 = toomcook3((a2+a1+a0) , (b2+b1+b0));
76  S2 = toomcook3(((a2<<2)+(a1<<1)+a0) , ((b2<<2)+(b1<<1)+b0));
77  S3 = toomcook3((a2-a1+a0) , (b2-b1+b0));
78  S4 = toomcook3(a2,b2);
79  T1 = (S3<<1) + S2;
80  T1 /= 3;
81  T1 += S0;
82  T1 >>= 1;
83  T1 -= S4<<1;
84  T2 = (S1 + S3)>>1;
85  S1 -= T1;
86  S2 = T2 - S0 - S4;
87  S3 = T1 - T2;
88
89  ret =  (S4<<(4*m)) + (S3<<(3*m)) + (S2<<(2*m)) + (S1<<(1*m)) + S0;
90
91
92  ret = sign *ret;
93
94  return ret;
95}
96
97define toomcook3square(a){
98  local alen  a0 a1 a2  m tmp tmp2 ret sign S0 S1 S2 S3 S4 T1 mask;
99
100  if(!isint(a))return newerror("toomcook3square(a): a is not integer");
101
102  alen = digits(a,2);
103
104  sign = sgn(a) * sgn(a);
105  if(sign == 0) return 0;
106
107  m = alen//3;
108  mask = ~-(1<<m);
109  /*
110    Cut-off at about 5,000 dec. digits
111    TODO: check
112  */
113
114  if(isdefined("test8900")){
115    if(m < 20) return a^2;
116  }
117  else{
118    if(m < 5000 ) return a^2;
119  }
120
121  a = abs(a);
122
123  a0 = a & mask;
124  a1 = (a>>m) & mask;
125  a2 = (a>>(2*m));
126
127  /*
128   Bodrato/Zanoni
129   */
130  S0 = toomcook3square(a0);
131  S1 = toomcook3square(a2+a1+a0);
132  S2 = toomcook3square(a2-a1+a0);
133  S3 = toomcook3(a1<<1,a2);
134  S4 = toomcook3square(a2);
135
136  T1 = (S1 + S2)>>1;
137  S1 = S1 - T1 - S3;
138  S2 = T1 - S4 -S0;
139
140
141  S1 = S1<<(1*m);
142  S2 = S2<<(2*m);
143  S3 = S3<<(3*m);
144  S4 = S4<<(4*m);
145
146  ret = S0 + S1 + S2 + S3 + S4;
147  ret = sign *ret;
148
149  return ret;
150}
151
152define toomcook4(a,b)
153{
154
155  local  a0 a1 a2 a3 b0 b1 b2 b3 b4 ret tmp tmp2 tmp3 sign;
156  local  m alen blen mask;
157  local w1, w2, w3, w4, w5, w6, w7;
158
159  if(!isint(a) || !isint(b))
160    return newerror("toomcook4(a,b): a and/or b is not integer");
161
162  alen = digits(a,2);
163  blen = digits(b,2);
164
165  sign = sgn(a) * sgn(b);
166
167  if(sign == 0) return 0;
168
169  m = min(alen//4,blen//4);
170  mask = ~-(1<<m);
171
172  if(isdefined("test8900")){
173    if(m < 100) return toomcook3(a,b);
174  }
175  else{
176    if(m < 256*3072) return toomcook3(a,b);
177  }
178
179  a = abs(a);
180  b = abs(b);
181
182
183  a0 = a & mask;
184  a1 = (a>>m) & mask;
185  a2 = (a>>(2*m)) & mask;
186  a3 = (a>>(3*m));
187
188  b0 =  b         & mask;
189  b1 = (b>>m)     & mask;
190  b2 = (b>>(2*m)) & mask;
191  b3 = (b>>(3*m));
192
193  /*
194    Bodrato / Zanoni
195  */
196
197  w3 = a3 + (a1 + (a2 + a0));
198  w7 = b3 + (b1 + (b2 + b0));
199
200  w4 = -a3 + (-a1 + (a2 + a0));
201  w5 = -b3 + (-b1 + (b2 + b0));
202
203  w3 = toomcook4(w3, w7);
204  w4 = toomcook4(w4, w5);
205
206  w5 = a3 + ((a1<<2) + ((a2<<1) + (a0<<3)));
207  w2 = b3 + ((b1<<2) + ((b2<<1) + (b0<<3)));
208
209  w6 = -a3 + (-(a1<<2) + ((a2<<1) + (a0<<3)));
210  w7 = -b3 + (-(b1<<2) + ((b2<<1) + (b0<<3)));
211
212  w5 = toomcook4(w5, w2);
213  w6 = toomcook4(w6, w7);
214
215
216  w2 = (a3<<3) + ((a1<<1) + ((a2<<2) + a0));
217  w7 = (b3<<3) + ((b1<<1) + ((b2<<2) + b0));
218
219
220  w2 = toomcook4(w2, w7);
221
222  w1 = toomcook4(a3, b3);
223  w7 = toomcook4(a0, b0);
224
225  w2 = w2 + w5;
226  w6 = w5 - w6;
227  w4 = w3 - w4;
228  w5 = w5 - w1;
229  w5 -= w7 << 6;
230  w4 = w4>>1;
231  w3 = w3 - w4;
232  w5 = w5<<1;
233  w5 = w5 - w6;
234  w2 -= w3 * 65;
235  w3 = w3 - w7;
236  w3 = w3 - w1;
237  w2 += w3 * 45;
238  w5 -= w3<<3;
239  w5 = w5//24;
240  w6 =  w6 - w2;
241  w2 -= w4<<4;
242  w2 = w2//18;
243  w3 = w3 - w5;
244  w4 = w4 - w2;
245  w6 += w2 * 30;
246  w6 = w6//60;
247  w2 = w2 - w6;
248
249
250  ret = w7 + (w6<<m) + (w5<<(2*m)) + (w4<<(3*m))+ (w3<<(4*m))+
251        (w2<<(5*m))+ (w1<<(6*m));
252
253  ret = sign *ret;
254
255  return ret;
256}
257
258define toomcook4square(a){
259  local  a0 a1 a2 a3 ret S0 S1 S2 S3 S4 S5 S6 S7 tmp tmp2 tmp3;
260  local  sign m alen mask;
261  local T0 T1 T2 T3 T4 T5 T6 T7 T8;
262
263  if(!isint(a) )return newerror("toomcook3square(a): a is not integer");
264
265  alen = digits(a,2);
266
267  sign = sgn(a) * sgn(a);
268  /* sgn(x) returns 0 if x = 0 */
269  if(sign == 0) return 0;
270
271  m = (alen)//4;
272  mask = ~-( 1 << m );
273
274  /*
275    cut-off at about 2 mio. dec. digits
276    TODO: check!
277  */
278
279  if(isdefined("test8900")){
280    if(m < 100) return toomcook3square(a);
281  }
282  else{
283    if(m < 512*3072) return toomcook3square(a);
284  }
285
286  a = abs(a);
287
288  a0 = a          & mask;
289  a1 = (a>>m)     & mask;
290  a2 = (a>>(2*m)) & mask;
291  a3 = (a>>(3*m)) ;
292
293  /*
294    Bodrato / Zanoni
295  */
296
297  S1 = toomcook4square(a0);
298  S2 = toomcook4(a0<<1,a1);
299  S3 = toomcook4((a0 + a1 - a2 - a3 ) , (a0 - a1 - a2 + a3 ));
300  S4 = toomcook4square(a0 + a1 + a2 + a3 );
301  S5 = toomcook4( (a0 - a2 )<<1 , (a1 - a3 ));
302  S6 = toomcook4(a3<<1 , a2);
303  S7 = toomcook4square(a3);
304
305  T1 = S3 + S4;
306  T2 = (T1 + S5 )>>1;
307  T3 = S2 + S6;
308  T4 = T2 - T3;
309  T5 = T3 - S5;
310  T6 = T4 - S3;
311  T7 = T4 - S1;
312  T8 = T6 - S7;
313
314  ret =   (S7<<(6*m)) + (S6<<(5*m)) + (T7<<(4*m))
315        + (T5<<(3*m)) + (T8<<(2*m)) + (S2<<(1*m))  + S1;
316
317  ret = sign *ret;
318
319  return ret;
320}
321
322/*
323 * TODO: Implement the asymmetric variations
324 */
325
326/*
327 * produce_long_random_number(n) returns large pseudo-random numbers.
328 * Really large numbers, e.g.:
329 *     produce_long_random_number(16)
330 * is ca 4,128,561 bits (ca 1,242,821 dec. digits) large. Exact length is not
331 * pre-determinable because of the chaotic output of the function random().
332 */
333define __CZ__produce_long_random_number(n)
334{
335  local ret k;
336  ret = 1;
337  if(!isint(n) || n<1)
338    return newerror("__CZ__produce_long_random_number(n): "
339		    "n is not an integer >=1");
340  for(k=0;k<n;k++){
341    ret += random();
342    ret  = toomcook4square(ret);
343  }
344  return ret;
345}
346
347
348/*
349 * restore internal function from resource debugging
350 * report important interface functions
351 */
352config("resource_debug", resource_debug_level),;
353if (config("resource_debug") & 3) {
354    print "toomcook3(a,b)";
355    print "toomcook3square(a)";
356    print "toomcook4(a,b)";
357    print "toomcook4square(a)";
358}
359