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-2008, 2013 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 either: 16 17 * the GNU Lesser General Public License as published by the Free 18 Software Foundation; either version 3 of the License, or (at your 19 option) any later version. 20 21 or 22 23 * the GNU General Public License as published by the Free Software 24 Foundation; either version 2 of the License, or (at your option) any 25 later version. 26 27 or both in parallel, as here. 28 29 The GNU MP Library is distributed in the hope that it will be useful, but 30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 31 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 32 for more details. 33 34 You should have received copies of the GNU General Public License and the 35 GNU Lesser General Public License along with the GNU MP Library. If not, 36 see https://www.gnu.org/licenses/. */ 37 38 39 #include "gmp-impl.h" 40 41 /* Evaluate in: 0, +1, -1, +2, -2, 1/2, +inf 42 43 <-s--><--n--><--n--><--n--> 44 ____ ______ ______ ______ 45 |_a3_|___a2_|___a1_|___a0_| 46 |b3_|___b2_|___b1_|___b0_| 47 <-t-><--n--><--n--><--n--> 48 49 v0 = a0 * b0 # A(0)*B(0) 50 v1 = ( a0+ a1+ a2+ a3)*( b0+ b1+ b2+ b3) # A(1)*B(1) ah <= 3 bh <= 3 51 vm1 = ( a0- a1+ a2- a3)*( b0- b1+ b2- b3) # A(-1)*B(-1) |ah| <= 1 |bh| <= 1 52 v2 = ( a0+2a1+4a2+8a3)*( b0+2b1+4b2+8b3) # A(2)*B(2) ah <= 14 bh <= 14 53 vm2 = ( a0-2a1+4a2-8a3)*( b0-2b1+4b2-8b3) # A(2)*B(2) ah <= 9 |bh| <= 9 54 vh = (8a0+4a1+2a2+ a3)*(8b0+4b1+2b2+ b3) # A(1/2)*B(1/2) ah <= 14 bh <= 14 55 vinf= a3 * b2 # A(inf)*B(inf) 56 */ 57 58 #if TUNE_PROGRAM_BUILD 59 #define MAYBE_mul_basecase 1 60 #define MAYBE_mul_toom22 1 61 #define MAYBE_mul_toom44 1 62 #else 63 #define MAYBE_mul_basecase \ 64 (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM22_THRESHOLD) 65 #define MAYBE_mul_toom22 \ 66 (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM33_THRESHOLD) 67 #define MAYBE_mul_toom44 \ 68 (MUL_TOOM6H_THRESHOLD >= 4 * MUL_TOOM44_THRESHOLD) 69 #endif 70 71 #define TOOM44_MUL_N_REC(p, a, b, n, ws) \ 72 do { \ 73 if (MAYBE_mul_basecase \ 74 && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD)) \ 75 mpn_mul_basecase (p, a, n, b, n); \ 76 else if (MAYBE_mul_toom22 \ 77 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD)) \ 78 mpn_toom22_mul (p, a, n, b, n, ws); \ 79 else if (! MAYBE_mul_toom44 \ 80 || BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD)) \ 81 mpn_toom33_mul (p, a, n, b, n, ws); \ 82 else \ 83 mpn_toom44_mul (p, a, n, b, n, ws); \ 84 } while (0) 85 86 /* Use of scratch space. In the product area, we store 87 88 ___________________ 89 |vinf|____|_v1_|_v0_| 90 s+t 2n-1 2n+1 2n 91 92 The other recursive products, vm1, v2, vm2, vh are stored in the 93 scratch area. When computing them, we use the product area for 94 intermediate values. 95 96 Next, we compute v1. We can store the intermediate factors at v0 97 and at vh + 2n + 2. 98 99 Finally, for v0 and vinf, factors are parts of the input operands, 100 and we need scratch space only for the recursive multiplication. 101 102 In all, if S(an) is the scratch need, the needed space is bounded by 103 104 S(an) <= 4 (2*ceil(an/4) + 1) + 1 + S(ceil(an/4) + 1) 105 106 which should give S(n) = 8 n/3 + c log(n) for some constant c. 107 */ 108 109 void 110 mpn_toom44_mul (mp_ptr pp, 111 mp_srcptr ap, mp_size_t an, 112 mp_srcptr bp, mp_size_t bn, 113 mp_ptr scratch) 114 { 115 mp_size_t n, s, t; 116 mp_limb_t cy; 117 enum toom7_flags flags; 118 119 #define a0 ap 120 #define a1 (ap + n) 121 #define a2 (ap + 2*n) 122 #define a3 (ap + 3*n) 123 #define b0 bp 124 #define b1 (bp + n) 125 #define b2 (bp + 2*n) 126 #define b3 (bp + 3*n) 127 128 ASSERT (an >= bn); 129 130 n = (an + 3) >> 2; 131 132 s = an - 3 * n; 133 t = bn - 3 * n; 134 135 ASSERT (0 < s && s <= n); 136 ASSERT (0 < t && t <= n); 137 ASSERT (s >= t); 138 139 /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the 140 * following limb, so these must be computed in order, and we need a 141 * one limb gap to tp. */ 142 #define v0 pp /* 2n */ 143 #define v1 (pp + 2 * n) /* 2n+1 */ 144 #define vinf (pp + 6 * n) /* s+t */ 145 #define v2 scratch /* 2n+1 */ 146 #define vm2 (scratch + 2 * n + 1) /* 2n+1 */ 147 #define vh (scratch + 4 * n + 2) /* 2n+1 */ 148 #define vm1 (scratch + 6 * n + 3) /* 2n+1 */ 149 #define tp (scratch + 8*n + 5) 150 151 /* apx and bpx must not overlap with v1 */ 152 #define apx pp /* n+1 */ 153 #define amx (pp + n + 1) /* n+1 */ 154 #define bmx (pp + 2*n + 2) /* n+1 */ 155 #define bpx (pp + 4*n + 2) /* n+1 */ 156 157 /* Total scratch need: 8*n + 5 + scratch for recursive calls. This 158 gives roughly 32 n/3 + log term. */ 159 160 /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3. */ 161 flags = (enum toom7_flags) (toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (apx, amx, ap, n, s, tp)); 162 163 /* Compute bpx = b0 + 2 b1 + 4 b2 + 8 b3 and bmx = b0 - 2 b1 + 4 b2 - 8 b3. */ 164 flags = (enum toom7_flags) (flags ^ (toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (bpx, bmx, bp, n, t, tp))); 165 166 TOOM44_MUL_N_REC (v2, apx, bpx, n + 1, tp); /* v2, 2n+1 limbs */ 167 TOOM44_MUL_N_REC (vm2, amx, bmx, n + 1, tp); /* vm2, 2n+1 limbs */ 168 169 /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */ 170 #if HAVE_NATIVE_mpn_addlsh1_n 171 cy = mpn_addlsh1_n (apx, a1, a0, n); 172 cy = 2*cy + mpn_addlsh1_n (apx, a2, apx, n); 173 if (s < n) 174 { 175 mp_limb_t cy2; 176 cy2 = mpn_addlsh1_n (apx, a3, apx, s); 177 apx[n] = 2*cy + mpn_lshift (apx + s, apx + s, n - s, 1); 178 MPN_INCR_U (apx + s, n+1-s, cy2); 179 } 180 else 181 apx[n] = 2*cy + mpn_addlsh1_n (apx, a3, apx, n); 182 #else 183 cy = mpn_lshift (apx, a0, n, 1); 184 cy += mpn_add_n (apx, apx, a1, n); 185 cy = 2*cy + mpn_lshift (apx, apx, n, 1); 186 cy += mpn_add_n (apx, apx, a2, n); 187 cy = 2*cy + mpn_lshift (apx, apx, n, 1); 188 apx[n] = cy + mpn_add (apx, apx, n, a3, s); 189 #endif 190 191 /* Compute bpx = 8 b0 + 4 b1 + 2 b2 + b3 = (((2*b0 + b1) * 2 + b2) * 2 + b3 */ 192 #if HAVE_NATIVE_mpn_addlsh1_n 193 cy = mpn_addlsh1_n (bpx, b1, b0, n); 194 cy = 2*cy + mpn_addlsh1_n (bpx, b2, bpx, n); 195 if (t < n) 196 { 197 mp_limb_t cy2; 198 cy2 = mpn_addlsh1_n (bpx, b3, bpx, t); 199 bpx[n] = 2*cy + mpn_lshift (bpx + t, bpx + t, n - t, 1); 200 MPN_INCR_U (bpx + t, n+1-t, cy2); 201 } 202 else 203 bpx[n] = 2*cy + mpn_addlsh1_n (bpx, b3, bpx, n); 204 #else 205 cy = mpn_lshift (bpx, b0, n, 1); 206 cy += mpn_add_n (bpx, bpx, b1, n); 207 cy = 2*cy + mpn_lshift (bpx, bpx, n, 1); 208 cy += mpn_add_n (bpx, bpx, b2, n); 209 cy = 2*cy + mpn_lshift (bpx, bpx, n, 1); 210 bpx[n] = cy + mpn_add (bpx, bpx, n, b3, t); 211 #endif 212 213 ASSERT (apx[n] < 15); 214 ASSERT (bpx[n] < 15); 215 216 TOOM44_MUL_N_REC (vh, apx, bpx, n + 1, tp); /* vh, 2n+1 limbs */ 217 218 /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3. */ 219 flags = (enum toom7_flags) (flags | (toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (apx, amx, ap, n, s, tp))); 220 221 /* Compute bpx = b0 + b1 + b2 + b3 and bmx = b0 - b1 + b2 - b3. */ 222 flags = (enum toom7_flags) (flags ^ (toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (bpx, bmx, bp, n, t, tp))); 223 224 TOOM44_MUL_N_REC (vm1, amx, bmx, n + 1, tp); /* vm1, 2n+1 limbs */ 225 /* Clobbers amx, bmx. */ 226 TOOM44_MUL_N_REC (v1, apx, bpx, n + 1, tp); /* v1, 2n+1 limbs */ 227 228 TOOM44_MUL_N_REC (v0, a0, b0, n, tp); 229 if (s > t) 230 mpn_mul (vinf, a3, s, b3, t); 231 else 232 TOOM44_MUL_N_REC (vinf, a3, b3, s, tp); /* vinf, s+t limbs */ 233 234 mpn_toom_interpolate_7pts (pp, n, flags, vm2, vm1, v2, vh, s + t, tp); 235 } 236