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