## "Fossies" - the Fresh Open Source Software Archive

### Member "NZMATH-1.2.0/nzmath/arygcd.py" (19 Nov 2012, 4717 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "arygcd.py" see the Fossies "Dox" file reference documentation.

```    1 """
2 binary-like gcd algorithms for rational, gauss, and eisenstein integer
3 """
4
5 import math
6
7
8 def bit_num(a):
9     """
10     return the number of bits
11     """
12     if a == 0:
13         return 1
14     else:
15         return int(math.log(a, 2)) + 1
16
17
18 def binarygcd(a, b):
19     """
20     Return the greatest common divisor of 2 integers a and b
21     by binary gcd algorithm.
22     """
23     if a < b:
24         a, b = b, a
25     if b == 0:
26         return a
27     a, b = b, a % b
28     if b == 0:
29         return a
30     k = 0
31     while not a&1 and not b&1:
32         k += 1
33         a >>= 1
34         b >>= 1
35     while not a&1:
36         a >>= 1
37     while not b&1:
38         b >>= 1
39     while a:
40         while not a & 1:
41             a >>= 1
42         if abs(a) < abs(b):
43             a, b = b, a
44         a = (a - b) >> 1
45     return b << k
46
47
48 def arygcd_i(a1, a2, b1, b2):
49     """
50     Return the greatest common divisor of 2 gauss-integers a1+a2*i and b1+b2*i
51     by (1+i)-ary gcd algorithm.
52     """
53     if a1 == 0 and a2 == 0:
54        return b1, b2
55     elif b1 == 0 and b2 == 0:
56        return a1, a2
57     ap, bp = 0, 0
58     while (a1-a2)&1 == 0:
59         a1, a2 = (a1+a2) >> 1, (a2-a1) >> 1
60         ap += 1
61     while (b1-b2)&1 == 0:
62         b1, b2 = (b1+b2) >> 1, (b2-b1) >> 1
63         bp += 1
64     k = min(ap, bp)
65     while a1 != 0 or a2 != 0:
66         while (a1-a2) & 1 == 0:
67             a1, a2 = (a1+a2) >> 1, (a2-a1) >> 1
68         norm_a, norm_b = _ap_norm_i(a1, a2, b1, b2)
69         if norm_a < norm_b:
70             a1, a2, b1, b2 = b1, b2, a1, a2
71         a1, a2 = _FormAdj_i(a1, a2)
72         b1, b2 = _FormAdj_i(b1, b2)
73         a1, a2 = a1-b1, a2-b2
74     if 0 in (ap, bp):
75         return b1, b2
76     else:
77         s, t = _n_pow_i(1, 1, k)
78         return (b1*s)-(b2*t), (b1*t)+(b2*s)
79
80 def _ap_norm_i(a, b, c, d):
81     """
82     Return approximately norm of 2 gaussInteger a+b*i and c+d*i
83     """
84     a, b, c, d = abs(a), abs(b), abs(c), abs(d)
85     max_dig = max(bit_num(a), bit_num(b), bit_num(c), bit_num(d))
86     if max_dig > 6:
87         sft_num = max_dig - 6
88         a, b, c, d = a >> sft_num, b >> sft_num, c >> sft_num, d >> sft_num
89         return a*a + b*b, c*c + d*d
90     else:
91         return a*a + b*b, c*c + d*d
92
93 def _n_pow_i(a, b, n):
94     """
95     return (1+i)**k
96     """
97     x = a
98     y = b
99
100     for i in range(1, n):
101         x1 = (x*a) - (y*b)
102         y1 = (y*a) + (x*b)
103         x = x1
104         y = y1
105     return x, y
106
108     """
109     transform gaussInteger a+b*i ->  form 1+2(1+i)*(x+y*i)
110     """
111     if a & 1 == 0:
112         a, b = -b, a
113     if (b - a + 1) & 3 == 0:
114         return a, b
115     else:
116         return -a, -b
117
118
119 def arygcd_w(a1, a2, b1, b2):
120     """
121     Return the greatest common divisor of 2 eisensteinIntegers a1+a2*w and b1+b2*w
122     by (1-w)-ary gcd algorithm.
123     """
124     if a1 == 0 and a2 == 0:
125        return b1, b2
126     elif b1 == 0 and b2 == 0:
127        return a1, a2
128     ap, bp = 0, 0
129     while (a1 + a2) % 3 == 0:
130         a1, a2 = (a1 + a1 - a2) / 3, (a1 + a2) / 3
131         ap += 1
132     while (b1 + b2) % 3 == 0:
133         b1, b2 = (b1 + b1 - b2) / 3, (b1 + b2) / 3
134         bp += 1
135     k = min(ap, bp)
136     while a1 != 0 or a2 != 0:
137         while (a1 + a2) % 3 == 0:
138             a1, a2 = (a1 + a1 - a2) / 3, (a1 + a2) / 3
139         nrm_a, nrm_b = _ap_norm_w(a1, a2, b1, b2)
140         if nrm_a < nrm_b:
141             a1, a2, b1, b2 = b1, b2, a1, a2
142         a1, a2 = _FormAdj_w(a1, a2)
143         b1, b2 = _FormAdj_w(b1, b2)
144         a1, a2 = a1 - b1, a2 - b2
145     k1, k2 = _n_pow_w(k)
146     return k1*b1 - k2*b2, k1*b2 + k2*b1 - k2*b2
147
148 def _ap_norm_w(a, b, c, d):
149     """
150     Return approximately norm of 2 eisensteinInteger a+b*w and c+d*w
151     """
152     a, b, c, d = abs(a), abs(b), abs(c), abs(d)
153     max_dig = max(bit_num(a), bit_num(b), bit_num(c), bit_num(d))
154     if max_dig > 6:
155         sft_num = max_dig - 6
156         a, b, c, d = a >> sft_num, b >> sft_num, c >> sft_num, d >> sft_num
157         subst1, subst2 = (a - b) >> sft_num, (c - d) >> sft_num
158         return a*a + b*b + subst1*subst1, c*c + d*d + subst2*subst2
159     else:
160         return a*a + b*b + (a - b)*(a - b), c*c + d*d + (c -d)*(c -d)
161
162 def _n_pow_w(n):
163     """
164     return (1-w)**k
165     """
166     x, y = divmod(n, 2)
167     k1 = 3**x
168     if y == 1:
169         k2 = -k1
170     else:
171         k2 = 0
172     return k1, k2
173
175     """
176     transform eisensteinInteger a1+a2*w ->  form 1+3*(x+y*w)
177     """
178     if a1 % 3 == 0:
179         if a2 % 3 == -1 or a2 % 3 == 2:
180             return a1 - a2, a1
181         else:
182             return a2 - a1, -a1
183     elif a1 % 3 == 1:
184         if a2 % 3 == 1:
185             return a2, a2 -a1
186         else:
187             return a1, a2
188     else:
189         if a2 % 3 == -1 or a2 % 3 == 2:
190             return -a2, a1 -a2
191         else:
192             return -a1, -a2
```