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 19281390 : resiu(GEN y, ulong x)
19 : {
20 19281390 : long i, ly = lgefint(y);
21 19281390 : ulong xi = get_Fl_red(x);
22 : LOCAL_HIREMAINDER;
23 :
24 19281390 : hiremainder = 0;
25 97851456 : for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
26 19281390 : 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 588872964 : gcd_plus_minus(GEN x, GEN y, GEN res)
32 : {
33 588872964 : pari_sp av = avma;
34 588872964 : long lx = lgefint(x)-1;
35 588872964 : long ly = lgefint(y)-1, lt,m,i;
36 : GEN t;
37 :
38 588872964 : if ((x[lx]^y[ly]) & 3) /* x != y mod 4*/
39 295453581 : t = addiispec(x+2,y+2,lx-1,ly-1);
40 : else
41 293419383 : t = subiispec(x+2,y+2,lx-1,ly-1);
42 :
43 589635087 : lt = lgefint(t)-1; while (!t[lt]) lt--;
44 588872964 : m = vals(t[lt]); lt++;
45 588872964 : if (m == 0) /* 2^32 | t */
46 : {
47 28329 : for (i = 2; i < lt; i++) res[i] = t[i];
48 : }
49 588860457 : else if (t[2] >> m)
50 : {
51 555644493 : shift_right(res,t, 2,lt, 0,m);
52 : }
53 : else
54 : {
55 33215964 : lt--; t++;
56 33215964 : shift_right(res,t, 2,lt, t[1],m);
57 : }
58 588872964 : res[1] = evalsigne(1)|evallgefint(lt);
59 588872964 : set_avma(av);
60 588872964 : }
61 :
62 : /* uses modified right-shift binary algorithm now --GN 1998Jul23 */
63 : static GEN
64 640744311 : gcdii_basecase(GEN a, GEN b)
65 : {
66 : long v, w;
67 : pari_sp av;
68 : GEN t, p1;
69 :
70 640744311 : switch (abscmpii(a,b))
71 : {
72 177977445 : case 0: return absi(a);
73 221266956 : case -1: swap(a,b);
74 : }
75 462766866 : if (!signe(b)) return absi(a);
76 : /* here |a|>|b|>0. Try single precision first */
77 271918743 : if (lgefint(a)==3)
78 239953956 : return igcduu((ulong)a[2], (ulong)b[2]);
79 31964787 : if (lgefint(b)==3)
80 : {
81 19281390 : ulong u = resiu(a,(ulong)b[2]);
82 19281390 : if (!u) return absi(b);
83 11638125 : return igcduu((ulong)b[2], u);
84 : }
85 :
86 : /* larger than gcd: GC using set_avma(av) (erasing t) is valid */
87 12683397 : av = avma; (void)new_chunk(lgefint(b)); /* HACK */
88 12683397 : t = remii(a,b);
89 12683397 : if (!signe(t)) { set_avma(av); return absi(b); }
90 :
91 8841357 : a = b; b = t;
92 8841357 : v = vali(a); a = shifti(a,-v); setabssign(a);
93 8841357 : w = vali(b); b = shifti(b,-w); setabssign(b);
94 8841357 : if (w < v) v = w;
95 8841357 : switch(abscmpii(a,b))
96 : {
97 170136 : case 0: set_avma(av); a = shifti(a,v); return a;
98 1723092 : case -1: swap(a,b);
99 : }
100 8671221 : 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 595897719 : 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 588872964 : gcd_plus_minus(a,b, t);
111 588872964 : if (is_pm1(t)) { set_avma(av); return int2n(v); }
112 588594837 : switch(abscmpii(t,b))
113 : {
114 382715307 : case -1: p1 = a; a = b; b = t; t = p1; break;
115 204572748 : case 1: swap(a,t); break;
116 1306782 : case 0: set_avma(av); b = shifti(b,v); return b;
117 : }
118 : }
119 : {
120 7024755 : long r[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
121 7024755 : r[2] = (long) gcduodd((ulong)b[2], (ulong)a[2]);
122 7024755 : set_avma(av); return shifti(r,v);
123 : }
124 : }
125 :
126 : GEN
127 640744311 : gcdii(GEN x, GEN y)
128 : {
129 640744311 : pari_sp av=avma;
130 650114799 : while (lgefint(y)-2 >= GCD_HALFGCD_LIMIT)
131 : {
132 9370488 : GEN M = HGCD0(x,y);
133 9370488 : x = gel(M,2); y = gel(M,3);
134 9370488 : if (signe(y) && expi(y)<magic_threshold(x))
135 : {
136 7193490 : swap(x,y);
137 7193490 : y = remii(y,x);
138 : }
139 9370488 : if (gc_needed(av, 1)) (void)gc_all(av,2,&x,&y);
140 : }
141 640744311 : return gc_INT(av, gcdii_basecase(x,y));
142 : }
|
