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_SALLOC_LIMBS (n + 1); 120 asm1 = TMP_SALLOC_LIMBS (n + 1); 121 as2 = TMP_SALLOC_LIMBS (n + 1); 122 ash = TMP_SALLOC_LIMBS (n + 1); 123 asmh = TMP_SALLOC_LIMBS (n + 1); 124 125 bs1 = TMP_SALLOC_LIMBS (n + 1); 126 bsm1 = TMP_SALLOC_LIMBS (n + 1); 127 bs2 = TMP_SALLOC_LIMBS (n + 1); 128 bsh = TMP_SALLOC_LIMBS (n + 1); 129 bsmh = TMP_SALLOC_LIMBS (n + 1); 130 131 gp = pp; 132 hp = pp + n + 1; 133 134 flags = 0; 135 136 /* Compute as1 and asm1. */ 137 gp[n] = mpn_add_n (gp, a0, a2, n); 138 hp[n] = mpn_add (hp, a1, n, a3, s); 139 #if HAVE_NATIVE_mpn_addsub_n 140 if (mpn_cmp (gp, hp, n + 1) < 0) 141 { 142 mpn_addsub_n (as1, asm1, hp, gp, n + 1); 143 flags ^= toom4_w3_neg; 144 } 145 else 146 { 147 mpn_addsub_n (as1, asm1, gp, hp, n + 1); 148 } 149 #else 150 mpn_add_n (as1, gp, hp, n + 1); 151 if (mpn_cmp (gp, hp, n + 1) < 0) 152 { 153 mpn_sub_n (asm1, hp, gp, n + 1); 154 flags ^= toom4_w3_neg; 155 } 156 else 157 { 158 mpn_sub_n (asm1, gp, hp, n + 1); 159 } 160 #endif 161 162 /* Compute as2. */ 163 #if HAVE_NATIVE_mpn_addlsh1_n 164 cy = mpn_addlsh1_n (as2, a2, a3, s); 165 if (s != n) 166 cy = mpn_add_1 (as2 + s, a2 + s, n - s, cy); 167 cy = 2 * cy + mpn_addlsh1_n (as2, a1, as2, n); 168 cy = 2 * cy + mpn_addlsh1_n (as2, a0, as2, n); 169 #else 170 cy = mpn_lshift (as2, a3, s, 1); 171 cy += mpn_add_n (as2, a2, as2, s); 172 if (s != n) 173 cy = mpn_add_1 (as2 + s, a2 + s, n - s, cy); 174 cy = 2 * cy + mpn_lshift (as2, as2, n, 1); 175 cy += mpn_add_n (as2, a1, as2, n); 176 cy = 2 * cy + mpn_lshift (as2, as2, n, 1); 177 cy += mpn_add_n (as2, a0, as2, n); 178 #endif 179 as2[n] = cy; 180 181 /* Compute ash and asmh. */ 182 cy = mpn_lshift (gp, a0, n, 3); /* 8a0 */ 183 #if HAVE_NATIVE_mpn_addlsh1_n 184 gp[n] = cy + mpn_addlsh1_n (gp, gp, a2, n); /* 8a0 + 2a2 */ 185 #else 186 cy += mpn_lshift (hp, a2, n, 1); /* 2a2 */ 187 gp[n] = cy + mpn_add_n (gp, gp, hp, n); /* 8a0 + 2a2 */ 188 #endif 189 cy = mpn_lshift (hp, a1, n, 2); /* 4a1 */ 190 hp[n] = cy + mpn_add (hp, hp, n, a3, s); /* 4a1 + a3 */ 191 #if HAVE_NATIVE_mpn_addsub_n 192 if (mpn_cmp (gp, hp, n + 1) < 0) 193 { 194 mpn_addsub_n (ash, asmh, hp, gp, n + 1); 195 flags ^= toom4_w1_neg; 196 } 197 else 198 { 199 mpn_addsub_n (ash, asmh, gp, hp, n + 1); 200 } 201 #else 202 mpn_add_n (ash, gp, hp, n + 1); 203 if (mpn_cmp (gp, hp, n + 1) < 0) 204 { 205 mpn_sub_n (asmh, hp, gp, n + 1); 206 flags ^= toom4_w1_neg; 207 } 208 else 209 { 210 mpn_sub_n (asmh, gp, hp, n + 1); 211 } 212 #endif 213 214 /* Compute bs1 and bsm1. */ 215 gp[n] = mpn_add_n (gp, b0, b2, n); 216 hp[n] = mpn_add (hp, b1, n, b3, t); 217 #if HAVE_NATIVE_mpn_addsub_n 218 if (mpn_cmp (gp, hp, n + 1) < 0) 219 { 220 mpn_addsub_n (bs1, bsm1, hp, gp, n + 1); 221 flags ^= toom4_w3_neg; 222 } 223 else 224 { 225 mpn_addsub_n (bs1, bsm1, gp, hp, n + 1); 226 } 227 #else 228 mpn_add_n (bs1, gp, hp, n + 1); 229 if (mpn_cmp (gp, hp, n + 1) < 0) 230 { 231 mpn_sub_n (bsm1, hp, gp, n + 1); 232 flags ^= toom4_w3_neg; 233 } 234 else 235 { 236 mpn_sub_n (bsm1, gp, hp, n + 1); 237 } 238 #endif 239 240 /* Compute bs2. */ 241 #if HAVE_NATIVE_mpn_addlsh1_n 242 cy = mpn_addlsh1_n (bs2, b2, b3, t); 243 if (t != n) 244 cy = mpn_add_1 (bs2 + t, b2 + t, n - t, cy); 245 cy = 2 * cy + mpn_addlsh1_n (bs2, b1, bs2, n); 246 cy = 2 * cy + mpn_addlsh1_n (bs2, b0, bs2, n); 247 #else 248 cy = mpn_lshift (bs2, b3, t, 1); 249 cy += mpn_add_n (bs2, b2, bs2, t); 250 if (t != n) 251 cy = mpn_add_1 (bs2 + t, b2 + t, n - t, cy); 252 cy = 2 * cy + mpn_lshift (bs2, bs2, n, 1); 253 cy += mpn_add_n (bs2, b1, bs2, n); 254 cy = 2 * cy + mpn_lshift (bs2, bs2, n, 1); 255 cy += mpn_add_n (bs2, b0, bs2, n); 256 #endif 257 bs2[n] = cy; 258 259 /* Compute bsh and bsmh. */ 260 cy = mpn_lshift (gp, b0, n, 3); /* 8b0 */ 261 #if HAVE_NATIVE_mpn_addlsh1_n 262 gp[n] = cy + mpn_addlsh1_n (gp, gp, b2, n); /* 8b0 + 2b2 */ 263 #else 264 cy += mpn_lshift (hp, b2, n, 1); /* 2b2 */ 265 gp[n] = cy + mpn_add_n (gp, gp, hp, n); /* 8b0 + 2b2 */ 266 #endif 267 cy = mpn_lshift (hp, b1, n, 2); /* 4b1 */ 268 hp[n] = cy + mpn_add (hp, hp, n, b3, t); /* 4b1 + b3 */ 269 #if HAVE_NATIVE_mpn_addsub_n 270 if (mpn_cmp (gp, hp, n + 1) < 0) 271 { 272 mpn_addsub_n (bsh, bsmh, hp, gp, n + 1); 273 flags ^= toom4_w1_neg; 274 } 275 else 276 { 277 mpn_addsub_n (bsh, bsmh, gp, hp, n + 1); 278 } 279 #else 280 mpn_add_n (bsh, gp, hp, n + 1); 281 if (mpn_cmp (gp, hp, n + 1) < 0) 282 { 283 mpn_sub_n (bsmh, hp, gp, n + 1); 284 flags ^= toom4_w1_neg; 285 } 286 else 287 { 288 mpn_sub_n (bsmh, gp, hp, n + 1); 289 } 290 #endif 291 292 ASSERT (as1[n] <= 3); 293 ASSERT (bs1[n] <= 3); 294 ASSERT (asm1[n] <= 1); 295 ASSERT (bsm1[n] <= 1); 296 ASSERT (as2[n] <= 14); 297 ASSERT (bs2[n] <= 14); 298 ASSERT (ash[n] <= 14); 299 ASSERT (bsh[n] <= 14); 300 ASSERT (asmh[n] <= 9); 301 ASSERT (bsmh[n] <= 9); 302 303 #define v0 pp /* 2n */ 304 #define v1 (scratch + 6 * n + 6) /* 2n+1 */ 305 #define vm1 scratch /* 2n+1 */ 306 #define v2 (scratch + 2 * n + 2) /* 2n+1 */ 307 #define vinf (pp + 6 * n) /* s+t */ 308 #define vh (pp + 2 * n) /* 2n+1 */ 309 #define vmh (scratch + 4 * n + 4) 310 #define scratch_out (scratch + 8 * n + 8) 311 312 /* vm1, 2n+1 limbs */ 313 TOOM44_MUL_N_REC (vm1, asm1, bsm1, n + 1, scratch_out); /* vm1, 2n+1 limbs */ 314 315 TOOM44_MUL_N_REC (v2 , as2 , bs2 , n + 1, scratch_out); /* v2, 2n+1 limbs */ 316 317 if (s > t) mpn_mul (vinf, a3, s, b3, t); 318 else TOOM44_MUL_N_REC (vinf, a3, b3, s, scratch_out); /* vinf, s+t limbs */ 319 320 TOOM44_MUL_N_REC (v1 , as1 , bs1 , n + 1, scratch_out); /* v1, 2n+1 limbs */ 321 322 TOOM44_MUL_N_REC (vh , ash , bsh , n + 1, scratch_out); 323 324 TOOM44_MUL_N_REC (vmh, asmh, bsmh, n + 1, scratch_out); 325 326 TOOM44_MUL_N_REC (v0 , ap , bp , n , scratch_out); /* v0, 2n limbs */ 327 328 mpn_toom_interpolate_7pts (pp, n, flags, vmh, vm1, v1, v2, s + t, scratch_out); 329 330 TMP_FREE; 331 } 332