/src/libgmp/mpn/toom_eval_pm2rexp.c
Line | Count | Source |
1 | | /* mpn_toom_eval_pm2rexp -- Evaluate a polynomial in +2^-k and -2^-k |
2 | | |
3 | | Contributed to the GNU project by 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 | | |
38 | | #include "gmp-impl.h" |
39 | | |
40 | | #if HAVE_NATIVE_mpn_addlsh_n |
41 | 18.5k | #define DO_mpn_addlsh_n(dst,src,n,s,ws) mpn_addlsh_n(dst,dst,src,n,s) |
42 | | #else |
43 | | static mp_limb_t |
44 | | DO_mpn_addlsh_n(mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws) |
45 | | { |
46 | | #if USE_MUL_1 && 0 |
47 | | return mpn_addmul_1(dst,src,n,CNST_LIMB(1) <<(s)); |
48 | | #else |
49 | | mp_limb_t __cy; |
50 | | __cy = mpn_lshift(ws,src,n,s); |
51 | | return __cy + mpn_add_n(dst,dst,ws,n); |
52 | | #endif |
53 | | } |
54 | | #endif |
55 | | |
56 | | /* Evaluates a polynomial of degree k >= 3. */ |
57 | | /* It returns 0 or ~0, depending on the sign of the result rm. */ |
58 | | unsigned |
59 | | mpn_toom_eval_pm2rexp (mp_ptr rp, mp_ptr rm, |
60 | | unsigned int q, mp_srcptr ap, mp_size_t n, mp_size_t t, |
61 | | unsigned int s, mp_ptr ws) |
62 | 4.11k | { |
63 | 4.11k | unsigned int i; |
64 | 4.11k | unsigned neg; |
65 | | /* {ap,q*n+t} -> {rp,n+1} {rm,n+1} , with {ws, n+1}*/ |
66 | 4.11k | ASSERT (n >= t); |
67 | 4.11k | ASSERT (s != 0); /* or _eval_pm1 should be used */ |
68 | 4.11k | ASSERT (q > 1); |
69 | 4.11k | ASSERT (s*q < GMP_NUMB_BITS); |
70 | 4.11k | rp[n] = mpn_lshift(rp, ap, n, s*q); |
71 | 4.11k | ws[n] = mpn_lshift(ws, ap+n, n, s*(q-1)); |
72 | 4.11k | if( (q & 1) != 0) { |
73 | 4.09k | ASSERT_NOCARRY(mpn_add(ws,ws,n+1,ap+n*q,t)); |
74 | 4.09k | rp[n] += DO_mpn_addlsh_n(rp, ap+n*(q-1), n, s, rm); |
75 | 4.09k | } else { |
76 | 23 | ASSERT_NOCARRY(mpn_add(rp,rp,n+1,ap+n*q,t)); |
77 | 23 | } |
78 | 11.3k | for(i=2; i<q-1; i++) |
79 | 7.22k | { |
80 | 7.22k | rp[n] += DO_mpn_addlsh_n(rp, ap+n*i, n, s*(q-i), rm); |
81 | 7.22k | i++; |
82 | 7.22k | ws[n] += DO_mpn_addlsh_n(ws, ap+n*i, n, s*(q-i), rm); |
83 | 7.22k | }; |
84 | | |
85 | 4.11k | neg = - (unsigned) (mpn_cmp (rp, ws, n + 1) < 0); |
86 | | |
87 | | #if HAVE_NATIVE_mpn_add_n_sub_n |
88 | | if (neg) |
89 | | mpn_add_n_sub_n (rp, rm, ws, rp, n + 1); |
90 | | else |
91 | | mpn_add_n_sub_n (rp, rm, rp, ws, n + 1); |
92 | | #else /* !HAVE_NATIVE_mpn_add_n_sub_n */ |
93 | 4.11k | if (neg) |
94 | 564 | mpn_sub_n (rm, ws, rp, n + 1); |
95 | 3.55k | else |
96 | 3.55k | mpn_sub_n (rm, rp, ws, n + 1); |
97 | | |
98 | 4.11k | ASSERT_NOCARRY (mpn_add_n (rp, rp, ws, n + 1)); |
99 | 4.11k | #endif /* !HAVE_NATIVE_mpn_add_n_sub_n */ |
100 | | |
101 | 4.11k | return neg; |
102 | 4.11k | } |