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 14853612 : invmod(GEN a, GEN b, GEN *res)
24 : {
25 14853612 : if (!signe(b)) { *res=absi(a); return 0; }
26 14853612 : if (NLIMBS(b) < INVMOD_GMP_LIMIT)
27 12587675 : return invmod_pari(a,b,res);
28 : { /* General case: use gcdext(a+b, b) since mpn_gcdext require S1>=S2 */
29 2265937 : pari_sp av = avma;
30 : GEN ca, cb, u, d;
31 2265937 : long l, su, sa = signe(a), lb,lna;
32 : mp_size_t lu;
33 : GEN na;
34 2265937 : if (!sa) { set_avma(av); *res = absi(b); return 0; }
35 2265933 : if (signe(b) < 0) b = negi(b);
36 2265933 : if (abscmpii(a, b) < 0)
37 2259063 : na = sa > 0? addii(a, b): subii(a, b);
38 : else
39 6940 : 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 for GC at the end. */
43 2265996 : lb = lgefint(b); lna = lgefint(na);
44 2265996 : ca = icopy_ef(na,lna+1);
45 2265993 : cb = icopy_ef( b,lb+1);
46 : /* must create u first else final icopy could fail. */
47 2265993 : u = cgeti(lna+1);
48 2265995 : d = cgeti(lna+1);
49 : /* na >= b */
50 2265992 : l = mpn_gcdext(LIMBS(d), LIMBS(u), &lu, LIMBS(ca), NLIMBS(ca), LIMBS(cb), NLIMBS(cb));
51 2266012 : d[1] = evalsigne(1)|evallgefint(l+2);
52 2266012 : if (!is_pm1(d)) {set_avma(av); *res=icopy(d); return 0;}
53 2265841 : su = lu?((sa ^ lu) < 0)? -1: 1: 0;
54 2265841 : u[1] = evalsigne(su) | evallgefint(labs(lu)+2);
55 2265841 : if (su < 0) u = addii(u, b);
56 2265838 : 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 91551488 : 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 91551488 : s = abscmpii(a,b);
82 91529135 : if (s < 0) { swap(a,b); pswap(pu,pv); }
83 : /* now |a| >= |b| */
84 :
85 91529135 : sa = signe(a); sb = signe(b);
86 91529135 : if (!sb)
87 : {
88 1408957 : if (pv) *pv = gen_0;
89 1408957 : switch(sa)
90 : {
91 4 : case 0: if (pu) *pu = gen_0; return gen_0;
92 1405180 : case 1: if (pu) *pu = gen_1; return icopy(a);
93 3772 : case -1: if (pu) *pu = gen_m1; return negi(a);
94 : }
95 : }
96 90120179 : if (s == 0) /* |a| == |b| != 0 */
97 : {
98 8892248 : if (pu) *pu = gen_0;
99 8892248 : if (sb > 0)
100 8370228 : { if (pv) *pv = gen_1; return icopy(b); }
101 : else
102 522020 : { if (pv) *pv = gen_m1; return negi(b); }
103 : }
104 : /* now |a| > |b| > 0 */
105 :
106 81227931 : if (lgefint(a) == 3) /* single-word affair */
107 : {
108 78081759 : g = xxgcduu((ulong)a[2], (ulong)b[2], 0, &xu, &xu1, &xv, &xv1, &s);
109 78813842 : sa = s > 0 ? sa : -sa;
110 78813842 : sb = s > 0 ? -sb : sb;
111 78813842 : if (pu)
112 : {
113 40629013 : if (xu == 0) *pu = gen_0; /* can happen when b divides a */
114 15735706 : else if (xu == 1) *pu = sa < 0 ? gen_m1 : gen_1;
115 9714127 : else if (xu == 2) *pu = sa < 0 ? gen_m2 : gen_2;
116 : else
117 : {
118 8621979 : *pu = cgeti(3);
119 8620291 : (*pu)[1] = evalsigne(sa)|evallgefint(3);
120 8620291 : (*pu)[2] = xu;
121 : }
122 : }
123 78812154 : if (pv)
124 : {
125 70758914 : if (xv == 1) *pv = sb < 0 ? gen_m1 : gen_1;
126 32224239 : else if (xv == 2) *pv = sb < 0 ? gen_m2 : gen_2;
127 : else
128 : {
129 29989191 : *pv = cgeti(3);
130 29873752 : (*pv)[1] = evalsigne(sb)|evallgefint(3);
131 29873752 : (*pv)[2] = xv;
132 : }
133 : }
134 78696715 : if (g == 1) return gen_1;
135 25932995 : else if (g == 2) return gen_2;
136 18035704 : else return utoipos(g);
137 : }
138 : else
139 : { /* general case */
140 3146172 : 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 for GC at the end.
144 : * NOTE: we must put u before d else the final icopy could fail. */
145 3146172 : GEN ca = icopy_ef(a,lgefint(a)+1);
146 3284820 : GEN cb = icopy_ef(b,lgefint(b)+1);
147 3284816 : GEN u = cgeti(lgefint(a)+1), v = NULL;
148 3284805 : GEN d = cgeti(lgefint(a)+1);
149 : long su,l;
150 : mp_size_t lu;
151 3284804 : l = mpn_gcdext(LIMBS(d), LIMBS(u), &lu, LIMBS(ca), NLIMBS(ca), LIMBS(cb), NLIMBS(cb));
152 3284823 : if (lu<=0)
153 : {
154 2770412 : if (lu==0) su=0;
155 399780 : else {su=-1;lu=-lu;}
156 : }
157 : else
158 514411 : su=1;
159 3284823 : if (sa<0) su=-su;
160 3284823 : d[1] = evalsigne(1)|evallgefint(l+2);
161 3284823 : u[1] = evalsigne(su)|evallgefint(lu+2);
162 3284823 : if (pv) v=diviiexact(subii(d,mulii(u,a)),b);
163 3284817 : set_avma(av);
164 3284818 : if (pu) *pu=icopy(u);
165 3284813 : if (pv) *pv=icopy(v);
166 3284812 : return icopy(d);
167 : }
168 : }
|
