source/libs/gmp/gmp-src/mpn/generic/toom_eval_pm2.c

   1 /* mpn_toom_eval_pm2 -- Evaluate a polynomial in +2 and -2
   2
   3    Contributed to the GNU project by Niels Möller and 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 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 #include "gmp.h"
  38 #include "gmp-impl.h"
  39
  40 /* DO_addlsh2(d,a,b,n,cy) computes cy,{d,n} <- {a,n} + 4*(cy,{b,n}), it
  41    can be used as DO_addlsh2(d,a,d,n,d[n]), for accumulation on {d,n+1}. */
  42 #if HAVE_NATIVE_mpn_addlsh2_n
  43 #define DO_addlsh2(d, a, b, n, cy)      \
  44 do {                                    \
  45   (cy) <<= 2;                           \
  46   (cy) += mpn_addlsh2_n(d, a, b, n);    \
  47 } while (0)
  48 #else
  49 #if HAVE_NATIVE_mpn_addlsh_n
  50 #define DO_addlsh2(d, a, b, n, cy)      \
  51 do {                                    \
  52   (cy) <<= 2;                           \
  53   (cy) += mpn_addlsh_n(d, a, b, n, 2);  \
  54 } while (0)
  55 #else
  56 /* The following is not a general substitute for addlsh2.
  57    It is correct if d == b, but it is not if d == a.  */
  58 #define DO_addlsh2(d, a, b, n, cy)      \
  59 do {                                    \
  60   (cy) <<= 2;                           \
  61   (cy) += mpn_lshift(d, b, n, 2);       \
  62   (cy) += mpn_add_n(d, d, a, n);        \
  63 } while (0)
  64 #endif
  65 #endif
  66
  67 /* Evaluates a polynomial of degree 2 < k < GMP_NUMB_BITS, in the
  68    points +2 and -2. */
  69 int
  70 mpn_toom_eval_pm2 (mp_ptr xp2, mp_ptr xm2, unsigned k,
  71                    mp_srcptr xp, mp_size_t n, mp_size_t hn, mp_ptr tp)
  72 {
  73   int i;
  74   int neg;
  75   mp_limb_t cy;
  76
  77   ASSERT (k >= 3);
  78   ASSERT (k < GMP_NUMB_BITS);
  79
  80   ASSERT (hn > 0);
  81   ASSERT (hn <= n);
  82
  83   /* The degree k is also the number of full-size coefficients, so
  84    * that last coefficient, of size hn, starts at xp + k*n. */
  85
  86   cy = 0;
  87   DO_addlsh2 (xp2, xp + (k-2) * n, xp + k * n, hn, cy);
  88   if (hn != n)
  89     cy = mpn_add_1 (xp2 + hn, xp + (k-2) * n + hn, n - hn, cy);
  90   for (i = k - 4; i >= 0; i -= 2)
  91     DO_addlsh2 (xp2, xp + i * n, xp2, n, cy);
  92   xp2[n] = cy;
  93
  94   k--;
  95
  96   cy = 0;
  97   DO_addlsh2 (tp, xp + (k-2) * n, xp + k * n, n, cy);
  98   for (i = k - 4; i >= 0; i -= 2)
  99     DO_addlsh2 (tp, xp + i * n, tp, n, cy);
 100   tp[n] = cy;
 101
 102   if (k & 1)
 103     ASSERT_NOCARRY(mpn_lshift (tp , tp , n + 1, 1));
 104   else
 105     ASSERT_NOCARRY(mpn_lshift (xp2, xp2, n + 1, 1));
 106
 107   neg = (mpn_cmp (xp2, tp, n + 1) < 0) ? ~0 : 0;
 108
 109 #if HAVE_NATIVE_mpn_add_n_sub_n
 110   if (neg)
 111     mpn_add_n_sub_n (xp2, xm2, tp, xp2, n + 1);
 112   else
 113     mpn_add_n_sub_n (xp2, xm2, xp2, tp, n + 1);
 114 #else /* !HAVE_NATIVE_mpn_add_n_sub_n */
 115   if (neg)
 116     mpn_sub_n (xm2, tp, xp2, n + 1);
 117   else
 118     mpn_sub_n (xm2, xp2, tp, n + 1);
 119
 120   mpn_add_n (xp2, xp2, tp, n + 1);
 121 #endif /* !HAVE_NATIVE_mpn_add_n_sub_n */
 122
 123   ASSERT (xp2[n] < (1<<(k+2))-1);
 124   ASSERT (xm2[n] < ((1<<(k+3))-1 - (1^k&1))/3);
 125
 126   neg ^= ((k & 1) - 1);
 127
 128   return neg;
 129 }
 130
 131 #undef DO_addlsh2