/src/libgmp/mpn/toom_eval_pm2exp.c
Line | Count | Source |
1 | | /* mpn_toom_eval_pm2exp -- Evaluate a polynomial in +2^k and -2^k |
2 | | |
3 | | Contributed to the GNU project by Niels Möller |
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 | | |
38 | | #include "gmp-impl.h" |
39 | | |
40 | | /* Evaluates a polynomial of degree k > 2, in the points +2^shift and -2^shift. */ |
41 | | /* It returns 0 or ~0, depending on the sign of the result xm2. */ |
42 | | unsigned |
43 | | mpn_toom_eval_pm2exp (mp_ptr xp2, mp_ptr xm2, unsigned k, |
44 | | mp_srcptr xp, mp_size_t n, mp_size_t hn, unsigned shift, |
45 | | mp_ptr tp) |
46 | 2.72k | { |
47 | 2.72k | unsigned i; |
48 | 2.72k | unsigned neg; |
49 | 2.72k | #if HAVE_NATIVE_mpn_addlsh_n |
50 | 2.72k | mp_limb_t cy; |
51 | 2.72k | #endif |
52 | | |
53 | 2.72k | ASSERT (k >= 3); |
54 | 2.72k | ASSERT (shift*k < GMP_NUMB_BITS); |
55 | | |
56 | 2.72k | ASSERT (hn > 0); |
57 | 2.72k | ASSERT (hn <= n); |
58 | | |
59 | | /* The degree k is also the number of full-size coefficients, so |
60 | | * that last coefficient, of size hn, starts at xp + k*n. */ |
61 | | |
62 | 2.72k | #if HAVE_NATIVE_mpn_addlsh_n |
63 | 2.72k | xp2[n] = mpn_addlsh_n (xp2, xp, xp + 2*n, n, 2*shift); |
64 | 7.51k | for (i = 4; i < k; i += 2) |
65 | 4.78k | xp2[n] += mpn_addlsh_n (xp2, xp2, xp + i*n, n, i*shift); |
66 | | |
67 | 2.72k | tp[n] = mpn_lshift (tp, xp+n, n, shift); |
68 | 7.52k | for (i = 3; i < k; i+= 2) |
69 | 4.79k | tp[n] += mpn_addlsh_n (tp, tp, xp+i*n, n, i*shift); |
70 | | |
71 | 2.72k | if (k & 1) |
72 | 2.71k | { |
73 | 2.71k | cy = mpn_addlsh_n (tp, tp, xp+k*n, hn, k*shift); |
74 | 2.71k | MPN_INCR_U (tp + hn, n+1 - hn, cy); |
75 | 2.71k | } |
76 | 12 | else |
77 | 12 | { |
78 | 12 | cy = mpn_addlsh_n (xp2, xp2, xp+k*n, hn, k*shift); |
79 | 12 | MPN_INCR_U (xp2 + hn, n+1 - hn, cy); |
80 | 12 | } |
81 | | |
82 | | #else /* !HAVE_NATIVE_mpn_addlsh_n */ |
83 | | xp2[n] = mpn_lshift (tp, xp+2*n, n, 2*shift); |
84 | | xp2[n] += mpn_add_n (xp2, xp, tp, n); |
85 | | for (i = 4; i < k; i += 2) |
86 | | { |
87 | | xp2[n] += mpn_lshift (tp, xp + i*n, n, i*shift); |
88 | | xp2[n] += mpn_add_n (xp2, xp2, tp, n); |
89 | | } |
90 | | |
91 | | tp[n] = mpn_lshift (tp, xp+n, n, shift); |
92 | | for (i = 3; i < k; i+= 2) |
93 | | { |
94 | | tp[n] += mpn_lshift (xm2, xp + i*n, n, i*shift); |
95 | | tp[n] += mpn_add_n (tp, tp, xm2, n); |
96 | | } |
97 | | |
98 | | xm2[hn] = mpn_lshift (xm2, xp + k*n, hn, k*shift); |
99 | | if (k & 1) |
100 | | mpn_add (tp, tp, n+1, xm2, hn+1); |
101 | | else |
102 | | mpn_add (xp2, xp2, n+1, xm2, hn+1); |
103 | | #endif /* !HAVE_NATIVE_mpn_addlsh_n */ |
104 | | |
105 | 2.72k | neg = - (unsigned) (mpn_cmp (xp2, tp, n + 1) < 0); |
106 | | |
107 | | #if HAVE_NATIVE_mpn_add_n_sub_n |
108 | | if (neg) |
109 | | mpn_add_n_sub_n (xp2, xm2, tp, xp2, n + 1); |
110 | | else |
111 | | mpn_add_n_sub_n (xp2, xm2, xp2, tp, n + 1); |
112 | | #else /* !HAVE_NATIVE_mpn_add_n_sub_n */ |
113 | 2.72k | if (neg) |
114 | 355 | mpn_sub_n (xm2, tp, xp2, n + 1); |
115 | 2.37k | else |
116 | 2.37k | mpn_sub_n (xm2, xp2, tp, n + 1); |
117 | | |
118 | 2.72k | mpn_add_n (xp2, xp2, tp, n + 1); |
119 | 2.72k | #endif /* !HAVE_NATIVE_mpn_add_n_sub_n */ |
120 | | |
121 | | /* FIXME: the following asserts are useless if (k+1)*shift >= GMP_LIMB_BITS */ |
122 | 2.72k | ASSERT ((k+1)*shift >= GMP_LIMB_BITS || |
123 | 2.72k | xp2[n] < ((CNST_LIMB(1)<<((k+1)*shift))-1)/((CNST_LIMB(1)<<shift)-1)); |
124 | 2.72k | ASSERT ((k+2)*shift >= GMP_LIMB_BITS || |
125 | 2.72k | xm2[n] < ((CNST_LIMB(1)<<((k+2)*shift))-((k&1)?(CNST_LIMB(1)<<shift):1))/((CNST_LIMB(1)<<(2*shift))-1)); |
126 | | |
127 | 2.72k | return neg; |
128 | 2.72k | } |