toom_eval_pm2exp.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom_eval_pm2exp -- Evaluate a polynomial in +2^k and -2^k
  2.    Contributed to the GNU project by Niels M鰈ler
  3.    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
  4.    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  5.    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  6. Copyright 2009 Free Software Foundation, Inc.
  7. This file is part of the GNU MP Library.
  8. The GNU MP Library is free software; you can redistribute it and/or modify
  9. it under the terms of the GNU Lesser General Public License as published by
  10. the Free Software Foundation; either version 3 of the License, or (at your
  11. option) any later version.
  12. The GNU MP Library is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15. License for more details.
  16. You should have received a copy of the GNU Lesser General Public License
  17. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  18. #include "gmp.h"
  19. #include "gmp-impl.h"
  20. /* Evaluates a polynomial of degree k > 2, in the points +2^shift and -2^shift. */
  21. int
  22. mpn_toom_eval_pm2exp (mp_ptr xp2, mp_ptr xm2, unsigned k,
  23.       mp_srcptr xp, mp_size_t n, mp_size_t hn, unsigned shift,
  24.       mp_ptr tp)
  25. {
  26.   unsigned i;
  27.   int neg;
  28. #if HAVE_NATIVE_mpn_addlsh_n
  29.   mp_limb_t cy;
  30. #endif
  31.   ASSERT (k >= 3);
  32.   ASSERT (shift*k < GMP_NUMB_BITS);
  33.   ASSERT (hn > 0);
  34.   ASSERT (hn <= n);
  35.   /* The degree k is also the number of full-size coefficients, so
  36.    * that last coefficient, of size hn, starts at xp + k*n. */
  37. #if HAVE_NATIVE_mpn_addlsh_n
  38.   xp2[n] = mpn_addlsh_n (xp2, xp, xp + 2*n, n, 2*shift);
  39.   for (i = 4; i < k; i += 2)
  40.     xp2[n] += mpn_addlsh_n (xp2, xp2, xp + i*n, n, i*shift);
  41.   tp[n] = mpn_lshift (tp, xp+n, n, shift);
  42.   for (i = 3; i < k; i+= 2)
  43.     tp[n] += mpn_addlsh_n (tp, tp, xp+i*n, n, i*shift);
  44.   if (k & 1)
  45.     {
  46.       cy = mpn_addlsh_n (tp, tp, xp+k*n, hn, k*shift);
  47.       MPN_INCR_U (tp + hn, n+1 - hn, cy);
  48.     }
  49.   else
  50.     {
  51.       cy = mpn_addlsh_n (xp2, xp2, xp+k*n, hn, k*shift);
  52.       MPN_INCR_U (xp2 + hn, n+1 - hn, cy);
  53.     }
  54. #else /* !HAVE_NATIVE_mpn_addlsh_n */
  55.   xp2[n] = mpn_lshift (tp, xp+2*n, n, 2*shift);
  56.   xp2[n] += mpn_add_n (xp2, xp, tp, n);
  57.   for (i = 4; i < k; i += 2)
  58.     {
  59.       xp2[n] += mpn_lshift (tp, xp + i*n, n, i*shift);
  60.       xp2[n] += mpn_add_n (xp2, xp2, tp, n);
  61.     }
  62.   tp[n] = mpn_lshift (tp, xp+n, n, shift);
  63.   for (i = 3; i < k; i+= 2)
  64.     {
  65.       tp[n] += mpn_lshift (xm2, xp + i*n, n, i*shift);
  66.       tp[n] += mpn_add_n (tp, tp, xm2, n);
  67.     }
  68.   xm2[hn] = mpn_lshift (xm2, xp + k*n, hn, k*shift);
  69.   if (k & 1)
  70.     mpn_add (tp, tp, n+1, xm2, hn+1);
  71.   else
  72.     mpn_add (xp2, xp2, n+1, xm2, hn+1);
  73. #endif /* !HAVE_NATIVE_mpn_addlsh_n */
  74.   neg = (mpn_cmp (xp2, tp, n + 1) < 0) ? ~0 : 0;
  75. #if HAVE_NATIVE_mpn_add_n_sub_n
  76.   if (neg)
  77.     mpn_add_n_sub_n (xp2, xm2, tp, xp2, n + 1);
  78.   else
  79.     mpn_add_n_sub_n (xp2, xm2, xp2, tp, n + 1);
  80. #else /* !HAVE_NATIVE_mpn_add_n_sub_n */
  81.   if (neg)
  82.     mpn_sub_n (xm2, tp, xp2, n + 1);
  83.   else
  84.     mpn_sub_n (xm2, xp2, tp, n + 1);
  85.   mpn_add_n (xp2, xp2, tp, n + 1);
  86. #endif /* !HAVE_NATIVE_mpn_add_n_sub_n */
  87.   /* FIXME: the following asserts are useless if (k+1)*shift >= GMP_LIMB_BITS */
  88.   ASSERT ((k+1)*shift >= GMP_LIMB_BITS ||
  89.   xp2[n] < ((CNST_LIMB(1)<<((k+1)*shift))-1)/((CNST_LIMB(1)<<shift)-1));
  90.   ASSERT ((k+2)*shift >= GMP_LIMB_BITS ||
  91.   xm2[n] < ((CNST_LIMB(1)<<((k+2)*shift))-((k&1)?(CNST_LIMB(1)<<shift):1))/((CNST_LIMB(1)<<(2*shift))-1));
  92.   return neg;
  93. }