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 18817005 : resiu(GEN y, ulong x)
19 : {
20 18817005 : long i, ly = lgefint(y);
21 18817005 : ulong xi = get_Fl_red(x);
22 : LOCAL_HIREMAINDER;
23 :
24 18817005 : hiremainder = 0;
25 96278583 : for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
26 18817005 : 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 567837438 : gcd_plus_minus(GEN x, GEN y, GEN res)
32 : {
33 567837438 : pari_sp av = avma;
34 567837438 : long lx = lgefint(x)-1;
35 567837438 : long ly = lgefint(y)-1, lt,m,i;
36 : GEN t;
37 :
38 567837438 : if ((x[lx]^y[ly]) & 3) /* x != y mod 4*/
39 284903844 : t = addiispec(x+2,y+2,lx-1,ly-1);
40 : else
41 282933594 : t = subiispec(x+2,y+2,lx-1,ly-1);
42 :
43 568564485 : lt = lgefint(t)-1; while (!t[lt]) lt--;
44 567837438 : m = vals(t[lt]); lt++;
45 567837438 : if (m == 0) /* 2^32 | t */
46 : {
47 27381 : for (i = 2; i < lt; i++) res[i] = t[i];
48 : }
49 567825417 : else if (t[2] >> m)
50 : {
51 535736229 : shift_right(res,t, 2,lt, 0,m);
52 : }
53 : else
54 : {
55 32089188 : lt--; t++;
56 32089188 : shift_right(res,t, 2,lt, t[1],m);
57 : }
58 567837438 : res[1] = evalsigne(1)|evallgefint(lt);
59 567837438 : set_avma(av);
60 567837438 : }
61 :
62 : /* uses modified right-shift binary algorithm now --GN 1998Jul23 */
63 : static GEN
64 623796450 : gcdii_basecase(GEN a, GEN b)
65 : {
66 : long v, w;
67 : pari_sp av;
68 : GEN t, p1;
69 :
70 623796450 : switch (abscmpii(a,b))
71 : {
72 177336339 : case 0: return absi(a);
73 214413882 : case -1: swap(a,b);
74 : }
75 446460111 : if (!signe(b)) return absi(a);
76 : /* here |a|>|b|>0. Try single precision first */
77 253340358 : if (lgefint(a)==3)
78 222281919 : return igcduu((ulong)a[2], (ulong)b[2]);
79 31058439 : if (lgefint(b)==3)
80 : {
81 18817005 : ulong u = resiu(a,(ulong)b[2]);
82 18817005 : if (!u) return absi(b);
83 11179989 : return igcduu((ulong)b[2], u);
84 : }
85 :
86 : /* larger than gcd: "set_avma(av)" gerepile (erasing t) is valid */
87 12241434 : av = avma; (void)new_chunk(lgefint(b)); /* HACK */
88 12241434 : t = remii(a,b);
89 12241434 : if (!signe(t)) { set_avma(av); return absi(b); }
90 :
91 8478459 : a = b; b = t;
92 8478459 : v = vali(a); a = shifti(a,-v); setabssign(a);
93 8478459 : w = vali(b); b = shifti(b,-w); setabssign(b);
94 8478459 : if (w < v) v = w;
95 8478459 : switch(abscmpii(a,b))
96 : {
97 168000 : case 0: set_avma(av); a = shifti(a,v); return a;
98 1660506 : case -1: swap(a,b);
99 : }
100 8310459 : 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 574570500 : 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 567837438 : gcd_plus_minus(a,b, t);
111 567837438 : if (is_pm1(t)) { set_avma(av); return int2n(v); }
112 567591570 : switch(abscmpii(t,b))
113 : {
114 368821446 : case -1: p1 = a; a = b; b = t; t = p1; break;
115 197499975 : case 1: swap(a,t); break;
116 1270149 : case 0: set_avma(av); b = shifti(b,v); return b;
117 : }
118 : }
119 : {
120 6733062 : long r[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
121 6733062 : r[2] = (long) gcduodd((ulong)b[2], (ulong)a[2]);
122 6733062 : set_avma(av); return shifti(r,v);
123 : }
124 : }
125 :
126 : GEN
127 623796450 : gcdii(GEN x, GEN y)
128 : {
129 623796450 : pari_sp av=avma;
130 632475654 : while (lgefint(y)-2 >= GCD_HALFGCD_LIMIT)
131 : {
132 8679204 : GEN M = HGCD0(x,y);
133 8679204 : x = gel(M,2); y = gel(M,3);
134 8679204 : if (signe(y) && expi(y)<magic_threshold(x))
135 : {
136 6572583 : swap(x,y);
137 6572583 : y = remii(y,x);
138 : }
139 8679204 : if (gc_needed(av, 1)) gerepileall(av,2,&x,&y);
140 : }
141 623796450 : return gerepileuptoint(av, gcdii_basecase(x,y));
142 : }
|
