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 22674244 : resiu(GEN y, ulong x)
19 : {
20 22674244 : return mpn_mod_1(LIMBS(y), NLIMBS(y), x);
21 : }
22 :
23 : GEN
24 326556796 : gcdii(GEN a, GEN b)
25 : {
26 : long v, w;
27 : pari_sp av;
28 : GEN t;
29 :
30 326556796 : switch (abscmpii(a,b))
31 : {
32 61688280 : case 0: return absi(a);
33 133795146 : case -1: swap(a,b);
34 : }
35 264906309 : if (!signe(b)) return absi(a);
36 : /* here |a|>|b|>0. Try single precision first */
37 173288290 : if (lgefint(a)==3)
38 132482413 : return igcduu((ulong)a[2], (ulong)b[2]);
39 40805877 : if (lgefint(b)==3)
40 : {
41 22674254 : ulong u = resiu(a,(ulong)b[2]);
42 22674267 : if (!u) return absi(b);
43 11314520 : return igcduu((ulong)b[2], u);
44 : }
45 : /* larger than gcd: "set_avma(av)" gerepile (erasing t) is valid */
46 18131623 : av = avma; (void)new_chunk(lgefint(b)+1); /* HACK */
47 18134876 : t = remii(a,b);
48 18135791 : if (!signe(t)) { set_avma(av); return absi(b); }
49 :
50 11847107 : a = b; b = t;
51 11847107 : v = vali(a); a = shifti(a,-v); setabssign(a);
52 11846234 : w = vali(b); b = shifti(b,-w); setabssign(b);
53 11846140 : if (w < v) v = w;
54 11846140 : switch(abscmpii(a,b))
55 : {
56 239973 : case 0: set_avma(av); a=shifti(a,v); return a;
57 2411604 : case -1: swap(a,b);
58 : }
59 11606385 : if (is_pm1(b)) { set_avma(av); return int2n(v); }
60 : {
61 : /* general case */
62 : /*This serve two purposes: 1) mpn_gcd destroy its input and need an extra
63 : * limb 2) this allows us to use icopy instead of gerepile later. NOTE: we
64 : * must put u before d else the final icopy could fail.
65 : */
66 11490609 : GEN res= cgeti(lgefint(a)+1);
67 11489696 : GEN ca = icopy_ef(a,lgefint(a)+1);
68 11489084 : GEN cb = icopy_ef(b,lgefint(b)+1);
69 11489484 : long l = mpn_gcd(LIMBS(res), LIMBS(ca), NLIMBS(ca), LIMBS(cb), NLIMBS(cb));
70 11491256 : res[1] = evalsigne(1)|evallgefint(l+2);
71 11491256 : set_avma(av);
72 11491032 : return shifti(res,v);
73 : }
74 : }
|
