1 /*
2 mpcx_rem
3
4 computes the remainder of f divided by g
5
6 Copyright (C) 2009 Andreas Enge
7
8 This file is part of the MPFRCX Library.
9
10 The MPFRCX Library is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 3 of the License, or (at your
13 option) any later version.
14
15 The MPFRCX Library is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 License for more details.
19
20 You should have received a copy of the GNU Lesser General Public License
21 along with the MPFRCX library; see the file COPYING.LESSER. If not, write to
22 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
23 MA 02111-1307, USA.
24 */
25
26 #include "mpfrcx-impl.h"
27
28 static void mpcx_rem_naive (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g);
29 static void mpcx_rem_newton (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g);
30 static void mpcx_rev (mpcx_ptr h, mpcx_srcptr f, const int exp);
31 static void mpcx_smul (mpcx_ptr h, mpcx_ptr f, mpcx_ptr g, const int exp);
32 static void mpcx_sinv (mpcx_ptr h, mpcx_ptr f, const int exp);
33
34 /*****************************************************************************/
35
mpcx_rem(mpcx_ptr r,mpcx_srcptr f,mpcx_srcptr g)36 void mpcx_rem (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g)
37
38 {
39 if (f->deg < g->deg)
40 mpcx_set (r, f);
41 else
42 if (g->deg >= 32 && f->deg >= 63)
43 mpcx_rem_newton (r, f, g);
44 else
45 mpcx_rem_naive (r, f, g);
46 }
47
48 /*****************************************************************************/
49
mpcx_rem_naive(mpcx_ptr r,mpcx_srcptr f,mpcx_srcptr g)50 static void mpcx_rem_naive (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g)
51 /* computes the remainder of f divided by g by the naive algorithm for */
52 /* deg f >= deg g */
53 {
54 mpc_t tmp;
55 int i, j;
56
57 if (mpc_cmp_si (g->coeff [g->deg], 1) == 0) {
58 mpc_init2 (tmp, mpc_get_prec (f->coeff [0]));
59 /* the only case we actually use for multievaluation */
60 mpcx_set (r, f);
61 for (i = f->deg - g->deg; i >= 0; i--)
62 for (j = 0; j < g->deg; j++) {
63 mpc_mul (tmp, r->coeff [g->deg + i], g->coeff [j], MPC_RNDNN);
64 mpc_sub (r->coeff [i + j], r->coeff [i + j], tmp, MPC_RNDNN);
65 }
66 r->deg = g->deg - 1;
67 /* To save memory, one might want to do a realloc on r here:
68 mpcx_realloc (r, r->deg + 1);
69 */
70 mpc_clear (tmp);
71 }
72 else {
73 printf ("*** Houston, we have a problem!\n");
74 printf ("*** Calling mpcx_rem_naive with a non-monic divisor.\n");
75 printf ("*** Go back programming!\n");
76 exit (1);
77 }
78 }
79
80 /*****************************************************************************/
81
mpcx_rem_newton(mpcx_ptr r,mpcx_srcptr f,mpcx_srcptr g)82 static void mpcx_rem_newton (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g) {
83 /* computes the remainder of f divided by g by Newton iterations for */
84 /* deg f >= deg g */
85
86 mpcx_t q, tmp;
87
88 mpcx_init (q, f->deg - g->deg + 1, r->prec);
89 mpcx_init (tmp, f->deg + 1, r->prec);
90
91 mpcx_rev (tmp, f, 0);
92 mpcx_rev (q, g, 0);
93 mpcx_sinv (q, q, f->deg - g->deg + 1);
94 mpcx_smul (q, tmp, q, f->deg - g->deg + 1);
95 mpcx_rev (q, q, f->deg - g->deg);
96
97 mpcx_mul (tmp, q, g);
98 mpcx_sub (r, f, tmp);
99 if (r->deg >= g->deg)
100 r->deg = g->deg - 1;
101
102 mpcx_clear (q);
103 mpcx_clear (tmp);
104 }
105
106 /*****************************************************************************/
107
mpcx_rev(mpcx_ptr h,mpcx_srcptr f,const int exp)108 static void mpcx_rev (mpcx_ptr h, mpcx_srcptr f, const int exp) {
109 int overlap, exp_local, i;
110 mpcx_t h_local;
111
112 if (exp == 0)
113 exp_local = f->deg;
114 else
115 exp_local = exp;
116
117 if (exp_local < f->deg) {
118 printf ("*** Computing a reverse polynomial of too low ");
119 printf ("order in 'mpcx_rev'.\n");
120 exit (1);
121 }
122
123 overlap = (h == f);
124 if (overlap)
125 mpcx_init (h_local, exp_local + 1, h->prec);
126 else {
127 mpcx_mv (h_local, h);
128 if (h_local->size < exp_local + 1)
129 mpcx_realloc (h_local, exp_local + 1);
130 }
131
132 h_local->deg = exp_local;
133 for (i = 0; i <= f->deg; i++)
134 mpc_set (h_local->coeff [exp_local - i], f->coeff [i], MPC_RNDNN);
135 for (i = exp_local - f->deg - 1; i >= 0; i--)
136 mpc_set_ui (h_local->coeff [i], 0, MPC_RNDNN);
137
138 while (h_local->deg >= 0
139 && mpc_cmp_si (h_local->coeff [h_local->deg], 0) == 0)
140 h_local->deg--;
141 if (overlap)
142 mpcx_clear (h);
143 mpcx_mv (h, h_local);
144 }
145
146 /*****************************************************************************/
147
mpcx_smul(mpcx_ptr h,mpcx_ptr f,mpcx_ptr g,const int exp)148 static void mpcx_smul (mpcx_ptr h, mpcx_ptr f, mpcx_ptr g, const int exp) {
149 /* computes the short product h = f*g mod X^exp; for the time being, */
150 /* by a complete product and a truncation */
151
152 int f_deg, g_deg, h_deg;
153
154 f_deg = f->deg;
155 g_deg = g->deg;
156 if (f_deg >= exp)
157 f->deg = exp - 1;
158 if (g_deg >= exp)
159 g->deg = exp - 1;
160 mpcx_mul (h, f, g);
161 h_deg = h->deg;
162 f->deg = f_deg;
163 g->deg = g_deg;
164 h->deg = h_deg;
165 /* The previous line is needed if h is one of f or g. */
166 if (h->deg >= exp) {
167 h->deg = exp - 1;
168 while (h->deg >= 0 && mpc_cmp_si (h->coeff [h->deg], 0) == 0)
169 h->deg--;
170 }
171 }
172
173 /*****************************************************************************/
174
mpcx_sinv(mpcx_ptr h,mpcx_ptr f,const int exp)175 static void mpcx_sinv (mpcx_ptr h, mpcx_ptr f, const int exp) {
176 /* computes h = f^{-1} mod X^exp by Newton iteration */
177
178 int overlap;
179 mpcx_t h_local, tmp;
180 int exp_local;
181
182 overlap = (h == f);
183 if (overlap)
184 mpcx_init (h_local, exp, h->prec);
185 else
186 mpcx_mv (h_local, h);
187 mpcx_init (tmp, exp, h->prec);
188
189 exp_local = 1;
190 h_local->deg = 0;
191 mpc_ui_div (h_local->coeff [0], 1, f->coeff [0], MPC_RNDNN);
192 while (exp_local < exp) {
193 exp_local = (2 * exp_local < exp ? 2 * exp_local : exp);
194 mpcx_smul (tmp, h_local, f, exp_local);
195 mpcx_si_sub (tmp, 2, tmp);
196 mpcx_smul (h_local, h_local, tmp, exp_local);
197 }
198
199 if (overlap)
200 mpcx_clear (h);
201 mpcx_clear (tmp);
202 mpcx_mv (h, h_local);
203 }
204