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