Line data Source code
1 : #line 2 "../src/kernel/none/gcdext.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 : * bezout(a,b,pu,pv)
18 : *==================================
19 : * Return g = gcd(a,b) >= 0, and assign GENs u,v through pointers pu,pv
20 : * such that g = u*a + v*b.
21 : * Special cases:
22 : * a == b == 0 ==> pick u=1, v=0
23 : * a != 0 == b ==> keep v=0
24 : * a == 0 != b ==> keep u=0
25 : * |a| == |b| != 0 ==> keep u=0, set v=+-1
26 : * Assignments through pu,pv will be suppressed when the corresponding
27 : * pointer is NULL (but the computations will happen nonetheless).
28 : */
29 :
30 : static GEN
31 57845820 : bezout_lehmer(GEN a, GEN b, GEN *pu, GEN *pv)
32 : {
33 : GEN t,u,u1,v,v1,d,d1,q,r;
34 : GEN *pt;
35 : pari_sp av, av1;
36 : long s, sa, sb;
37 : ulong g;
38 : ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
39 : int lhmres; /* Lehmer stage return value */
40 :
41 57845820 : s = abscmpii(a,b);
42 57845820 : if (s < 0)
43 : {
44 36416943 : t=b; b=a; a=t;
45 36416943 : pt=pu; pu=pv; pv=pt;
46 : }
47 : /* now |a| >= |b| */
48 :
49 57845820 : sa = signe(a); sb = signe(b);
50 57845820 : if (!sb)
51 : {
52 985722 : if (pv) *pv = gen_0;
53 985722 : switch(sa)
54 : {
55 3 : case 0: if (pu) *pu = gen_0; return gen_0;
56 983724 : case 1: if (pu) *pu = gen_1; return icopy(a);
57 1995 : case -1: if (pu) *pu = gen_m1; return(negi(a));
58 : }
59 : }
60 56860098 : if (s == 0) /* |a| == |b| != 0 */
61 : {
62 6842808 : if (pu) *pu = gen_0;
63 6842808 : if (sb > 0)
64 6469863 : { if (pv) *pv = gen_1; return icopy(b); }
65 : else
66 372945 : { if (pv) *pv = gen_m1; return(negi(b)); }
67 : }
68 : /* now |a| > |b| > 0 */
69 :
70 50017290 : if (lgefint(a) == 3) /* single-word affair */
71 : {
72 48103659 : g = xxgcduu(uel(a,2), uel(b,2), 0, &xu, &xu1, &xv, &xv1, &s);
73 48103659 : sa = s > 0 ? sa : -sa;
74 48103659 : sb = s > 0 ? -sb : sb;
75 48103659 : if (pu)
76 : {
77 25931601 : if (xu == 0) *pu = gen_0; /* can happen when b divides a */
78 9024567 : else if (xu == 1) *pu = sa < 0 ? gen_m1 : gen_1;
79 5099460 : else if (xu == 2) *pu = sa < 0 ? gen_m2 : gen_2;
80 : else
81 : {
82 4551504 : *pu = cgeti(3);
83 4551504 : (*pu)[1] = evalsigne(sa)|evallgefint(3);
84 4551504 : (*pu)[2] = xu;
85 : }
86 : }
87 48103659 : if (pv)
88 : {
89 44067969 : if (xv == 1) *pv = sb < 0 ? gen_m1 : gen_1;
90 18047160 : else if (xv == 2) *pv = sb < 0 ? gen_m2 : gen_2;
91 : else
92 : {
93 16459320 : *pv = cgeti(3);
94 16459320 : (*pv)[1] = evalsigne(sb)|evallgefint(3);
95 16459320 : (*pv)[2] = xv;
96 : }
97 : }
98 48103659 : if (g == 1) return gen_1;
99 17146050 : else if (g == 2) return gen_2;
100 11904756 : else return utoipos(g);
101 : }
102 :
103 : /* general case */
104 1913631 : av = avma;
105 1913631 : (void)new_chunk(lgefint(b) + (lgefint(a)<<1)); /* room for u,v,gcd */
106 : /* if a is significantly larger than b, calling lgcdii() is not the best
107 : * way to start -- reduce a mod b first
108 : */
109 1913631 : if (lgefint(a) > lgefint(b))
110 : {
111 1446738 : d = absi_shallow(b);
112 1446738 : q = dvmdii(absi_shallow(a), d, &d1);
113 1446738 : if (!signe(d1)) /* a == qb */
114 : {
115 1224951 : set_avma(av);
116 1224951 : if (pu) *pu = gen_0;
117 1224951 : if (pv) *pv = sb < 0 ? gen_m1 : gen_1;
118 1224951 : return icopy(d);
119 : }
120 : else
121 : {
122 221787 : u = gen_0;
123 221787 : u1 = v = gen_1;
124 221787 : v1 = negi(q);
125 : }
126 : /* if this results in lgefint(d) == 3, will fall past main loop */
127 : }
128 : else
129 : {
130 466893 : d = absi_shallow(a);
131 466893 : d1 = absi_shallow(b);
132 466893 : u = v1 = gen_1; u1 = v = gen_0;
133 : }
134 688680 : av1 = avma;
135 :
136 : /* main loop is almost identical to that of invmod() */
137 2195247 : while (lgefint(d) > 3 && signe(d1))
138 : {
139 1506567 : lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, ULONG_MAX);
140 1506567 : if (lhmres != 0) /* check progress */
141 : { /* apply matrix */
142 1142235 : if ((lhmres == 1) || (lhmres == -1))
143 : {
144 60705 : if (xv1 == 1)
145 : {
146 28218 : r = subii(d,d1); d=d1; d1=r;
147 28218 : a = subii(u,u1); u=u1; u1=a;
148 28218 : a = subii(v,v1); v=v1; v1=a;
149 : }
150 : else
151 : {
152 32487 : r = subii(d, mului(xv1,d1)); d=d1; d1=r;
153 32487 : a = subii(u, mului(xv1,u1)); u=u1; u1=a;
154 32487 : a = subii(v, mului(xv1,v1)); v=v1; v1=a;
155 : }
156 : }
157 : else
158 : {
159 1081530 : r = subii(muliu(d,xu), muliu(d1,xv));
160 1081530 : d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
161 1081530 : a = subii(muliu(u,xu), muliu(u1,xv));
162 1081530 : u1 = subii(muliu(u,xu1), muliu(u1,xv1)); u = a;
163 1081530 : a = subii(muliu(v,xu), muliu(v1,xv));
164 1081530 : v1 = subii(muliu(v,xu1), muliu(v1,xv1)); v = a;
165 1081530 : if (lhmres&1) { togglesign(d); togglesign(u); togglesign(v); }
166 548028 : else { togglesign(d1); togglesign(u1); togglesign(v1); }
167 : }
168 : }
169 1506567 : if (lhmres <= 0 && signe(d1))
170 : {
171 390213 : q = dvmdii(d,d1,&r);
172 390213 : a = subii(u,mulii(q,u1));
173 390213 : u=u1; u1=a;
174 390213 : a = subii(v,mulii(q,v1));
175 390213 : v=v1; v1=a;
176 390213 : d=d1; d1=r;
177 : }
178 1506567 : if (gc_needed(av,1))
179 : {
180 0 : if(DEBUGMEM>1) pari_warn(warnmem,"bezout");
181 0 : gerepileall(av1,6, &d,&d1,&u,&u1,&v,&v1);
182 : }
183 : } /* end while */
184 :
185 : /* Postprocessing - final sprint */
186 688680 : if (signe(d1))
187 : {
188 : /* Assertions: lgefint(d)==lgefint(d1)==3, and
189 : * gcd(d,d1) is nonzero and fits into one word
190 : */
191 324873 : g = xxgcduu(uel(d,2), uel(d1,2), 0, &xu, &xu1, &xv, &xv1, &s);
192 324873 : u = subii(muliu(u,xu), muliu(u1, xv));
193 324873 : v = subii(muliu(v,xu), muliu(v1, xv));
194 324873 : if (s < 0) { sa = -sa; sb = -sb; }
195 324873 : set_avma(av);
196 324873 : if (pu) *pu = sa < 0 ? negi(u) : icopy(u);
197 324873 : if (pv) *pv = sb < 0 ? negi(v) : icopy(v);
198 324873 : if (g == 1) return gen_1;
199 258483 : else if (g == 2) return gen_2;
200 242292 : else return utoipos(g);
201 : }
202 : /* get here when the final sprint was skipped (d1 was zero already).
203 : * Now the matrix is final, and d contains the gcd.
204 : */
205 363807 : set_avma(av);
206 363807 : if (pu) *pu = sa < 0 ? negi(u) : icopy(u);
207 363807 : if (pv) *pv = sb < 0 ? negi(v) : icopy(v);
208 363807 : return icopy(d);
209 : }
210 :
211 : static GEN
212 390 : addmulmul(GEN u, GEN v, GEN x, GEN y)
213 390 : { return addii(mulii(u, x),mulii(v, y)); }
214 :
215 : static GEN
216 222 : bezout_halfgcd(GEN x, GEN y, GEN *ptu, GEN *ptv)
217 : {
218 222 : pari_sp av=avma;
219 222 : GEN u, v, R = matid2();
220 564 : while (lgefint(y)-2 >= EXTGCD_HALFGCD_LIMIT)
221 : {
222 342 : GEN M = HGCD0(x,y);
223 342 : R = ZM2_mul(R, gel(M, 1));
224 342 : x = gel(M,2); y = gel(M,3);
225 342 : if (signe(y) && expi(y)<magic_threshold(x))
226 : {
227 138 : GEN r, q = dvmdii(x, y, &r);
228 138 : x = y; y = r;
229 138 : R = mulq(R, q);
230 : }
231 342 : if (gc_needed(av, 1))
232 0 : gerepileall(av,3,&x,&y,&R);
233 : }
234 222 : y = bezout_lehmer(x,y,&u,&v);
235 222 : R = ZM_inv2(R);
236 222 : if (ptu) *ptu = addmulmul(u,v,gcoeff(R,1,1),gcoeff(R,2,1));
237 222 : if (ptv) *ptv = addmulmul(u,v,gcoeff(R,1,2),gcoeff(R,2,2));
238 222 : return y;
239 : }
240 :
241 : GEN
242 57845820 : bezout(GEN x, GEN y, GEN *ptu, GEN *ptv)
243 : {
244 57845820 : pari_sp av = avma;
245 : GEN d;
246 57845820 : if (lgefint(y)-2 >= EXTGCD_HALFGCD_LIMIT)
247 222 : d = bezout_halfgcd(x, y, ptu, ptv);
248 : else
249 57845598 : d = bezout_lehmer(x, y, ptu, ptv);
250 57845820 : if (ptv) return gc_all(av,ptu?3:2, &d, ptv, ptu);
251 29465922 : return gc_all(av, ptu?2:1, &d, ptu);
252 : }
|
