contrib/gmp/mpn/generic/toom_interpolate_7pts.c

   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 }