Line data Source code
1 : #line 2 "../src/kernel/none/gcd.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 : /* assume y > x > 0. return y mod x */
17 : static ulong
18 18862842 : resiu(GEN y, ulong x)
19 : {
20 18862842 : long i, ly = lgefint(y);
21 18862842 : ulong xi = get_Fl_red(x);
22 : LOCAL_HIREMAINDER;
23 :
24 18862842 : hiremainder = 0;
25 131130408 : for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
26 18862842 : return hiremainder;
27 : }
28 :
29 : /* Assume x>y>0, both of them odd. return x-y if x=y mod 4, x+y otherwise */
30 : static void
31 561306918 : gcd_plus_minus(GEN x, GEN y, GEN res)
32 : {
33 561306918 : pari_sp av = avma;
34 561306918 : long lx = lgefint(x)-1;
35 561306918 : long ly = lgefint(y)-1, lt,m,i;
36 : GEN t;
37 :
38 561306918 : if ((x[lx]^y[ly]) & 3) /* x != y mod 4*/
39 281638029 : t = addiispec(x+2,y+2,lx-1,ly-1);
40 : else
41 279668889 : t = subiispec(x+2,y+2,lx-1,ly-1);
42 :
43 562033800 : lt = lgefint(t)-1; while (!t[lt]) lt--;
44 561306918 : m = vals(t[lt]); lt++;
45 561306918 : if (m == 0) /* 2^32 | t */
46 : {
47 27330 : for (i = 2; i < lt; i++) res[i] = t[i];
48 : }
49 561294909 : else if (t[2] >> m)
50 : {
51 529556286 : shift_right(res,t, 2,lt, 0,m);
52 : }
53 : else
54 : {
55 31738623 : lt--; t++;
56 31738623 : shift_right(res,t, 2,lt, t[1],m);
57 : }
58 561306918 : res[1] = evalsigne(1)|evallgefint(lt);
59 561306918 : set_avma(av);
60 561306918 : }
61 :
62 : /* uses modified right-shift binary algorithm now --GN 1998Jul23 */
63 : static GEN
64 619975284 : gcdii_basecase(GEN a, GEN b)
65 : {
66 : long v, w;
67 : pari_sp av;
68 : GEN t, p1;
69 :
70 619975284 : switch (abscmpii(a,b))
71 : {
72 177518046 : case 0: return absi(a);
73 213632895 : case -1: swap(a,b);
74 : }
75 442457238 : if (!signe(b)) return absi(a);
76 : /* here |a|>|b|>0. Try single precision first */
77 249772254 : if (lgefint(a)==3)
78 218717373 : return igcduu((ulong)a[2], (ulong)b[2]);
79 31054881 : if (lgefint(b)==3)
80 : {
81 18862842 : ulong u = resiu(a,(ulong)b[2]);
82 18862842 : if (!u) return absi(b);
83 11153931 : return igcduu((ulong)b[2], u);
84 : }
85 :
86 : /* larger than gcd: "set_avma(av)" gerepile (erasing t) is valid */
87 12192039 : av = avma; (void)new_chunk(lgefint(b)); /* HACK */
88 12192039 : t = remii(a,b);
89 12192039 : if (!signe(t)) { set_avma(av); return absi(b); }
90 :
91 8448672 : a = b; b = t;
92 8448672 : v = vali(a); a = shifti(a,-v); setabssign(a);
93 8448672 : w = vali(b); b = shifti(b,-w); setabssign(b);
94 8448672 : if (w < v) v = w;
95 8448672 : switch(abscmpii(a,b))
96 : {
97 176409 : case 0: set_avma(av); a = shifti(a,v); return a;
98 1650846 : case -1: swap(a,b);
99 : }
100 8272263 : if (is_pm1(b)) { set_avma(av); return int2n(v); }
101 :
102 : /* we have three consecutive memory locations: a,b,t.
103 : * All computations are done in place */
104 :
105 : /* a and b are odd now, and a>b>1 */
106 567993837 : while (lgefint(a) > 3)
107 : {
108 : /* if a=b mod 4 set t=a-b, otherwise t=a+b, then strip powers of 2 */
109 : /* so that t <= (a+b)/4 < a/2 */
110 561306918 : gcd_plus_minus(a,b, t);
111 561306918 : if (is_pm1(t)) { set_avma(av); return int2n(v); }
112 561061062 : switch(abscmpii(t,b))
113 : {
114 364511679 : case -1: p1 = a; a = b; b = t; t = p1; break;
115 195271266 : case 1: swap(a,t); break;
116 1278117 : case 0: set_avma(av); b = shifti(b,v); return b;
117 : }
118 : }
119 : {
120 6686919 : long r[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
121 6686919 : r[2] = (long) gcduodd((ulong)b[2], (ulong)a[2]);
122 6686919 : set_avma(av); return shifti(r,v);
123 : }
124 : }
125 :
126 : GEN
127 619975284 : gcdii(GEN x, GEN y)
128 : {
129 619975284 : pari_sp av=avma;
130 628584693 : while (lgefint(y)-2 >= GCD_HALFGCD_LIMIT)
131 : {
132 8609409 : GEN M = HGCD0(x,y);
133 8609409 : x = gel(M,2); y = gel(M,3);
134 8609409 : if (signe(y) && expi(y)<magic_threshold(x))
135 : {
136 6543096 : swap(x,y);
137 6543096 : y = remii(y,x);
138 : }
139 8609409 : if (gc_needed(av, 1)) gerepileall(av,2,&x,&y);
140 : }
141 619975284 : return gerepileuptoint(av, gcdii_basecase(x,y));
142 : }
|
