Line data Source code
1 : #line 2 "../src/kernel/none/ratlift.c"
2 : /* Copyright (C) 2003 The PARI group.
3 :
4 : This file is part of the PARI/GP package.
5 :
6 : PARI/GP is free software; you can redistribute it and/or modify it under the
7 : terms of the GNU General Public License as published by the Free Software
8 : Foundation; either version 2 of the License, or (at your option) any later
9 : version. It is distributed in the hope that it will be useful, but WITHOUT
10 : ANY WARRANTY WHATSOEVER.
11 :
12 : Check the License for details. You should have received a copy of it, along
13 : with the package; see the file 'COPYING'. If not, write to the Free Software
14 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
15 :
16 : /*==========================================================
17 : * Fp_ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)
18 : *==========================================================
19 : * Reconstruct rational number from its residue x mod m
20 : * Given t_INT x, m, amax>=0, bmax>0 such that
21 : * 0 <= x < m; 2*amax*bmax < m
22 : * attempts to find t_INT a, b such that
23 : * (1) a = b*x (mod m)
24 : * (2) |a| <= amax, 0 < b <= bmax
25 : * (3) gcd(m, b) = gcd(a, b)
26 : * If unsuccessful, it will return 0 and leave a,b unchanged (and
27 : * caller may deduce no such a,b exist). If successful, sets a,b
28 : * and returns 1. If there exist a,b satisfying (1), (2), and
29 : * (3') gcd(m, b) = 1
30 : * then they are uniquely determined subject to (1),(2) and
31 : * (3'') gcd(a, b) = 1,
32 : * and will be returned by the routine. (The caller may wish to
33 : * check gcd(a,b)==1, either directly or based on known prime
34 : * divisors of m, depending on the application.)
35 : * Reference:
36 : @article {MR97c:11116,
37 : AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.},
38 : TITLE = {Efficient rational number reconstruction},
39 : JOURNAL = {J. Symbolic Comput.},
40 : VOLUME = {20},
41 : YEAR = {1995},
42 : NUMBER = {3},
43 : PAGES = {287--297},
44 : }
45 : * Preprint available from:
46 : * ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz */
47 : static ulong
48 109760 : get_vmax(GEN r, long lb, long lbb)
49 : {
50 109760 : long lr = lb - lgefint(r);
51 : ulong vmax;
52 109760 : if (lr > 1) /* still more than a word's worth to go */
53 0 : vmax = ULONG_MAX; /* (cannot in fact happen) */
54 : else
55 : { /* take difference of bit lengths */
56 109760 : long lbr = bfffo(*int_MSW(r));
57 109760 : lr = lr*BITS_IN_LONG - lbb + lbr;
58 109760 : if ((ulong)lr > BITS_IN_LONG)
59 7678 : vmax = ULONG_MAX;
60 102082 : else if (lr == 0)
61 8663 : vmax = 1UL;
62 : else
63 93419 : vmax = 1UL << (lr-1); /* pessimistic but faster than a division */
64 : }
65 109760 : return vmax;
66 : }
67 :
68 : /* assume bmax <= sqrt(m), fast if amax <=sqrt(m) */
69 : static int
70 4817392 : Fp_ratlift_hgcd(GEN n, GEN m, GEN amax, GEN bmax, GEN *pa, GEN *pb)
71 : {
72 4817392 : pari_sp av = avma;
73 : GEN x, y, a, b;
74 4817392 : GEN H = halfgcdii(n, m), M = gel(H,1), V = gel(H,2);
75 4817777 : x = gel(V,1); a = gel(V,2); y = gcoeff(M,1,1); b = gcoeff(M,2,1);
76 5006152 : while(abscmpii(b, bmax)<=0)
77 : {
78 : GEN q, r, u;
79 4689635 : if (abscmpii(a, amax)<=0)
80 : {
81 4501143 : if (signe(b)<0) { a = negi(a); b = negi(b); }
82 4500938 : *pa =a; *pb = b;
83 4500938 : gerepileall(av, 2, pb, pa); return 1;
84 : }
85 188367 : q = dvmdii(x, a, &r); x = a; a = r;
86 188376 : u = subii(y, mulii(b, q));
87 188375 : y = b; b = u;
88 : }
89 316395 : return gc_bool(av, 0);
90 : }
91 :
92 : /* Assume x,m,amax >= 0,bmax > 0 are t_INTs, 0 <= x < m, 2 amax * bmax < m */
93 : int
94 21177384 : Fp_ratlift(GEN x, GEN m, GEN amax, GEN bmax, GEN *a, GEN *b)
95 : {
96 : GEN d, d1, v, v1, q, r;
97 21177384 : pari_sp av = avma, av1;
98 : long lb, lbb, s, s0;
99 : ulong vmax;
100 : ulong xu, xu1, xv, xv1; /* Lehmer stage recurrence matrix */
101 : int lhmres; /* Lehmer stage return value */
102 :
103 : /* special cases x=0 and/or amax=0 */
104 21177384 : if (!signe(x)) { *a = gen_0; *b = gen_1; return 1; }
105 8880062 : if (!signe(amax)) return 0;
106 : /* assert: m > x > 0, amax > 0 */
107 :
108 : /* check whether a=x, b=1 is a solution */
109 8880062 : if (cmpii(x,amax) <= 0) { *a = icopy(x); *b = gen_1; return 1; }
110 :
111 4903630 : if (amax == bmax || equalii(amax, bmax))
112 4817111 : return Fp_ratlift_hgcd(x, m, amax, bmax, a, b);
113 :
114 : /* There is no special case for single-word numbers since this is
115 : * mainly meant to be used with large moduli. */
116 86519 : (void)new_chunk(lgefint(bmax) + lgefint(amax)); /* room for a,b */
117 86518 : d = m; d1 = x;
118 86518 : v = gen_0; v1 = gen_1;
119 : /* assert d1 > amax, v1 <= bmax here */
120 86518 : lb = lgefint(bmax);
121 86518 : lbb = bfffo(*int_MSW(bmax));
122 86518 : s = 1;
123 86518 : av1 = avma;
124 :
125 : /* General case: Euclidean division chain starting with m div x, and
126 : * with bounds on the sequence of convergents' denoms v_j.
127 : * Just to be different from what invmod and bezout are doing, we work
128 : * here with the all-nonnegative matrices [u,u1;v,v1]=prod_j([0,1;1,q_j]).
129 : * Loop invariants:
130 : * (a) (sign)*[-v,v1]*x = [d,d1] (mod m) (componentwise)
131 : * (sign initially +1, changes with each Euclidean step)
132 : * so [a,b] will be obtained in the form [-+d,v] or [+-d1,v1];
133 : * this congruence is a consequence of
134 : *
135 : * (b) [x,m]~ = [u,u1;v,v1]*[d1,d]~,
136 : * where u,u1 is the usual numerator sequence starting with 1,0
137 : * instead of 0,1 (just multiply the eqn on the left by the inverse
138 : * matrix, which is det*[v1,-u1;-v,u], where "det" is the same as the
139 : * "(sign)" in (a)). From m = v*d1 + v1*d and
140 : *
141 : * (c) d > d1 >= 0, 0 <= v < v1,
142 : * we have d >= m/(2*v1), so while v1 remains smaller than m/(2*amax),
143 : * the pair [-(sign)*d,v] satisfies (1) but violates (2) (d > amax).
144 : * Conversely, v1 > bmax indicates that no further solutions will be
145 : * forthcoming; [-(sign)*d,v] will be the last, and first, candidate.
146 : * Thus there's at most one point in the chain division where a solution
147 : * can live: v < bmax, v1 >= m/(2*amax) > bmax, and this is acceptable
148 : * iff in fact d <= amax (e.g. m=221, x=34 or 35, amax=bmax=10 fail on
149 : * this count while x=32,33,36,37 succeed). However, a division may leave
150 : * a zero residue before we ever reach this point (consider m=210, x=35,
151 : * amax=bmax=10), and our caller may find that gcd(d,v) > 1 (Examples:
152 : * keep m=210 and consider any of x=29,31,32,33,34,36,37,38,39,40,41).
153 : * Furthermore, at the start of the loop body we have in fact
154 : *
155 : * (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0,
156 : * (and are never done already).
157 : *
158 : * Main loop is similar to those of invmod() and bezout(), except for
159 : * having to determine appropriate vmax bounds, and checking termination
160 : * conditions. The signe(d1) condition is only for paranoia */
161 109622 : while (lgefint(d) > 3 && signe(d1))
162 : {
163 : /* determine vmax for lgcdii so as to ensure v won't overshoot.
164 : * If v+v1 > bmax, the next step would take v1 beyond the limit, so
165 : * since [+-d1,v1] is not a solution, we give up. Otherwise if v+v1
166 : * is way shorter than bmax, use vmax=ULONG_MAX. Otherwise, set vmax
167 : * to a crude lower approximation of bmax/(v+v1), or to 1, which will
168 : * allow the inner loop to do one step */
169 68509 : r = addii(v,v1);
170 68509 : if (cmpii(r,bmax) > 0) return gc_long(av, 0); /* done, not found */
171 68016 : vmax = get_vmax(r, lb, lbb);
172 : /* do a Lehmer-Jebelean round */
173 68016 : lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, vmax);
174 68016 : if (lhmres) /* check progress */
175 : { /* apply matrix */
176 68013 : if (lhmres == 1 || lhmres == -1)
177 : {
178 767 : s = -s;
179 767 : if (xv1 == 1)
180 : { /* re-use v+v1 computed above */
181 2 : v = v1; v1 = r;
182 2 : r = subii(d,d1); d = d1; d1 = r;
183 : }
184 : else
185 : {
186 765 : r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
187 765 : r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
188 : }
189 : }
190 : else
191 : {
192 67246 : r = subii(muliu(d,xu), muliu(d1,xv));
193 67246 : d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
194 67246 : r = addii(muliu(v,xu), muliu(v1,xv));
195 67246 : v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
196 67246 : if (lhmres&1) { togglesign(d); s = -s; } else togglesign(d1);
197 : }
198 : /* check whether we're done. Assert v <= bmax here. Examine v1:
199 : * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
200 : * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed*/
201 68013 : if (cmpii(v1,bmax) > 0)
202 : {
203 31791 : set_avma(av);
204 31791 : if (cmpii(d,amax) > 0) return 0; /* done, not found */
205 : /* done, found */
206 16857 : *a = icopy(d); setsigne(*a,-s);
207 16857 : *b = icopy(v); return 1;
208 : }
209 36222 : if (cmpii(d1,amax) <= 0)
210 : { /* done, found */
211 13121 : set_avma(av);
212 13121 : if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
213 13121 : *b = icopy(v1); return 1;
214 : }
215 : } /* lhmres != 0 */
216 :
217 23104 : if (lhmres <= 0 && signe(d1))
218 : {
219 5 : q = dvmdii(d,d1,&r);
220 5 : d = d1; d1 = r;
221 5 : r = addii(v, mulii(q,v1));
222 5 : v = v1; v1 = r;
223 5 : s = -s;
224 : /* check whether we are done now. Since we weren't before the div, it
225 : * suffices to examine v1 and d1 -- the new d (former d1) cannot cut it */
226 5 : if (cmpii(v1,bmax) > 0) return gc_long(av, 0); /* done, not found */
227 5 : if (cmpii(d1,amax) <= 0) /* done, found */
228 : {
229 0 : set_avma(av);
230 0 : if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
231 0 : *b = icopy(v1); return 1;
232 : }
233 : }
234 :
235 23104 : if (gc_needed(av,1))
236 : {
237 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ratlift");
238 0 : gerepileall(av1, 4, &d, &d1, &v, &v1);
239 : }
240 : } /* end while */
241 :
242 : /* Postprocessing - final sprint. Since we usually underestimate vmax,
243 : * this function needs a loop here instead of a simple conditional.
244 : * Note we can only get here when amax fits into one word (which will
245 : * typically not be the case!). The condition is bogus -- d1 is never
246 : * zero at the start of the loop. There will be at most a few iterations,
247 : * so we don't bother collecting garbage */
248 46352 : while (signe(d1))
249 : {
250 : /* Assertions: lgefint(d)==lgefint(d1)==3.
251 : * Moreover, we aren't done already, or we would have returned by now.
252 : * Recompute vmax */
253 46352 : r = addii(v,v1);
254 46352 : if (cmpii(r,bmax) > 0) return gc_long(av, 0); /* done, not found */
255 41744 : vmax = get_vmax(r, lb, lbb);
256 : /* single-word "Lehmer", discarding the gcd or whatever it returns */
257 41744 : (void)rgcduu((ulong)*int_MSW(d), (ulong)*int_MSW(d1), vmax, &xu, &xu1, &xv, &xv1, &s0);
258 41746 : if (xv1 == 1) /* avoid multiplications */
259 : { /* re-use r = v+v1 computed above */
260 0 : v = v1; v1 = r;
261 0 : r = subii(d,d1); d = d1; d1 = r;
262 0 : s = -s;
263 : }
264 41746 : else if (xu == 0) /* and xv==1, xu1==1, xv1 > 1 */
265 : {
266 7547 : r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
267 7547 : r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
268 7546 : s = -s;
269 : }
270 : else
271 : {
272 34199 : r = subii(muliu(d,xu), muliu(d1,xv));
273 34199 : d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
274 34199 : r = addii(muliu(v,xu), muliu(v1,xv));
275 34197 : v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
276 34199 : if (s0 < 0) { togglesign(d); s = -s; } else togglesign(d1);
277 : }
278 : /* check whether we're done, as above. Assert v <= bmax.
279 : * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
280 : * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed.
281 : */
282 41745 : if (cmpii(v1,bmax) > 0)
283 : {
284 31702 : set_avma(av);
285 31702 : if (cmpii(d,amax) > 0) return 0; /* done, not found */
286 : /* done, found */
287 21634 : *a = icopy(d); setsigne(*a,-s);
288 21634 : *b = icopy(v); return 1;
289 : }
290 10042 : if (cmpii(d1,amax) <= 0)
291 : { /* done, found */
292 4804 : set_avma(av);
293 4805 : if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
294 4805 : *b = icopy(v1); return 1;
295 : }
296 : } /* while */
297 :
298 : /* We have run into d1 == 0 before returning. This cannot happen */
299 0 : pari_err_BUG("ratlift failed to catch d1 == 0");
300 : return 0; /* LCOV_EXCL_LINE */
301 : }
|
