1 /* Implementation of the multiplication algorithm for Toom-Cook 6.5-way. 2 3 Contributed to the GNU project by Marco Bodrato. 4 5 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY 6 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 7 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. 8 9 Copyright 2009, 2010, 2012 Free Software Foundation, Inc. 10 11 This file is part of the GNU MP Library. 12 13 The GNU MP Library is free software; you can redistribute it and/or modify 14 it under the terms of either: 15 16 * the GNU Lesser General Public License as published by the Free 17 Software Foundation; either version 3 of the License, or (at your 18 option) any later version. 19 20 or 21 22 * the GNU General Public License as published by the Free Software 23 Foundation; either version 2 of the License, or (at your option) any 24 later version. 25 26 or both in parallel, as here. 27 28 The GNU MP Library is distributed in the hope that it will be useful, but 29 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 30 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 31 for more details. 32 33 You should have received copies of the GNU General Public License and the 34 GNU Lesser General Public License along with the GNU MP Library. If not, 35 see https://www.gnu.org/licenses/. */ 36 37 38 #include "gmp-impl.h" 39 40 41 #if GMP_NUMB_BITS < 21 42 #error Not implemented. 43 #endif 44 45 #if TUNE_PROGRAM_BUILD 46 #define MAYBE_mul_basecase 1 47 #define MAYBE_mul_toom22 1 48 #define MAYBE_mul_toom33 1 49 #define MAYBE_mul_toom6h 1 50 #else 51 #define MAYBE_mul_basecase \ 52 (MUL_TOOM6H_THRESHOLD < 6 * MUL_TOOM22_THRESHOLD) 53 #define MAYBE_mul_toom22 \ 54 (MUL_TOOM6H_THRESHOLD < 6 * MUL_TOOM33_THRESHOLD) 55 #define MAYBE_mul_toom33 \ 56 (MUL_TOOM6H_THRESHOLD < 6 * MUL_TOOM44_THRESHOLD) 57 #define MAYBE_mul_toom6h \ 58 (MUL_FFT_THRESHOLD >= 6 * MUL_TOOM6H_THRESHOLD) 59 #endif 60 61 #define TOOM6H_MUL_N_REC(p, a, b, f, p2, a2, b2, 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 if (f) \ 67 mpn_mul_basecase (p2, a2, n, b2, n); \ 68 } else if (MAYBE_mul_toom22 \ 69 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD)) { \ 70 mpn_toom22_mul (p, a, n, b, n, ws); \ 71 if (f) \ 72 mpn_toom22_mul (p2, a2, n, b2, n, ws); \ 73 } else if (MAYBE_mul_toom33 \ 74 && BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD)) { \ 75 mpn_toom33_mul (p, a, n, b, n, ws); \ 76 if (f) \ 77 mpn_toom33_mul (p2, a2, n, b2, n, ws); \ 78 } else if (! MAYBE_mul_toom6h \ 79 || BELOW_THRESHOLD (n, MUL_TOOM6H_THRESHOLD)) { \ 80 mpn_toom44_mul (p, a, n, b, n, ws); \ 81 if (f) \ 82 mpn_toom44_mul (p2, a2, n, b2, n, ws); \ 83 } else { \ 84 mpn_toom6h_mul (p, a, n, b, n, ws); \ 85 if (f) \ 86 mpn_toom6h_mul (p2, a2, n, b2, n, ws); \ 87 } \ 88 } while (0) 89 90 #define TOOM6H_MUL_REC(p, a, na, b, nb, ws) \ 91 do { mpn_mul (p, a, na, b, nb); \ 92 } while (0) 93 94 /* Toom-6.5 , compute the product {pp,an+bn} <- {ap,an} * {bp,bn} 95 With: an >= bn >= 46, an*6 < bn * 17. 96 It _may_ work with bn<=46 and bn*17 < an*6 < bn*18 97 98 Evaluate in: infinity, +4, -4, +2, -2, +1, -1, +1/2, -1/2, +1/4, -1/4, 0. 99 */ 100 /* Estimate on needed scratch: 101 S(n) <= (n+5)\6*10+4+MAX(S((n+5)\6),1+2*(n+5)\6), 102 since n>42; S(n) <= ceil(log(n)/log(6))*(10+4)+n*12\6 < n*2 + lg2(n)*6 103 */ 104 105 void 106 mpn_toom6h_mul (mp_ptr pp, 107 mp_srcptr ap, mp_size_t an, 108 mp_srcptr bp, mp_size_t bn, mp_ptr scratch) 109 { 110 mp_size_t n, s, t; 111 int p, q, half; 112 int sign; 113 114 /***************************** decomposition *******************************/ 115 116 ASSERT (an >= bn); 117 /* Can not handle too much unbalancement */ 118 ASSERT (bn >= 42); 119 /* Can not handle too much unbalancement */ 120 ASSERT ((an*3 < bn * 8) || (bn >= 46 && an * 6 < bn * 17)); 121 122 /* Limit num/den is a rational number between 123 (12/11)^(log(4)/log(2*4-1)) and (12/11)^(log(6)/log(2*6-1)) */ 124 #define LIMIT_numerator (18) 125 #define LIMIT_denominat (17) 126 127 if (LIKELY (an * LIMIT_denominat < LIMIT_numerator * bn)) /* is 6*... < 6*... */ 128 { 129 n = 1 + (an - 1) / (size_t) 6; 130 p = q = 5; 131 half = 0; 132 133 s = an - 5 * n; 134 t = bn - 5 * n; 135 } 136 else { 137 if (an * 5 * LIMIT_numerator < LIMIT_denominat * 7 * bn) 138 { p = 7; q = 6; } 139 else if (an * 5 * LIMIT_denominat < LIMIT_numerator * 7 * bn) 140 { p = 7; q = 5; } 141 else if (an * LIMIT_numerator < LIMIT_denominat * 2 * bn) /* is 4*... < 8*... */ 142 { p = 8; q = 5; } 143 else if (an * LIMIT_denominat < LIMIT_numerator * 2 * bn) /* is 4*... < 8*... */ 144 { p = 8; q = 4; } 145 else 146 { p = 9; q = 4; } 147 148 half = (p ^ q) & 1; 149 n = 1 + (q * an >= p * bn ? (an - 1) / (size_t) p : (bn - 1) / (size_t) q); 150 p--; q--; 151 152 s = an - p * n; 153 t = bn - q * n; 154 155 /* With LIMIT = 16/15, the following recover is needed only if bn<=73*/ 156 if (half) { /* Recover from badly chosen splitting */ 157 if (UNLIKELY (s<1)) {p--; s+=n; half=0;} 158 else if (UNLIKELY (t<1)) {q--; t+=n; half=0;} 159 } 160 } 161 #undef LIMIT_numerator 162 #undef LIMIT_denominat 163 164 ASSERT (0 < s && s <= n); 165 ASSERT (0 < t && t <= n); 166 ASSERT (half || s + t > 3); 167 ASSERT (n > 2); 168 169 #define r4 (pp + 3 * n) /* 3n+1 */ 170 #define r2 (pp + 7 * n) /* 3n+1 */ 171 #define r0 (pp +11 * n) /* s+t <= 2*n */ 172 #define r5 (scratch) /* 3n+1 */ 173 #define r3 (scratch + 3 * n + 1) /* 3n+1 */ 174 #define r1 (scratch + 6 * n + 2) /* 3n+1 */ 175 #define v0 (pp + 7 * n) /* n+1 */ 176 #define v1 (pp + 8 * n+1) /* n+1 */ 177 #define v2 (pp + 9 * n+2) /* n+1 */ 178 #define v3 (scratch + 9 * n + 3) /* n+1 */ 179 #define wsi (scratch + 9 * n + 3) /* 3n+1 */ 180 #define wse (scratch +10 * n + 4) /* 2n+1 */ 181 182 /* Alloc also 3n+1 limbs for wsi... toom_interpolate_12pts may 183 need all of them */ 184 /* if (scratch == NULL) */ 185 /* scratch = TMP_SALLOC_LIMBS(mpn_toom6_sqr_itch(n * 6)); */ 186 ASSERT (12 * n + 6 <= mpn_toom6h_mul_itch(an,bn)); 187 ASSERT (12 * n + 6 <= mpn_toom6_sqr_itch(n * 6)); 188 189 /********************** evaluation and recursive calls *********************/ 190 /* $\pm1/2$ */ 191 sign = mpn_toom_eval_pm2rexp (v2, v0, p, ap, n, s, 1, pp) ^ 192 mpn_toom_eval_pm2rexp (v3, v1, q, bp, n, t, 1, pp); 193 /* A(-1/2)*B(-1/2)*2^. */ /* A(+1/2)*B(+1/2)*2^. */ 194 TOOM6H_MUL_N_REC(pp, v0, v1, 2, r5, v2, v3, n + 1, wse); 195 mpn_toom_couple_handling (r5, 2 * n + 1, pp, sign, n, 1+half , half); 196 197 /* $\pm1$ */ 198 sign = mpn_toom_eval_pm1 (v2, v0, p, ap, n, s, pp); 199 if (UNLIKELY (q == 3)) 200 sign ^= mpn_toom_eval_dgr3_pm1 (v3, v1, bp, n, t, pp); 201 else 202 sign ^= mpn_toom_eval_pm1 (v3, v1, q, bp, n, t, pp); 203 /* A(-1)*B(-1) */ /* A(1)*B(1) */ 204 TOOM6H_MUL_N_REC(pp, v0, v1, 2, r3, v2, v3, n + 1, wse); 205 mpn_toom_couple_handling (r3, 2 * n + 1, pp, sign, n, 0, 0); 206 207 /* $\pm4$ */ 208 sign = mpn_toom_eval_pm2exp (v2, v0, p, ap, n, s, 2, pp) ^ 209 mpn_toom_eval_pm2exp (v3, v1, q, bp, n, t, 2, pp); 210 /* A(-4)*B(-4) */ 211 TOOM6H_MUL_N_REC(pp, v0, v1, 2, r1, v2, v3, n + 1, wse); /* A(+4)*B(+4) */ 212 mpn_toom_couple_handling (r1, 2 * n + 1, pp, sign, n, 2, 4); 213 214 /* $\pm1/4$ */ 215 sign = mpn_toom_eval_pm2rexp (v2, v0, p, ap, n, s, 2, pp) ^ 216 mpn_toom_eval_pm2rexp (v3, v1, q, bp, n, t, 2, pp); 217 /* A(-1/4)*B(-1/4)*4^. */ /* A(+1/4)*B(+1/4)*4^. */ 218 TOOM6H_MUL_N_REC(pp, v0, v1, 2, r4, v2, v3, n + 1, wse); 219 mpn_toom_couple_handling (r4, 2 * n + 1, pp, sign, n, 2*(1+half), 2*(half)); 220 221 /* $\pm2$ */ 222 sign = mpn_toom_eval_pm2 (v2, v0, p, ap, n, s, pp) ^ 223 mpn_toom_eval_pm2 (v3, v1, q, bp, n, t, pp); 224 /* A(-2)*B(-2) */ /* A(+2)*B(+2) */ 225 TOOM6H_MUL_N_REC(pp, v0, v1, 2, r2, v2, v3, n + 1, wse); 226 mpn_toom_couple_handling (r2, 2 * n + 1, pp, sign, n, 1, 2); 227 228 #undef v0 229 #undef v1 230 #undef v2 231 #undef v3 232 #undef wse 233 234 /* A(0)*B(0) */ 235 TOOM6H_MUL_N_REC(pp, ap, bp, 0, pp, ap, bp, n, wsi); 236 237 /* Infinity */ 238 if (UNLIKELY (half != 0)) { 239 if (s > t) { 240 TOOM6H_MUL_REC(r0, ap + p * n, s, bp + q * n, t, wsi); 241 } else { 242 TOOM6H_MUL_REC(r0, bp + q * n, t, ap + p * n, s, wsi); 243 }; 244 }; 245 246 mpn_toom_interpolate_12pts (pp, r1, r3, r5, n, s+t, half, wsi); 247 248 #undef r0 249 #undef r1 250 #undef r2 251 #undef r3 252 #undef r4 253 #undef r5 254 #undef wsi 255 } 256 257 #undef TOOM6H_MUL_N_REC 258 #undef TOOM6H_MUL_REC 259 #undef MAYBE_mul_basecase 260 #undef MAYBE_mul_toom22 261 #undef MAYBE_mul_toom33 262 #undef MAYBE_mul_toom6h 263