1 /* mpn_toom_interpolate_7pts -- Interpolate for toom44, 53, 62. 2 3 Contributed to the GNU project by Niels M�ller. 4 Improvements by 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, 2009 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 #include "gmp.h" 28 #include "gmp-impl.h" 29 30 #define BINVERT_3 MODLIMB_INVERSE_3 31 32 #define BINVERT_9 \ 33 ((((GMP_NUMB_MAX / 9) << (6 - GMP_NUMB_BITS % 6)) * 8 & GMP_NUMB_MAX) | 0x39) 34 35 #define BINVERT_15 \ 36 ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 4)) / 15) * 14 * 16 & GMP_NUMB_MAX) + 15)) 37 38 /* For the various mpn_divexact_byN here, fall back to using either 39 mpn_pi1_bdiv_q_1 or mpn_divexact_1. The former has less overhead and is 40 many faster if it is native. For now, since mpn_divexact_1 is native on 41 several platforms where mpn_pi1_bdiv_q_1 does not yet exist, do not use 42 mpn_pi1_bdiv_q_1 unconditionally. FIXME. */ 43 44 /* For odd divisors, mpn_divexact_1 works fine with two's complement. */ 45 #ifndef mpn_divexact_by3 46 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 47 #define mpn_divexact_by3(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,3,BINVERT_3,0) 48 #else 49 #define mpn_divexact_by3(dst,src,size) mpn_divexact_1(dst,src,size,3) 50 #endif 51 #endif 52 53 #ifndef mpn_divexact_by9 54 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 55 #define mpn_divexact_by9(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,9,BINVERT_9,0) 56 #else 57 #define mpn_divexact_by9(dst,src,size) mpn_divexact_1(dst,src,size,9) 58 #endif 59 #endif 60 61 #ifndef mpn_divexact_by15 62 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 63 #define mpn_divexact_by15(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,15,BINVERT_15,0) 64 #else 65 #define mpn_divexact_by15(dst,src,size) mpn_divexact_1(dst,src,size,15) 66 #endif 67 #endif 68 69 /* Interpolation for toom4, using the evaluation points 0, infinity, 70 1, -1, 2, -2, 1/2. More precisely, we want to compute 71 f(2^(GMP_NUMB_BITS * n)) for a polynomial f of degree 6, given the 72 seven values 73 74 w0 = f(0), 75 w1 = f(-2), 76 w2 = f(1), 77 w3 = f(-1), 78 w4 = f(2) 79 w5 = 64 * f(1/2) 80 w6 = limit at infinity of f(x) / x^6, 81 82 The result is 6*n + w6n limbs. At entry, w0 is stored at {rp, 2n }, 83 w2 is stored at { rp + 2n, 2n+1 }, and w6 is stored at { rp + 6n, 84 w6n }. The other values are 2n + 1 limbs each (with most 85 significant limbs small). f(-1) and f(-1/2) may be negative, signs 86 determined by the flag bits. Inputs are destroyed. 87 88 Needs (2*n + 1) limbs of temporary storage. 89 */ 90 91 void 92 mpn_toom_interpolate_7pts (mp_ptr rp, mp_size_t n, enum toom7_flags flags, 93 mp_ptr w1, mp_ptr w3, mp_ptr w4, mp_ptr w5, 94 mp_size_t w6n, mp_ptr tp) 95 { 96 mp_size_t m; 97 mp_limb_t cy; 98 99 m = 2*n + 1; 100 #define w0 rp 101 #define w2 (rp + 2*n) 102 #define w6 (rp + 6*n) 103 104 ASSERT (w6n > 0); 105 ASSERT (w6n <= 2*n); 106 107 /* Using formulas similar to Marco Bodrato's 108 109 W5 = W5 + W4 110 W1 =(W4 - W1)/2 111 W4 = W4 - W0 112 W4 =(W4 - W1)/4 - W6*16 113 W3 =(W2 - W3)/2 114 W2 = W2 - W3 115 116 W5 = W5 - W2*65 May be negative. 117 W2 = W2 - W6 - W0 118 W5 =(W5 + W2*45)/2 Now >= 0 again. 119 W4 =(W4 - W2)/3 120 W2 = W2 - W4 121 122 W1 = W5 - W1 May be negative. 123 W5 =(W5 - W3*8)/9 124 W3 = W3 - W5 125 W1 =(W1/15 + W5)/2 Now >= 0 again. 126 W5 = W5 - W1 127 128 where W0 = f(0), W1 = f(-2), W2 = f(1), W3 = f(-1), 129 W4 = f(2), W5 = f(1/2), W6 = f(oo), 130 131 Note that most intermediate results are positive; the ones that 132 may be negative are represented in two's complement. We must 133 never shift right a value that may be negative, since that would 134 invalidate the sign bit. On the other hand, divexact by odd 135 numbers work fine with two's complement. 136 */ 137 138 mpn_add_n (w5, w5, w4, m); 139 if (flags & toom7_w1_neg) 140 { 141 #ifdef HAVE_NATIVE_mpn_rsh1add_n 142 mpn_rsh1add_n (w1, w1, w4, m); 143 #else 144 mpn_add_n (w1, w1, w4, m); ASSERT (!(w1[0] & 1)); 145 mpn_rshift (w1, w1, m, 1); 146 #endif 147 } 148 else 149 { 150 #ifdef HAVE_NATIVE_mpn_rsh1sub_n 151 mpn_rsh1sub_n (w1, w4, w1, m); 152 #else 153 mpn_sub_n (w1, w4, w1, m); ASSERT (!(w1[0] & 1)); 154 mpn_rshift (w1, w1, m, 1); 155 #endif 156 } 157 mpn_sub (w4, w4, m, w0, 2*n); 158 mpn_sub_n (w4, w4, w1, m); ASSERT (!(w4[0] & 3)); 159 mpn_rshift (w4, w4, m, 2); /* w4>=0 */ 160 161 tp[w6n] = mpn_lshift (tp, w6, w6n, 4); 162 mpn_sub (w4, w4, m, tp, w6n+1); 163 164 if (flags & toom7_w3_neg) 165 { 166 #ifdef HAVE_NATIVE_mpn_rsh1add_n 167 mpn_rsh1add_n (w3, w3, w2, m); 168 #else 169 mpn_add_n (w3, w3, w2, m); ASSERT (!(w3[0] & 1)); 170 mpn_rshift (w3, w3, m, 1); 171 #endif 172 } 173 else 174 { 175 #ifdef HAVE_NATIVE_mpn_rsh1sub_n 176 mpn_rsh1sub_n (w3, w2, w3, m); 177 #else 178 mpn_sub_n (w3, w2, w3, m); ASSERT (!(w3[0] & 1)); 179 mpn_rshift (w3, w3, m, 1); 180 #endif 181 } 182 183 mpn_sub_n (w2, w2, w3, m); 184 185 mpn_submul_1 (w5, w2, m, 65); 186 mpn_sub (w2, w2, m, w6, w6n); 187 mpn_sub (w2, w2, m, w0, 2*n); 188 189 mpn_addmul_1 (w5, w2, m, 45); ASSERT (!(w5[0] & 1)); 190 mpn_rshift (w5, w5, m, 1); 191 mpn_sub_n (w4, w4, w2, m); 192 193 mpn_divexact_by3 (w4, w4, m); 194 mpn_sub_n (w2, w2, w4, m); 195 196 mpn_sub_n (w1, w5, w1, m); 197 mpn_lshift (tp, w3, m, 3); 198 mpn_sub_n (w5, w5, tp, m); 199 mpn_divexact_by9 (w5, w5, m); 200 mpn_sub_n (w3, w3, w5, m); 201 202 mpn_divexact_by15 (w1, w1, m); 203 mpn_add_n (w1, w1, w5, m); ASSERT (!(w1[0] & 1)); 204 mpn_rshift (w1, w1, m, 1); /* w1>=0 now */ 205 mpn_sub_n (w5, w5, w1, m); 206 207 /* These bounds are valid for the 4x4 polynomial product of toom44, 208 * and they are conservative for toom53 and toom62. */ 209 ASSERT (w1[2*n] < 2); 210 ASSERT (w2[2*n] < 3); 211 ASSERT (w3[2*n] < 4); 212 ASSERT (w4[2*n] < 3); 213 ASSERT (w5[2*n] < 2); 214 215 /* Addition chain. Note carries and the 2n'th limbs that need to be 216 * added in. 217 * 218 * Special care is needed for w2[2n] and the corresponding carry, 219 * since the "simple" way of adding it all together would overwrite 220 * the limb at wp[2*n] and rp[4*n] (same location) with the sum of 221 * the high half of w3 and the low half of w4. 222 * 223 * 7 6 5 4 3 2 1 0 224 * | | | | | | | | | 225 * ||w3 (2n+1)| 226 * ||w4 (2n+1)| 227 * ||w5 (2n+1)| ||w1 (2n+1)| 228 * + | w6 (w6n)| ||w2 (2n+1)| w0 (2n) | (share storage with r) 229 * ----------------------------------------------- 230 * r | | | | | | | | | 231 * c7 c6 c5 c4 c3 Carries to propagate 232 */ 233 234 cy = mpn_add_n (rp + n, rp + n, w1, m); 235 MPN_INCR_U (w2 + n + 1, n , cy); 236 cy = mpn_add_n (rp + 3*n, rp + 3*n, w3, n); 237 MPN_INCR_U (w3 + n, n + 1, w2[2*n] + cy); 238 cy = mpn_add_n (rp + 4*n, w3 + n, w4, n); 239 MPN_INCR_U (w4 + n, n + 1, w3[2*n] + cy); 240 cy = mpn_add_n (rp + 5*n, w4 + n, w5, n); 241 MPN_INCR_U (w5 + n, n + 1, w4[2*n] + cy); 242 if (w6n > n + 1) 243 ASSERT_NOCARRY (mpn_add (rp + 6*n, rp + 6*n, w6n, w5 + n, n + 1)); 244 else 245 { 246 ASSERT_NOCARRY (mpn_add_n (rp + 6*n, rp + 6*n, w5 + n, w6n)); 247 #if WANT_ASSERT 248 { 249 mp_size_t i; 250 for (i = w6n; i <= n; i++) 251 ASSERT (w5[n + i] == 0); 252 } 253 #endif 254 } 255 } 256