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 59691672 : 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 59691672 : s = abscmpii(a,b);
42 59691672 : if (s < 0)
43 : {
44 37837596 : t=b; b=a; a=t;
45 37837596 : pt=pu; pu=pv; pv=pt;
46 : }
47 : /* now |a| >= |b| */
48 :
49 59691672 : sa = signe(a); sb = signe(b);
50 59691672 : if (!sb)
51 : {
52 1011960 : if (pv) *pv = gen_0;
53 1011960 : switch(sa)
54 : {
55 3 : case 0: if (pu) *pu = gen_0; return gen_0;
56 1009128 : case 1: if (pu) *pu = gen_1; return icopy(a);
57 2829 : case -1: if (pu) *pu = gen_m1; return(negi(a));
58 : }
59 : }
60 58679712 : if (s == 0) /* |a| == |b| != 0 */
61 : {
62 6582576 : if (pu) *pu = gen_0;
63 6582576 : if (sb > 0)
64 6191214 : { if (pv) *pv = gen_1; return icopy(b); }
65 : else
66 391362 : { if (pv) *pv = gen_m1; return(negi(b)); }
67 : }
68 : /* now |a| > |b| > 0 */
69 :
70 52097136 : if (lgefint(a) == 3) /* single-word affair */
71 : {
72 50114613 : g = xxgcduu(uel(a,2), uel(b,2), 0, &xu, &xu1, &xv, &xv1, &s);
73 50114613 : sa = s > 0 ? sa : -sa;
74 50114613 : sb = s > 0 ? -sb : sb;
75 50114613 : if (pu)
76 : {
77 27405714 : if (xu == 0) *pu = gen_0; /* can happen when b divides a */
78 9355401 : else if (xu == 1) *pu = sa < 0 ? gen_m1 : gen_1;
79 5295024 : else if (xu == 2) *pu = sa < 0 ? gen_m2 : gen_2;
80 : else
81 : {
82 4718001 : *pu = cgeti(3);
83 4718001 : (*pu)[1] = evalsigne(sa)|evallgefint(3);
84 4718001 : (*pu)[2] = xu;
85 : }
86 : }
87 50114613 : if (pv)
88 : {
89 45906006 : if (xv == 1) *pv = sb < 0 ? gen_m1 : gen_1;
90 18370068 : else if (xv == 2) *pv = sb < 0 ? gen_m2 : gen_2;
91 : else
92 : {
93 16792947 : *pv = cgeti(3);
94 16792947 : (*pv)[1] = evalsigne(sb)|evallgefint(3);
95 16792947 : (*pv)[2] = xv;
96 : }
97 : }
98 50114613 : if (g == 1) return gen_1;
99 17541129 : else if (g == 2) return gen_2;
100 11972787 : else return utoipos(g);
101 : }
102 :
103 : /* general case */
104 1982523 : av = avma;
105 1982523 : (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 1982523 : if (lgefint(a) > lgefint(b))
110 : {
111 1497603 : d = absi_shallow(b);
112 1497603 : q = dvmdii(absi_shallow(a), d, &d1);
113 1497603 : if (!signe(d1)) /* a == qb */
114 : {
115 1273812 : set_avma(av);
116 1273812 : if (pu) *pu = gen_0;
117 1273812 : if (pv) *pv = sb < 0 ? gen_m1 : gen_1;
118 1273812 : return icopy(d);
119 : }
120 : else
121 : {
122 223791 : u = gen_0;
123 223791 : u1 = v = gen_1;
124 223791 : v1 = negi(q);
125 : }
126 : /* if this results in lgefint(d) == 3, will fall past main loop */
127 : }
128 : else
129 : {
130 484920 : d = absi_shallow(a);
131 484920 : d1 = absi_shallow(b);
132 484920 : u = v1 = gen_1; u1 = v = gen_0;
133 : }
134 708711 : av1 = avma;
135 :
136 : /* main loop is almost identical to that of invmod() */
137 2238156 : while (lgefint(d) > 3 && signe(d1))
138 : {
139 1529445 : lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, ULONG_MAX);
140 1529445 : if (lhmres != 0) /* check progress */
141 : { /* apply matrix */
142 1153026 : if ((lhmres == 1) || (lhmres == -1))
143 : {
144 66000 : if (xv1 == 1)
145 : {
146 31761 : r = subii(d,d1); d=d1; d1=r;
147 31761 : a = subii(u,u1); u=u1; u1=a;
148 31761 : a = subii(v,v1); v=v1; v1=a;
149 : }
150 : else
151 : {
152 34239 : r = subii(d, mului(xv1,d1)); d=d1; d1=r;
153 34239 : a = subii(u, mului(xv1,u1)); u=u1; u1=a;
154 34239 : a = subii(v, mului(xv1,v1)); v=v1; v1=a;
155 : }
156 : }
157 : else
158 : {
159 1087026 : r = subii(muliu(d,xu), muliu(d1,xv));
160 1087026 : d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
161 1087026 : a = subii(muliu(u,xu), muliu(u1,xv));
162 1087026 : u1 = subii(muliu(u,xu1), muliu(u1,xv1)); u = a;
163 1087026 : a = subii(muliu(v,xu), muliu(v1,xv));
164 1087026 : v1 = subii(muliu(v,xu1), muliu(v1,xv1)); v = a;
165 1087026 : if (lhmres&1) { togglesign(d); togglesign(u); togglesign(v); }
166 552459 : else { togglesign(d1); togglesign(u1); togglesign(v1); }
167 : }
168 : }
169 1529445 : if (lhmres <= 0 && signe(d1))
170 : {
171 405303 : q = dvmdii(d,d1,&r);
172 405303 : a = subii(u,mulii(q,u1));
173 405303 : u=u1; u1=a;
174 405303 : a = subii(v,mulii(q,v1));
175 405303 : v=v1; v1=a;
176 405303 : d=d1; d1=r;
177 : }
178 1529445 : if (gc_needed(av,1))
179 : {
180 0 : if(DEBUGMEM>1) pari_warn(warnmem,"bezout");
181 0 : (void)gc_all(av1,6, &d,&d1,&u,&u1,&v,&v1);
182 : }
183 : } /* end while */
184 :
185 : /* Postprocessing - final sprint */
186 708711 : 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 330570 : g = xxgcduu(uel(d,2), uel(d1,2), 0, &xu, &xu1, &xv, &xv1, &s);
192 330570 : u = subii(muliu(u,xu), muliu(u1, xv));
193 330570 : v = subii(muliu(v,xu), muliu(v1, xv));
194 330570 : if (s < 0) { sa = -sa; sb = -sb; }
195 330570 : set_avma(av);
196 330570 : if (pu) *pu = sa < 0 ? negi(u) : icopy(u);
197 330570 : if (pv) *pv = sb < 0 ? negi(v) : icopy(v);
198 330570 : if (g == 1) return gen_1;
199 261054 : else if (g == 2) return gen_2;
200 244140 : 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 378141 : set_avma(av);
206 378141 : if (pu) *pu = sa < 0 ? negi(u) : icopy(u);
207 378141 : if (pv) *pv = sb < 0 ? negi(v) : icopy(v);
208 378141 : return icopy(d);
209 : }
210 :
211 : static GEN
212 747 : addmulmul(GEN u, GEN v, GEN x, GEN y)
213 747 : { return addii(mulii(u, x),mulii(v, y)); }
214 :
215 : static GEN
216 408 : bezout_halfgcd(GEN x, GEN y, GEN *ptu, GEN *ptv)
217 : {
218 408 : pari_sp av=avma;
219 408 : GEN u, v, R = matid2();
220 942 : while (lgefint(y)-2 >= EXTGCD_HALFGCD_LIMIT)
221 : {
222 534 : GEN M = HGCD0(x,y);
223 534 : R = ZM2_mul(R, gel(M, 1));
224 534 : x = gel(M,2); y = gel(M,3);
225 534 : if (signe(y) && expi(y)<magic_threshold(x))
226 : {
227 297 : GEN r, q = dvmdii(x, y, &r);
228 297 : x = y; y = r;
229 297 : R = mulq(R, q);
230 : }
231 534 : if (gc_needed(av, 1))
232 0 : (void)gc_all(av,3,&x,&y,&R);
233 : }
234 408 : y = bezout_lehmer(x,y,&u,&v);
235 408 : R = ZM_inv2(R);
236 408 : if (ptu) *ptu = addmulmul(u,v,gcoeff(R,1,1),gcoeff(R,2,1));
237 408 : if (ptv) *ptv = addmulmul(u,v,gcoeff(R,1,2),gcoeff(R,2,2));
238 408 : return y;
239 : }
240 :
241 : GEN
242 59691672 : bezout(GEN x, GEN y, GEN *ptu, GEN *ptv)
243 : {
244 59691672 : pari_sp av = avma;
245 : GEN d;
246 59691672 : if (lgefint(y)-2 >= EXTGCD_HALFGCD_LIMIT)
247 408 : d = bezout_halfgcd(x, y, ptu, ptv);
248 : else
249 59691264 : d = bezout_lehmer(x, y, ptu, ptv);
250 59691672 : if (ptv) return gc_all(av,ptu?3:2, &d, ptv, ptu);
251 29556978 : return gc_all(av, ptu?2:1, &d, ptu);
252 : }
|
