1 /* mpn_toom44_mul -- Multiply {ap,an} and {bp,bn} where an and bn are close in 2 size. Or more accurately, bn <= an < (4/3)bn. 3 4 Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato. 5 6 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY 7 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 8 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. 9 10 Copyright 2006, 2007, 2008 Free Software Foundation, Inc. 11 12 This file is part of the GNU MP Library. 13 14 The GNU MP Library is free software; you can redistribute it and/or modify 15 it under the terms of the GNU Lesser General Public License as published by 16 the Free Software Foundation; either version 3 of the License, or (at your 17 option) any later version. 18 19 The GNU MP Library is distributed in the hope that it will be useful, but 20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 21 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 22 License for more details. 23 24 You should have received a copy of the GNU Lesser General Public License 25 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 26 27 28 /* 29 Things to work on: 30 31 1. Trim allocation. The allocations for as1, asm1, bs1, and bsm1 could be 32 avoided by instead reusing the pp area and the scratch area. 33 2. Use new toom functions for the recursive calls. 34 */ 35 36 #include "gmp.h" 37 #include "gmp-impl.h" 38 39 /* Evaluate in: -1, -1/2, 0, +1/2, +1, +2, +inf 40 41 <-s--><--n--><--n--><--n--> 42 ____ ______ ______ ______ 43 |_a3_|___a2_|___a1_|___a0_| 44 |b3_|___b2_|___b1_|___b0_| 45 <-t-><--n--><--n--><--n--> 46 47 v0 = a0 * b0 # A(0)*B(0) 48 v1 = ( a0+ a1+ a2+ a3)*( b0+ b1+ b2+ b3) # A(1)*B(1) ah <= 3 bh <= 3 49 vm1 = ( a0- a1+ a2- a3)*( b0- b1+ b2- b3) # A(-1)*B(-1) |ah| <= 1 |bh| <= 1 50 v2 = ( a0+2a1+4a2+8a3)*( b0+2b1+4b2+8b3) # A(2)*B(2) ah <= 14 bh <= 14 51 vh = (8a0+4a1+2a2+ a3)*(8b0+4b1+2b2+ b3) # A(1/2)*B(1/2) ah <= 14 bh <= 14 52 vmh = (8a0-4a1+2a2- a3)*(8b0-4b1+2b2- b3) # A(-1/2)*B(-1/2) -4<=ah<=9 -4<=bh<=9 53 vinf= a3 * b2 # A(inf)*B(inf) 54 */ 55 56 #if TUNE_PROGRAM_BUILD 57 #define MAYBE_mul_basecase 1 58 #define MAYBE_mul_toom22 1 59 #define MAYBE_mul_toom44 1 60 #else 61 #define MAYBE_mul_basecase \ 62 (MUL_TOOM44_THRESHOLD < 4 * MUL_KARATSUBA_THRESHOLD) 63 #define MAYBE_mul_toom22 \ 64 (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM33_THRESHOLD) 65 #define MAYBE_mul_toom44 \ 66 (MUL_FFT_THRESHOLD >= 4 * MUL_TOOM44_THRESHOLD) 67 #endif 68 69 #define TOOM44_MUL_N_REC(p, a, b, n, ws) \ 70 do { \ 71 if (MAYBE_mul_basecase \ 72 && BELOW_THRESHOLD (n, MUL_KARATSUBA_THRESHOLD)) \ 73 mpn_mul_basecase (p, a, n, b, n); \ 74 else if (MAYBE_mul_toom22 \ 75 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD)) \ 76 mpn_kara_mul_n (p, a, b, n, ws); \ 77 else if (! MAYBE_mul_toom44 \ 78 || BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD)) \ 79 mpn_toom3_mul_n (p, a, b, n, ws); \ 80 else \ 81 mpn_toom44_mul (p, a, n, b, n, ws); \ 82 } while (0) 83 84 void 85 mpn_toom44_mul (mp_ptr pp, 86 mp_srcptr ap, mp_size_t an, 87 mp_srcptr bp, mp_size_t bn, 88 mp_ptr scratch) 89 { 90 mp_size_t n, s, t; 91 mp_limb_t cy; 92 mp_ptr gp, hp; 93 mp_ptr as1, asm1, as2, ash, asmh; 94 mp_ptr bs1, bsm1, bs2, bsh, bsmh; 95 enum toom4_flags flags; 96 TMP_DECL; 97 98 #define a0 ap 99 #define a1 (ap + n) 100 #define a2 (ap + 2*n) 101 #define a3 (ap + 3*n) 102 #define b0 bp 103 #define b1 (bp + n) 104 #define b2 (bp + 2*n) 105 #define b3 (bp + 3*n) 106 107 n = (an + 3) >> 2; 108 109 s = an - 3 * n; 110 t = bn - 3 * n; 111 112 ASSERT (an >= bn); 113 114 ASSERT (0 < s && s <= n); 115 ASSERT (0 < t && t <= n); 116 117 TMP_MARK; 118 119 as1 = TMP_ALLOC_LIMBS (10 * n + 10); 120 asm1 = as1 + n + 1; 121 as2 = asm1 + n + 1; 122 ash = as2 + n + 1; 123 asmh = ash + n + 1; 124 bs1 = asmh + n + 1; 125 bsm1 = bs1 + n + 1; 126 bs2 = bsm1 + n + 1; 127 bsh = bs2 + n + 1; 128 bsmh = bsh + n + 1; 129 130 gp = pp; 131 hp = pp + n + 1; 132 133 flags = 0; 134 135 /* Compute as1 and asm1. */ 136 gp[n] = mpn_add_n (gp, a0, a2, n); 137 hp[n] = mpn_add (hp, a1, n, a3, s); 138 #if HAVE_NATIVE_mpn_addsub_n 139 if (mpn_cmp (gp, hp, n + 1) < 0) 140 { 141 mpn_addsub_n (as1, asm1, hp, gp, n + 1); 142 flags ^= toom4_w3_neg; 143 } 144 else 145 { 146 mpn_addsub_n (as1, asm1, gp, hp, n + 1); 147 } 148 #else 149 mpn_add_n (as1, gp, hp, n + 1); 150 if (mpn_cmp (gp, hp, n + 1) < 0) 151 { 152 mpn_sub_n (asm1, hp, gp, n + 1); 153 flags ^= toom4_w3_neg; 154 } 155 else 156 { 157 mpn_sub_n (asm1, gp, hp, n + 1); 158 } 159 #endif 160 161 /* Compute as2. */ 162 #if HAVE_NATIVE_mpn_addlsh1_n 163 cy = mpn_addlsh1_n (as2, a2, a3, s); 164 if (s != n) 165 cy = mpn_add_1 (as2 + s, a2 + s, n - s, cy); 166 cy = 2 * cy + mpn_addlsh1_n (as2, a1, as2, n); 167 cy = 2 * cy + mpn_addlsh1_n (as2, a0, as2, n); 168 #else 169 cy = mpn_lshift (as2, a3, s, 1); 170 cy += mpn_add_n (as2, a2, as2, s); 171 if (s != n) 172 cy = mpn_add_1 (as2 + s, a2 + s, n - s, cy); 173 cy = 2 * cy + mpn_lshift (as2, as2, n, 1); 174 cy += mpn_add_n (as2, a1, as2, n); 175 cy = 2 * cy + mpn_lshift (as2, as2, n, 1); 176 cy += mpn_add_n (as2, a0, as2, n); 177 #endif 178 as2[n] = cy; 179 180 /* Compute ash and asmh. */ 181 cy = mpn_lshift (gp, a0, n, 3); /* 8a0 */ 182 #if HAVE_NATIVE_mpn_addlsh1_n 183 gp[n] = cy + mpn_addlsh1_n (gp, gp, a2, n); /* 8a0 + 2a2 */ 184 #else 185 cy += mpn_lshift (hp, a2, n, 1); /* 2a2 */ 186 gp[n] = cy + mpn_add_n (gp, gp, hp, n); /* 8a0 + 2a2 */ 187 #endif 188 cy = mpn_lshift (hp, a1, n, 2); /* 4a1 */ 189 hp[n] = cy + mpn_add (hp, hp, n, a3, s); /* 4a1 + a3 */ 190 #if HAVE_NATIVE_mpn_addsub_n 191 if (mpn_cmp (gp, hp, n + 1) < 0) 192 { 193 mpn_addsub_n (ash, asmh, hp, gp, n + 1); 194 flags ^= toom4_w1_neg; 195 } 196 else 197 { 198 mpn_addsub_n (ash, asmh, gp, hp, n + 1); 199 } 200 #else 201 mpn_add_n (ash, gp, hp, n + 1); 202 if (mpn_cmp (gp, hp, n + 1) < 0) 203 { 204 mpn_sub_n (asmh, hp, gp, n + 1); 205 flags ^= toom4_w1_neg; 206 } 207 else 208 { 209 mpn_sub_n (asmh, gp, hp, n + 1); 210 } 211 #endif 212 213 /* Compute bs1 and bsm1. */ 214 gp[n] = mpn_add_n (gp, b0, b2, n); 215 hp[n] = mpn_add (hp, b1, n, b3, t); 216 #if HAVE_NATIVE_mpn_addsub_n 217 if (mpn_cmp (gp, hp, n + 1) < 0) 218 { 219 mpn_addsub_n (bs1, bsm1, hp, gp, n + 1); 220 flags ^= toom4_w3_neg; 221 } 222 else 223 { 224 mpn_addsub_n (bs1, bsm1, gp, hp, n + 1); 225 } 226 #else 227 mpn_add_n (bs1, gp, hp, n + 1); 228 if (mpn_cmp (gp, hp, n + 1) < 0) 229 { 230 mpn_sub_n (bsm1, hp, gp, n + 1); 231 flags ^= toom4_w3_neg; 232 } 233 else 234 { 235 mpn_sub_n (bsm1, gp, hp, n + 1); 236 } 237 #endif 238 239 /* Compute bs2. */ 240 #if HAVE_NATIVE_mpn_addlsh1_n 241 cy = mpn_addlsh1_n (bs2, b2, b3, t); 242 if (t != n) 243 cy = mpn_add_1 (bs2 + t, b2 + t, n - t, cy); 244 cy = 2 * cy + mpn_addlsh1_n (bs2, b1, bs2, n); 245 cy = 2 * cy + mpn_addlsh1_n (bs2, b0, bs2, n); 246 #else 247 cy = mpn_lshift (bs2, b3, t, 1); 248 cy += mpn_add_n (bs2, b2, bs2, t); 249 if (t != n) 250 cy = mpn_add_1 (bs2 + t, b2 + t, n - t, cy); 251 cy = 2 * cy + mpn_lshift (bs2, bs2, n, 1); 252 cy += mpn_add_n (bs2, b1, bs2, n); 253 cy = 2 * cy + mpn_lshift (bs2, bs2, n, 1); 254 cy += mpn_add_n (bs2, b0, bs2, n); 255 #endif 256 bs2[n] = cy; 257 258 /* Compute bsh and bsmh. */ 259 cy = mpn_lshift (gp, b0, n, 3); /* 8b0 */ 260 #if HAVE_NATIVE_mpn_addlsh1_n 261 gp[n] = cy + mpn_addlsh1_n (gp, gp, b2, n); /* 8b0 + 2b2 */ 262 #else 263 cy += mpn_lshift (hp, b2, n, 1); /* 2b2 */ 264 gp[n] = cy + mpn_add_n (gp, gp, hp, n); /* 8b0 + 2b2 */ 265 #endif 266 cy = mpn_lshift (hp, b1, n, 2); /* 4b1 */ 267 hp[n] = cy + mpn_add (hp, hp, n, b3, t); /* 4b1 + b3 */ 268 #if HAVE_NATIVE_mpn_addsub_n 269 if (mpn_cmp (gp, hp, n + 1) < 0) 270 { 271 mpn_addsub_n (bsh, bsmh, hp, gp, n + 1); 272 flags ^= toom4_w1_neg; 273 } 274 else 275 { 276 mpn_addsub_n (bsh, bsmh, gp, hp, n + 1); 277 } 278 #else 279 mpn_add_n (bsh, gp, hp, n + 1); 280 if (mpn_cmp (gp, hp, n + 1) < 0) 281 { 282 mpn_sub_n (bsmh, hp, gp, n + 1); 283 flags ^= toom4_w1_neg; 284 } 285 else 286 { 287 mpn_sub_n (bsmh, gp, hp, n + 1); 288 } 289 #endif 290 291 ASSERT (as1[n] <= 3); 292 ASSERT (bs1[n] <= 3); 293 ASSERT (asm1[n] <= 1); 294 ASSERT (bsm1[n] <= 1); 295 ASSERT (as2[n] <= 14); 296 ASSERT (bs2[n] <= 14); 297 ASSERT (ash[n] <= 14); 298 ASSERT (bsh[n] <= 14); 299 ASSERT (asmh[n] <= 9); 300 ASSERT (bsmh[n] <= 9); 301 302 #define v0 pp /* 2n */ 303 #define v1 (scratch + 6 * n + 6) /* 2n+1 */ 304 #define vm1 scratch /* 2n+1 */ 305 #define v2 (scratch + 2 * n + 2) /* 2n+1 */ 306 #define vinf (pp + 6 * n) /* s+t */ 307 #define vh (pp + 2 * n) /* 2n+1 */ 308 #define vmh (scratch + 4 * n + 4) 309 #define scratch_out (scratch + 8 * n + 8) 310 311 /* vm1, 2n+1 limbs */ 312 TOOM44_MUL_N_REC (vm1, asm1, bsm1, n + 1, scratch_out); /* vm1, 2n+1 limbs */ 313 314 TOOM44_MUL_N_REC (v2 , as2 , bs2 , n + 1, scratch_out); /* v2, 2n+1 limbs */ 315 316 if (s > t) mpn_mul (vinf, a3, s, b3, t); 317 else TOOM44_MUL_N_REC (vinf, a3, b3, s, scratch_out); /* vinf, s+t limbs */ 318 319 TOOM44_MUL_N_REC (v1 , as1 , bs1 , n + 1, scratch_out); /* v1, 2n+1 limbs */ 320 321 TOOM44_MUL_N_REC (vh , ash , bsh , n + 1, scratch_out); 322 323 TOOM44_MUL_N_REC (vmh, asmh, bsmh, n + 1, scratch_out); 324 325 TOOM44_MUL_N_REC (v0 , ap , bp , n , scratch_out); /* v0, 2n limbs */ 326 327 mpn_toom_interpolate_7pts (pp, n, flags, vmh, vm1, v1, v2, s + t, scratch_out); 328 329 TMP_FREE; 330 } 331