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 18480990 : resiu(GEN y, ulong x)
19 : {
20 18480990 : long i, ly = lgefint(y);
21 18480990 : ulong xi = get_Fl_red(x);
22 : LOCAL_HIREMAINDER;
23 :
24 18480990 : hiremainder = 0;
25 143380728 : for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
26 18480990 : 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 558606939 : gcd_plus_minus(GEN x, GEN y, GEN res)
32 : {
33 558606939 : pari_sp av = avma;
34 558606939 : long lx = lgefint(x)-1;
35 558606939 : long ly = lgefint(y)-1, lt,m,i;
36 : GEN t;
37 :
38 558606939 : if ((x[lx]^y[ly]) & 3) /* x != y mod 4*/
39 280283814 : t = addiispec(x+2,y+2,lx-1,ly-1);
40 : else
41 278323125 : t = subiispec(x+2,y+2,lx-1,ly-1);
42 :
43 559334046 : lt = lgefint(t)-1; while (!t[lt]) lt--;
44 558606939 : m = vals(t[lt]); lt++;
45 558606939 : if (m == 0) /* 2^32 | t */
46 : {
47 27321 : for (i = 2; i < lt; i++) res[i] = t[i];
48 : }
49 558594933 : else if (t[2] >> m)
50 : {
51 527018259 : shift_right(res,t, 2,lt, 0,m);
52 : }
53 : else
54 : {
55 31576674 : lt--; t++;
56 31576674 : shift_right(res,t, 2,lt, t[1],m);
57 : }
58 558606939 : res[1] = evalsigne(1)|evallgefint(lt);
59 558606939 : set_avma(av);
60 558606939 : }
61 :
62 : /* uses modified right-shift binary algorithm now --GN 1998Jul23 */
63 : static GEN
64 611053446 : gcdii_basecase(GEN a, GEN b)
65 : {
66 : long v, w;
67 : pari_sp av;
68 : GEN t, p1;
69 :
70 611053446 : switch (abscmpii(a,b))
71 : {
72 175817895 : case 0: return absi(a);
73 206359299 : case -1: swap(a,b);
74 : }
75 435235551 : if (!signe(b)) return absi(a);
76 : /* here |a|>|b|>0. Try single precision first */
77 243439374 : if (lgefint(a)==3)
78 212914146 : return igcduu((ulong)a[2], (ulong)b[2]);
79 30525228 : if (lgefint(b)==3)
80 : {
81 18480990 : ulong u = resiu(a,(ulong)b[2]);
82 18480990 : if (!u) return absi(b);
83 11121342 : return igcduu((ulong)b[2], u);
84 : }
85 :
86 : /* larger than gcd: "set_avma(av)" gerepile (erasing t) is valid */
87 12044238 : av = avma; (void)new_chunk(lgefint(b)); /* HACK */
88 12044238 : t = remii(a,b);
89 12044238 : if (!signe(t)) { set_avma(av); return absi(b); }
90 :
91 8393322 : a = b; b = t;
92 8393322 : v = vali(a); a = shifti(a,-v); setabssign(a);
93 8393322 : w = vali(b); b = shifti(b,-w); setabssign(b);
94 8393322 : if (w < v) v = w;
95 8393322 : switch(abscmpii(a,b))
96 : {
97 175224 : case 0: set_avma(av); a = shifti(a,v); return a;
98 1644312 : case -1: swap(a,b);
99 : }
100 8218098 : 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 565246500 : 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 558606939 : gcd_plus_minus(a,b, t);
111 558606939 : if (is_pm1(t)) { set_avma(av); return int2n(v); }
112 558361098 : switch(abscmpii(t,b))
113 : {
114 362719770 : case -1: p1 = a; a = b; b = t; t = p1; break;
115 194370816 : case 1: swap(a,t); break;
116 1270512 : case 0: set_avma(av); b = shifti(b,v); return b;
117 : }
118 : }
119 : {
120 6639561 : long r[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
121 6639561 : r[2] = (long) gcduodd((ulong)b[2], (ulong)a[2]);
122 6639561 : set_avma(av); return shifti(r,v);
123 : }
124 : }
125 :
126 : GEN
127 611053446 : gcdii(GEN x, GEN y)
128 : {
129 611053446 : pari_sp av=avma;
130 619781478 : while (lgefint(y)-2 >= GCD_HALFGCD_LIMIT)
131 : {
132 8728032 : GEN M = HGCD0(x,y);
133 8728032 : x = gel(M,2); y = gel(M,3);
134 8728032 : if (signe(y) && expi(y)<magic_threshold(x))
135 : {
136 6703932 : swap(x,y);
137 6703932 : y = remii(y,x);
138 : }
139 8728032 : if (gc_needed(av, 1)) gerepileall(av,2,&x,&y);
140 : }
141 611053446 : return gerepileuptoint(av, gcdii_basecase(x,y));
142 : }
|
