1 /* mpn_toom_eval_pm2exp -- Evaluate a polynomial in +2^k and -2^k
3 Contributed to the GNU project by Niels Möller
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.
9 Copyright 2009 Free Software Foundation, Inc.
11 This file is part of the GNU MP Library.
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of the GNU Lesser General Public License as published by
15 the Free Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
18 The GNU MP Library is distributed in the hope that it will be useful, but
19 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
20 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21 License for more details.
23 You should have received a copy of the GNU Lesser General Public License
24 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
30 /* Evaluates a polynomial of degree k > 2, in the points +2^shift and -2^shift. */
32 mpn_toom_eval_pm2exp (mp_ptr xp2, mp_ptr xm2, unsigned k,
33 mp_srcptr xp, mp_size_t n, mp_size_t hn, unsigned shift,
38 #if HAVE_NATIVE_mpn_addlsh_n
43 ASSERT (shift*k < GMP_NUMB_BITS);
48 /* The degree k is also the number of full-size coefficients, so
49 * that last coefficient, of size hn, starts at xp + k*n. */
51 #if HAVE_NATIVE_mpn_addlsh_n
52 xp2[n] = mpn_addlsh_n (xp2, xp, xp + 2*n, n, 2*shift);
53 for (i = 4; i < k; i += 2)
54 xp2[n] += mpn_addlsh_n (xp2, xp2, xp + i*n, n, i*shift);
56 tp[n] = mpn_lshift (tp, xp+n, n, shift);
57 for (i = 3; i < k; i+= 2)
58 tp[n] += mpn_addlsh_n (tp, tp, xp+i*n, n, i*shift);
62 cy = mpn_addlsh_n (tp, tp, xp+k*n, hn, k*shift);
63 MPN_INCR_U (tp + hn, n+1 - hn, cy);
67 cy = mpn_addlsh_n (xp2, xp2, xp+k*n, hn, k*shift);
68 MPN_INCR_U (xp2 + hn, n+1 - hn, cy);
71 #else /* !HAVE_NATIVE_mpn_addlsh_n */
72 xp2[n] = mpn_lshift (tp, xp+2*n, n, 2*shift);
73 xp2[n] += mpn_add_n (xp2, xp, tp, n);
74 for (i = 4; i < k; i += 2)
76 xp2[n] += mpn_lshift (tp, xp + i*n, n, i*shift);
77 xp2[n] += mpn_add_n (xp2, xp2, tp, n);
80 tp[n] = mpn_lshift (tp, xp+n, n, shift);
81 for (i = 3; i < k; i+= 2)
83 tp[n] += mpn_lshift (xm2, xp + i*n, n, i*shift);
84 tp[n] += mpn_add_n (tp, tp, xm2, n);
87 xm2[hn] = mpn_lshift (xm2, xp + k*n, hn, k*shift);
89 mpn_add (tp, tp, n+1, xm2, hn+1);
91 mpn_add (xp2, xp2, n+1, xm2, hn+1);
92 #endif /* !HAVE_NATIVE_mpn_addlsh_n */
94 neg = (mpn_cmp (xp2, tp, n + 1) < 0) ? ~0 : 0;
96 #if HAVE_NATIVE_mpn_add_n_sub_n
98 mpn_add_n_sub_n (xp2, xm2, tp, xp2, n + 1);
100 mpn_add_n_sub_n (xp2, xm2, xp2, tp, n + 1);
101 #else /* !HAVE_NATIVE_mpn_add_n_sub_n */
103 mpn_sub_n (xm2, tp, xp2, n + 1);
105 mpn_sub_n (xm2, xp2, tp, n + 1);
107 mpn_add_n (xp2, xp2, tp, n + 1);
108 #endif /* !HAVE_NATIVE_mpn_add_n_sub_n */
110 /* FIXME: the following asserts are useless if (k+1)*shift >= GMP_LIMB_BITS */
111 ASSERT ((k+1)*shift >= GMP_LIMB_BITS ||
112 xp2[n] < ((CNST_LIMB(1)<<((k+1)*shift))-1)/((CNST_LIMB(1)<<shift)-1));
113 ASSERT ((k+2)*shift >= GMP_LIMB_BITS ||
114 xm2[n] < ((CNST_LIMB(1)<<((k+2)*shift))-((k&1)?(CNST_LIMB(1)<<shift):1))/((CNST_LIMB(1)<<(2*shift))-1));