Line data Source code
1 : #line 2 "../src/kernel/gmp/gcdext.c"
2 : /* Copyright (C) 2000-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 : * invmod(a,b,res)
18 : *==================================
19 : * If a is invertible, return 1, and set res = a^{ -1 }
20 : * Otherwise, return 0, and set res = gcd(a,b)
21 : */
22 : int
23 22145191 : invmod(GEN a, GEN b, GEN *res)
24 : {
25 22145191 : if (!signe(b)) { *res=absi(a); return 0; }
26 22145191 : if (NLIMBS(b) < INVMOD_GMP_LIMIT)
27 19959765 : return invmod_pari(a,b,res);
28 : { /* General case: use gcdext(a+b, b) since mpn_gcdext require S1>=S2 */
29 2185426 : pari_sp av = avma;
30 : GEN ca, cb, u, d;
31 2185426 : long l, su, sa = signe(a), lb,lna;
32 : mp_size_t lu;
33 : GEN na;
34 2185426 : if (!sa) { set_avma(av); *res = absi(b); return 0; }
35 2185422 : if (signe(b) < 0) b = negi(b);
36 2185422 : if (abscmpii(a, b) < 0)
37 2178670 : na = sa > 0? addii(a, b): subii(a, b);
38 : else
39 6798 : na = a;
40 : /* Copy serves two purposes:
41 : * 1) mpn_gcdext destroys its input and needs an extra limb
42 : * 2) allows us to use icopy instead of gerepile later. */
43 2185466 : lb = lgefint(b); lna = lgefint(na);
44 2185466 : ca = icopy_ef(na,lna+1);
45 2185467 : cb = icopy_ef( b,lb+1);
46 : /* must create u first else final icopy could fail. */
47 2185468 : u = cgeti(lna+1);
48 2185467 : d = cgeti(lna+1);
49 : /* na >= b */
50 2185464 : l = mpn_gcdext(LIMBS(d), LIMBS(u), &lu, LIMBS(ca), NLIMBS(ca), LIMBS(cb), NLIMBS(cb));
51 2185468 : d[1] = evalsigne(1)|evallgefint(l+2);
52 2185468 : if (!is_pm1(d)) {set_avma(av); *res=icopy(d); return 0;}
53 2185321 : su = lu?((sa ^ lu) < 0)? -1: 1: 0;
54 2185321 : u[1] = evalsigne(su) | evallgefint(labs(lu)+2);
55 2185321 : if (su < 0) u = addii(u, b);
56 2185319 : set_avma(av); *res=icopy(u); return 1;
57 : }
58 : }
59 :
60 : /*==================================
61 : * bezout(a,b,pu,pv)
62 : *==================================
63 : * Return g = gcd(a,b) >= 0, and assign GENs u,v through pointers pu,pv
64 : * such that g = u*a + v*b.
65 : * Special cases:
66 : * a == b == 0 ==> pick u=1, v=0
67 : * a != 0 == b ==> keep v=0
68 : * a == 0 != b ==> keep u=0
69 : * |a| == |b| != 0 ==> keep u=0, set v=+-1
70 : * Assignments through pu,pv will be suppressed when the corresponding
71 : * pointer is NULL (but the computations will happen nonetheless).
72 : */
73 :
74 : GEN
75 77634373 : bezout(GEN a, GEN b, GEN *pu, GEN *pv)
76 : {
77 : long s, sa, sb;
78 : ulong g;
79 : ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
80 :
81 77634373 : s = abscmpii(a,b);
82 77631876 : if (s < 0) { swap(a,b); pswap(pu,pv); }
83 : /* now |a| >= |b| */
84 :
85 77631876 : sa = signe(a); sb = signe(b);
86 77631876 : if (!sb)
87 : {
88 1340828 : if (pv) *pv = gen_0;
89 1340828 : switch(sa)
90 : {
91 4 : case 0: if (pu) *pu = gen_0; return gen_0;
92 1337068 : case 1: if (pu) *pu = gen_1; return icopy(a);
93 3756 : case -1: if (pu) *pu = gen_m1; return negi(a);
94 : }
95 : }
96 76291048 : if (s == 0) /* |a| == |b| != 0 */
97 : {
98 9216011 : if (pu) *pu = gen_0;
99 9216011 : if (sb > 0)
100 8711326 : { if (pv) *pv = gen_1; return icopy(b); }
101 : else
102 504685 : { if (pv) *pv = gen_m1; return negi(b); }
103 : }
104 : /* now |a| > |b| > 0 */
105 :
106 67075037 : if (lgefint(a) == 3) /* single-word affair */
107 : {
108 63905816 : g = xxgcduu((ulong)a[2], (ulong)b[2], 0, &xu, &xu1, &xv, &xv1, &s);
109 64150382 : sa = s > 0 ? sa : -sa;
110 64150382 : sb = s > 0 ? -sb : sb;
111 64150382 : if (pu)
112 : {
113 34646182 : if (xu == 0) *pu = gen_0; /* can happen when b divides a */
114 12024787 : else if (xu == 1) *pu = sa < 0 ? gen_m1 : gen_1;
115 6814871 : else if (xu == 2) *pu = sa < 0 ? gen_m2 : gen_2;
116 : else
117 : {
118 6077870 : *pu = cgeti(3);
119 6079240 : (*pu)[1] = evalsigne(sa)|evallgefint(3);
120 6079240 : (*pu)[2] = xu;
121 : }
122 : }
123 64151752 : if (pv)
124 : {
125 58802181 : if (xv == 1) *pv = sb < 0 ? gen_m1 : gen_1;
126 23987819 : else if (xv == 2) *pv = sb < 0 ? gen_m2 : gen_2;
127 : else
128 : {
129 21860408 : *pv = cgeti(3);
130 21824958 : (*pv)[1] = evalsigne(sb)|evallgefint(3);
131 21824958 : (*pv)[2] = xv;
132 : }
133 : }
134 64116302 : if (g == 1) return gen_1;
135 22348229 : else if (g == 2) return gen_2;
136 15436335 : else return utoipos(g);
137 : }
138 : else
139 : { /* general case */
140 3169221 : pari_sp av = avma;
141 : /*Copy serves two purposes:
142 : * 1) mpn_gcdext destroys its input and needs an extra limb
143 : * 2) allows us to use icopy instead of gerepile later.
144 : * NOTE: we must put u before d else the final icopy could fail. */
145 3169221 : GEN ca = icopy_ef(a,lgefint(a)+1);
146 3235507 : GEN cb = icopy_ef(b,lgefint(b)+1);
147 3235505 : GEN u = cgeti(lgefint(a)+1), v = NULL;
148 3235495 : GEN d = cgeti(lgefint(a)+1);
149 : long su,l;
150 : mp_size_t lu;
151 3235493 : l = mpn_gcdext(LIMBS(d), LIMBS(u), &lu, LIMBS(ca), NLIMBS(ca), LIMBS(cb), NLIMBS(cb));
152 3235512 : if (lu<=0)
153 : {
154 2724642 : if (lu==0) su=0;
155 396635 : else {su=-1;lu=-lu;}
156 : }
157 : else
158 510870 : su=1;
159 3235512 : if (sa<0) su=-su;
160 3235512 : d[1] = evalsigne(1)|evallgefint(l+2);
161 3235512 : u[1] = evalsigne(su)|evallgefint(lu+2);
162 3235512 : if (pv) v=diviiexact(subii(d,mulii(u,a)),b);
163 3235503 : set_avma(av);
164 3235502 : if (pu) *pu=icopy(u);
165 3235501 : if (pv) *pv=icopy(v);
166 3235499 : return icopy(d);
167 : }
168 : }
|
