Line data Source code
1 : #line 2 "../src/kernel/gmp/gcd.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 : /* assume y > x > 0. return y mod x */
17 : static ulong
18 39101529 : resiu(GEN y, ulong x)
19 : {
20 39101529 : return mpn_mod_1(LIMBS(y), NLIMBS(y), x);
21 : }
22 :
23 : GEN
24 866236736 : gcdii(GEN a, GEN b)
25 : {
26 : long v, w;
27 : pari_sp av;
28 : GEN t;
29 :
30 866236736 : switch (abscmpii(a,b))
31 : {
32 238228500 : case 0: return absi(a);
33 307959889 : case -1: swap(a,b);
34 : }
35 627959287 : if (!signe(b)) return absi(a);
36 : /* here |a|>|b|>0. Try single precision first */
37 374512252 : if (lgefint(a)==3)
38 311231692 : return igcduu((ulong)a[2], (ulong)b[2]);
39 63280560 : if (lgefint(b)==3)
40 : {
41 39101600 : ulong u = resiu(a,(ulong)b[2]);
42 39102856 : if (!u) return absi(b);
43 23845280 : return igcduu((ulong)b[2], u);
44 : }
45 : /* larger than gcd: GC using set_avma(av) (erasing t) is valid */
46 24178960 : av = avma; (void)new_chunk(lgefint(b)+1); /* HACK */
47 24233841 : t = remii(a,b);
48 24240442 : if (!signe(t)) { set_avma(av); return absi(b); }
49 :
50 16313367 : a = b; b = t;
51 16313367 : v = vali(a); a = shifti(a,-v); setabssign(a);
52 16312530 : w = vali(b); b = shifti(b,-w); setabssign(b);
53 16312527 : if (w < v) v = w;
54 16312527 : switch(abscmpii(a,b))
55 : {
56 305455 : case 0: set_avma(av); a=shifti(a,v); return a;
57 3450494 : case -1: swap(a,b);
58 : }
59 16007087 : if (is_pm1(b)) { set_avma(av); return int2n(v); }
60 : {
61 : /* Serves two purposes: 1) mpn_gcd destroys its input and needs an extra limb
62 : * 2) allows us to use icopy for GC. */
63 15878778 : GEN res= cgeti(lgefint(a)+1);
64 15878455 : GEN ca = icopy_ef(a,lgefint(a)+1);
65 15877940 : GEN cb = icopy_ef(b,lgefint(b)+1);
66 15878403 : long l = mpn_gcd(LIMBS(res), LIMBS(ca), NLIMBS(ca), LIMBS(cb), NLIMBS(cb));
67 15879492 : res[1] = evalsigne(1)|evallgefint(l+2);
68 15879492 : set_avma(av); return shifti(res,v);
69 : }
70 : }
|
